***********************************************************************
*     表題 u-v-time-graph-4-3.f                                       *
*                                                                     *
*     k≠0.0 の慣性重力波の速度ベクトルとその成分, 圧力傾度力の       * 
*     一周期分時間変化の図示                                          *
*                                                                     *
*     u-v-time-graph-4.f, u-v-time-graph-4-2.f                        *
*     と異なる緯度のグラフ図示用                                      *
*                                                                     *
*     履歴 2002/01/28 H.Taniguchi                                     *
*          2002/02/12 H.Taniguchi                                     *
*          2002/02/19 H.Taniguchi u-v-time-graph-4-3.f 版作成         *
*     All rights reserved, Copyright (C) 2002, Hiroshi TANIGUCHI      *
***********************************************************************

      integer n
      integer xgr
      integer ygr
      integer tgr
      real k
      real pi                   !" 円周率
      real dx                   !" x グリッド幅(度)
      real dxlon                !" 一波長あたりの経度

      parameter (xgr=100,ygr=100,tgr=100,numv=34) !" グリッド数の設定
      parameter (ymin=-2.0, ymax=2.0)
      parameter (pi=3.1415927)

      real H(0:2,ygr)           !" エルミート多項式
      real Hdy(0:2,ygr)         !" エルミート多項式微分形

      real a                    !" 地球の半径 
      real dg                   !" 単位グリッド
      real t                    !" 時間
      real gamma                !" 南北シアー (du/dy)
      real beta                 !" df/dy
      real ys                   !" gamma/2*beta
      real yc                   !" 
      real dy                   !" 格子点間隔(m 単位)
      real dydg                 !" 格子点間隔(緯度単位)
      real y(1:ygr)             !" y 座標(u, phi)
      real y2(1:ygr)            !" y 座標(u, phi)
      real vy(0:ygr)            !" y 座標(v) 
      real vy2(0:ygr)           !" y 座標(v) 
      real x(1:xgr)             !" x 座標(ラジアン)
      real xx(1:xgr)            !" x 座標(経度)
      real xxx(1:xgr)           !" x 座標(経度)  ※ 無次元

      real t1(xgr), u1(xgr), v1(xgr), px1(xgr), py1(xgr)
      real dudt1(xgr), dvdt1(xgr), betayu1(xgr), betayv1(xgr)
c      real checkx1(xgr), checky1(xgr)
      real t2(xgr), u2(xgr), v2(xgr), px2(xgr), py2(xgr)
      real dudt2(xgr), dvdt2(xgr), betayu2(xgr), betayv2(xgr)
c      real checkx2(xgr), checky2(xgr)
      real t3(xgr), u3(xgr), v3(xgr), px3(xgr), py3(xgr)
      real dudt3(xgr), dvdt3(xgr), betayu3(xgr), betayv3(xgr)
c      real checkx3(xgr), checky3(xgr)
      real uve1(1:numv, 3)
      real vve1(1:numv, 3)
      real uve2(1:numv, 3)
      real vve2(1:numv, 3)
      real uve3(1:numv, 3)
      real vve3(1:numv, 3)
      real uve4(1:numv, 3)
      real vve4(1:numv, 3)
      real vecty(3)

      real urvect(1:ygr)        !" 固有ベクトル(u の実部) 
      real uivect(1:ygr)        !" 固有ベクトル(u の虚部) 
      real uve(xgr,ygr)         !" ベクトル表示用変数( u )

      real vrvect(1:ygr)        !" 固有ベクトル(v の実部)
      real vivect(1:ygr)        !" 固有ベクトル(v の虚部)
      real vve(xgr,ygr)         !" ベクトル表示用変数( v )

      real prvect(1:ygr)        !" 固有ベクトル(phi の実部)
      real pivect(1:ygr)        !" 固有ベクトル(phi の虚部)

      real omega
      real rval2                !" 固有値実数部分(固有ベクトル用)
      real ival2                !" 固有値虚数部分(固有ベクトル用)

      real ppr(xgr,ygr,tgr)     !" ジオポテンシャル(phi) ※一波長分データ作成用
      real uur(xgr,ygr,tgr)     !" 東西速度(u)           ※一波長分データ作成用
      real vvr(xgr,ygr,tgr)     !" 南北速度(v)           ※一波長分データ作成用

      real dudt(xgr,ygr,tgr)     
      real betayv(xgr,ygr,tgr)   
      real dpdx(xgr,ygr,tgr)      
      real dvdt(xgr,ygr,tgr)     
      real betayu(xgr,ygr,tgr)   
      real dpdy(xgr,ygr,tgr)     
      real dpdt(xgr,ygr,tgr)     
      real ghdudx(xgr,ygr,tgr)   
      real ghdvdy(xgr,ygr,tgr)   

      character csgi*1, aaa

*---- !" モード番号, 東西波数, 波番号入力

      write(*,*) 'n=?' 
      read(*,*) n
      write(*,*) 'k=?'
      read(*,*) k

*---- !" 東西波数の有次元化

      write(*,*) '有次元 k =',k

      if (k.ne.0.) then 
         dxlon = 360./k         !" 一波長あたりの経度
         write(*,*) '一波長あたりの有次元経度 dxlon = 360./k =', dxlon
      else
         dxlon = 360. 
         write(*,*) '一波長あたりの有次元経度 dxlon =', dxlon
      end if

      write(*,*) '何波長分の絵を描きますか?'
      read(*,*) kgrid

*---- x 座標データ

      if (k.ne.0.) then 
         dx = (dxlon*kgrid)/(xgr-1) !" kgrid 波長分描いた時の1grid あたり経度
      else
         dx = dxlon/(xgr-1)
      end if

      write(*,*) kgrid,'波長描かせた時の 1 グリッドあたりの経度 dx=',dx
      
      do 1230 ix=1,xgr
         xx(ix) = dx*(ix-1)               !" ix 番目のグリッドの経度
         write(*,*) 'xx(',ix,')=',xx(ix)
         x(ix)  = (pi*xx(ix))/180.        !" ix 番目のグリッドのラジアン 
         xxx(ix) = x(ix)                  !" 有次元 x 座標 (ラジアン)
         write(*,*) 'xxx(',ix,')=',xxx(ix)
 1230 continue

      write(*,*) '(1)東進慣性重力波, (2)西進慣性重力波' 
      read(*,*) nwave

*---- !" 振動数計算

      if (nwave.eq.1) omega =  sqrt(k**2 + 2.*n + 1.)
      if (nwave.eq.2) omega = -sqrt(k**2 + 2.*n + 1.)

      write(*,*) '振動数=', omega
      write(*,*) '周期=', 2.*pi/omega
      read(*,*) aaa

      rval2 = omega
      ival2 = 0.

*---- y 座標データ

      a = 6378.e3
      dg=2.D0*pi*a/360.D0                      !" 1度あたりの距離(m)

      gamma = 0.00001
      beta  = 2.29e-11

      ys=gamma/(2.D0*beta)                     !" 不安定領域中心位置
      write(*,*) 'ys=',ys

      dy=3.D0*ys/dble(ygr)
      dydg=dy/dg

      do ii=0,ygr

         vy(ii) =(-1.5D0*ys+dy*(ii))/dg        !" v の y座標 (緯度)
         vy2(ii)=(-1.5*ys+dy*(ii))             !" v の y座標

      end do

      do iii=1,ygr

         y(iii)=(vy(iii)+vy(iii-1))/2.D0       !" u,phi の y座標 (緯度)
         y2(iii)=(vy2(iii)+vy2(iii-1))/2.      !" u,phi の y座標

         write(*,*) 'y(',iii,')=',y(iii)

*---- !" エルミート多項式 H と微分形 Hdy 計算

         H(0,iii) = 1.
         Hdy(0,iii) = 0.
         H(1,iii) = 2.*vy(iii)
c         H(1,iii) = 2.*y(iii)
         Hdy(1,iii) = 2.
         H(2,iii) = 4.*vy(iii)**2 - 2. 
c         H(2,iii) = 4.*y(iii)**2 - 2. 
         Hdy(2,iii) = 8.*y(iii)

      end do

      do inny=1,ygr    

        vrvect(inny) = H(n,inny)*exp( -(1./2.)*(y(inny)**2) )
        vivect(inny) = 0.

        urvect(inny) = 0.
        uivect(inny) 
     & =  (1./(omega**2 - k**2))*( omega*vy(inny)*vrvect(inny) 
c     & =  (1./(omega**2 - k**2))*( omega*y(inny)*vrvect(inny) 
     &   - k*( 
     &         Hdy(n,inny)*exp( -(1./2.)*(vy(inny)**2) ) 
c     &         Hdy(n,inny)*exp( -(1./2.)*(y(inny)**2) ) 
     &       - vy(inny)*vrvect(inny) 
c     &       - y(inny)*vrvect(inny) 
     &        ) 
     &                            )

        prvect(inny) = 0.
        pivect(inny) 
     & =  (1./(omega**2 - k**2))*( k*vy(inny)*vrvect(inny) 
c     & =  (1./(omega**2 - k**2))*( k*y(inny)*vrvect(inny) 
     &   - omega*(
     &              Hdy(n,inny)*exp( -(1./2.)*(vy(inny)**2) ) 
c     &              Hdy(n,inny)*exp( -(1./2.)*(y(inny)**2) ) 
     &            - vy(inny)*vrvect(inny)
c     &            - y(inny)*vrvect(inny)
     &            )      
     &                            )

      end do

      do it = 1, xgr
         t = ((2.*pi/omega)/(xgr-1))*(it-1)
        do iix = 1,xgr
           do iiy = 1,ygr 
              ppr(iix,iiy,it) = ( prvect(iiy)*cos(k*xxx(iix)-rval2*t)
     &                          - pivect(iiy)*sin(k*xxx(iix)-rval2*t) )

              uur(iix,iiy,it) = ( urvect(iiy)*cos(k*xxx(iix)-rval2*t)
     &                          - uivect(iiy)*sin(k*xxx(iix)-rval2*t) )

              vvr(iix,iiy,it) = ( vrvect(iiy)*cos(k*xxx(iix)-rval2*t)
     &                          - vivect(iiy)*sin(k*xxx(iix)-rval2*t) )
           end do
        end do
      end do

*---- 各項を計算(位相関係を明らかにするため)

*---- 欠損値処理

      do iit=1,xgr
         do inx=1,xgr          

            dudt(inx,1,iit)     = RMISS
            dudt(inx,ygr,iit)   = RMISS
            betayv(inx,1,iit)   = RMISS
            betayv(inx,ygr,iit) = RMISS
            dpdx(inx,1,iit)     = RMISS 
            dpdx(inx,ygr,iit)   = RMISS

            dvdt(inx,1,iit)     = RMISS
            dvdt(inx,ygr,iit)   = RMISS
            betayu(inx,1,iit)   = RMISS
            betayu(inx,ygr,iit) = RMISS
            dpdy(inx,1,iit)     = RMISS
            dpdy(inx,ygr,iit)   = RMISS

            dpdt(inx,1,iit)     = RMISS
            dpdt(inx,ygr,iit)   = RMISS
            ghdudx(inx,1,iit)   = RMISS
            ghdudx(inx,ygr,iit) = RMISS
            ghdvdy(inx,1,iit)   = RMISS
            ghdvdy(inx,ygr,iit) = RMISS

         end do
      end do

      do itnx=1,xgr 

         t = ((2.*pi/omega)/(xgr-1))*(itnx-1)
         
         do inx=1,xgr
            do iny=2,ygr-1 

*-------- u の式 -------------------------------------------------------

*---- ∂u/∂t 

      dudt(inx,iny,itnx) 
     &           = (
     &               ( rval2*uivect(iny) + ival2*urvect(iny) )*( 
     &                 cos( k*xxx(inx) - rval2*t )              ) 
     &              +( rval2*urvect(iny) - ival2*uivect(iny) )*( 
     &                 sin( k*xxx(inx) - rval2*t )              ) 
     &              )*exp( ival2*t )

*---- βyv

      betayv(inx,iny,itnx) 
     &           = y(iny)*(
     &                           vrvect(iny)*cos( k*xxx(inx) - rval2*t )
     &                          -vivect(iny)*sin( k*xxx(inx) - rval2*t )
     &                          )*exp( ival2*t )

*---- -∂Φ/∂x

      dpdx(inx,iny,itnx)
     &           = k*(
     &                 prvect(iny)*sin( k*xxx(inx) - rval2*t )
     &                +pivect(iny)*cos( k*xxx(inx) - rval2*t ) 
     &                )*exp( ival2*t )

*-------- v の式 -------------------------------------------------------

*---- ∂v/∂t 

      dvdt(inx,iny,itnx) 
     &           = (
     &               ( rval2*vivect(iny) + ival2*vrvect(iny) )*( 
     &                 cos( k*xxx(inx) - rval2*t )              ) 
     &              +( rval2*vrvect(iny) - ival2*vivect(iny) )*( 
     &                 sin( k*xxx(inx) - rval2*t )              ) 
     &              )*exp( ival2*t )

*---- -βyu

      betayu(inx,iny,itnx) 
     &           = -vy(iny)*(
c     &           = -y(iny)*(
     &                       urvect(iny)*cos( k*xxx(inx) - rval2*t )
     &                      -uivect(iny)*sin( k*xxx(inx) - rval2*t )
     &                      )*exp( ival2*t )

*---- -∂Φ/∂y

      dpdy(inx,iny,itnx)
     &           = -(
     &                ( ( prvect(iny+1) - prvect(iny-1) )/(2.*dydg) )*
     &                    cos( k*xxx(inx) - rval2*t )
     &               -( ( pivect(iny+1) - pivect(iny-1) )/(2.*dydg) )* 
     &                    sin( k*xxx(inx) - rval2*t )
     &               )*exp( ival2*t )

*-------- Φ の式 ------------------------------------------------------

*---- ∂Φ/∂t 

      dpdt(inx,iny,itnx) 
     &           = (
     &               ( rval2*pivect(iny) + ival2*prvect(iny) )*( 
     &                 cos( k*xxx(inx) - rval2*t )              ) 
     &              +( rval2*prvect(iny) - ival2*pivect(iny) )*( 
     &                 sin( k*xxx(inx) - rval2*t )              ) 
     &              )*exp( ival2*t )

*---- -gH・∂u/∂x (ただし, gH = 1)

      ghdudx(inx,iny,itnx) 
     &           = k*(
     &                 urvect(iny)*sin( k*xxx(inx) - rval2*t )
     &                +uivect(iny)*cos( k*xxx(inx) - rval2*t )
     &                )*exp( ival2*t )

*---- -gH・∂v/∂y (ただし, gH = 1)

      ghdvdy(inx,iny,itnx)
     &           = -( 
     &               ( ( vrvect(iny+1) - vrvect(iny-1) )/(2.*dydg) )*
     &                    cos( k*xxx(inx) - rval2*t )
     &              -( ( vivect(iny+1) - vivect(iny-1) )/(2.*dydg) )* 
     &                    sin( k*xxx(inx) - rval2*t )
     &               )*exp( ival2*t )

      end do
      end do
      end do

      vecty(1) = -0.1
      vecty(2) = 0.
      vecty(3) = 0.1

c         write(*,*) '何番目の y 座標を用いますか'
c         read iy

         ii1x = 1
         ii1y = 25
c         ii1y = 66

      do i1 = 1, xgr
c      do i1 = 0, xgr
         t1(i1)      = ((2.*pi/omega)/(xgr-1))*(i1-1)
         u1(i1)      = uur(ii1x,ii1y,i1)
         v1(i1)      = vvr(ii1x,ii1y,i1)
         dudt1(i1)   = dudt(ii1x,ii1y,i1)
         dvdt1(i1)   = dvdt(ii1x,ii1y,i1)
         betayv1(i1) = betayv(ii1x,ii1y,i1) 
         betayu1(i1) = betayu(ii1x,ii1y,i1) 
         px1(i1)     = dpdx(ii1x,ii1y,i1)
         py1(i1)     = dpdy(ii1x,ii1y,i1)
c         checkx1(i1) = dudt1(i1)-betayv1(i1)-px1(i1)
c         checky1(i1) = dvdt1(i1)-betayu1(i1)-py1(i1)
      end do

      do ix1 = 1, numv
         uve1(ix1, 1) = 0.
         vve1(ix1, 1) = 0.
         uve1(ix1, 2) = u1(1+int((xgr-1)/(numv-1))*(ix1-1))
         vve1(ix1, 2) = v1(1+int((xgr-1)/(numv-1))*(ix1-1))
         uve1(ix1, 3) = 0.
         vve1(ix1, 3) = 0.
      end do

c         write(*,*) '(2) 何番目の y 座標を用いますか?'
c         read ii2y

         ii2x = 1 
         ii2y = 42
c         ii2y = 79

      do i2 = 1, xgr
c      do i2 = 0, xgr
         t2(i2)      = ((2.*pi/omega)/(xgr-1))*(i2-1)
         u2(i2)      = uur(ii2x,ii2y,i2)
         v2(i2)      = vvr(ii2x,ii2y,i2)
         dudt2(i2)   = dudt(ii2x,ii2y,i2)
         dvdt2(i2)   = dvdt(ii2x,ii2y,i2)
         betayv2(i2) = betayv(ii2x,ii2y,i2) 
         betayu2(i2) = betayu(ii2x,ii2y,i2) 
         px2(i2)     = dpdx(ii2x,ii2y,i2)
         py2(i2)     = dpdy(ii2x,ii2y,i2)
c         checkx2(i2) = dudt2(i2)-betayv2(i2)-px2(i2)
c         checky2(i2) = dvdt2(i2)-betayu2(i2)-py2(i2)
      end do

      do ix2 = 1, numv
         uve2(ix2, 1) = 0.
         vve2(ix2, 1) = 0.
         uve2(ix2, 2) = u2(1+int((xgr-1)/(numv-1))*(ix2-1))
         vve2(ix2, 2) = v2(1+int((xgr-1)/(numv-1))*(ix2-1))
         uve2(ix2, 3) = 0.
         vve2(ix2, 3) = 0.
      end do

c         write(*,*) '(3) 何番目の y 座標を用いますか?'
c         read ii3y

         ii3x = 1
         ii3y = 51
c         ii3y = 85

      do i3 = 1, xgr
c      do i3 = 0, xgr
         t3(i3)      = ((2.*pi/omega)/(xgr-1))*(i3-1)
         u3(i3)      = uur(ii3x,ii3y,i3)
         v3(i3)      = vvr(ii3x,ii3y,i3)
         dudt3(i3)   = dudt(ii3x,ii3y,i3)
         dvdt3(i3)   = dvdt(ii3x,ii3y,i3)
         betayv3(i3) = betayv(ii3x,ii3y,i3) 
         betayu3(i3) = betayu(ii3x,ii3y,i3) 
         px3(i3)     = dpdx(ii3x,ii3y,i3)
         py3(i3)     = dpdy(ii3x,ii3y,i3)
c         checkx3(i3) = dudt3(i3)-betayv3(i3)-px3(i3)
c         checky3(i3) = dvdt3(i3)-betayu3(i3)-py3(i3)
       end do

      do ix3 = 1, numv
         uve3(ix3, 1) = 0.
         vve3(ix3, 1) = 0.
         uve3(ix3, 2) = u3(1+int((xgr-1)/(numv-1))*(ix3-1))
         vve3(ix3, 2) = v3(1+int((xgr-1)/(numv-1))*(ix3-1))
         uve3(ix3, 3) = 0.
         vve3(ix3, 3) = 0.
      end do


      write(*,*) 't1(1)=',t1(1)
      write(*,*) 't2(1)=',t2(1)
      write(*,*) 't3(1)=',t3(1)
      write(*,*) 't1(xgr)=',t1(xgr)
      write(*,*) 't2(xgr)=',t2(xgr)
      write(*,*) 't3(xgr)=',t3(xgr)

      write(*,*) ' workstation id (i) ? ;'
      call sgpwsn
      read(*,*) iws

*--------- 出力装置オープン      
      call sgopn(-abs(iws))

*--------- コーナマーカを消す
      call sglset('lcorner', .FALSE.)

*--------- 第一レベルのマージン
      call slmgn( 0., 0., 0., 0. )

*--------- 画面分割 (縦'T' 横'Y', 縦の分割数, 横の分割数)
      call sldiv( 'S', 1, 1 )
      call slmgn( 0., 0., 0., 0. )
      call sgsvpt( 0., 0., 0., 0. )

*--------- 1フレーム目の設定
      call grfrm
      call glrget('RUNDEF',RUNDEF)
      call grswnd( t3(1), t3(xgr), ymin, ymax )
c      call grswnd( RUNDEF, RUNDEF, ymin, ymax )
      call grsvpt( 0.15, 0.90, 0.70, 0.90 )
*--------- 書きたいデータの中から最大値と最小値をみつける
      call usspnt(xgr,t3,u3)
      call usspnt(xgr,t3,v3)
      call usspnt(xgr,t3,dudt3)
      call usspnt(xgr,t3,dvdt3)
      call usspnt(xgr,t3,betayu3)
      call usspnt(xgr,t3,betayv3)
      call usspnt(xgr,t3,py3)
      call usspnt(xgr,t3,px3)
c      call usspnt(xgr,t3,checkx3)
c      call usspnt(xgr,t3,checky3)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 上下の座標軸を描く
      call usxaxs('B')
      call usxaxs('T')
*--------- 左右の座標軸を描く
      call usyaxs('l')
      call usyaxs('r')

*--------- X軸のタイトル指定
      call uzfact( 0.7 )
*--------- Y軸のタイトル指定
      call uysttl('l','Amp.',0.0)
      call uzfact( 1/0.7 )

*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .TRUE. )
*--------- 折れ線のラベルの設定
      call sgsplc( 'u' )
*--------- 折れ線を描く
      call uulinz(xgr,t3,u3,1,3)
*--------- 折れ線のラベルの設定
      call sgsplc( 'v' )
      call uulinz(xgr,t3,v3,3,3)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'dudt' )
*--------- 折れ線を描く
      call uulinz(xgr,t3,dudt3,1,43)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.2) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'dvdt' )
      call uulinz(xgr,t3,dvdt3,3,43)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.9) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'yu' )
*--------- 折れ線を描く
      call uulinz(xgr,t3,betayu3,3,33)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'yv' )
      call uulinz(xgr,t3,betayv3,1,33)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.7) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'py' )
      call uulinz(xgr,t3,py3,3,23)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'px' )
      call uulinz(xgr,t3,px3,1,23)
*--------- 折れ線のラベルの設定
c      call sgrset( 'FFCT', 0.1) !" 線分の始まり部分(全長に対して)
c      call sgsplc( 'checkx' )
c      call uulinz(xgr,t3,checkx3,3,73)
*--------- 折れ線のラベルの設定
c      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
c      call sgsplc( 'checky' )
c      call uulinz(xgr,t3,checky3,1,73)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .FALSE. )

*--------- 1.5フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grsvpt( 0.15, 0.90, 0.86, 0.98 )
      call grswnd( t3(1), t3(xgr), vecty(1), vecty(3))
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call uwsgxb( t1(1), t1(xgr), numv )
      call uwsgya( vecty, 3 )

*--------- ベクトル図を描くルーチン
      call uglset( 'LMSG', .FALSE.)
      call uglset( 'LUMSG', .FALSE.)
      call ugiset( 'ICENT', 0) !" ベクトルの始点の位置設定
      call ugvect( uve3, numv, vve3, numv, numv, 3 )

*--------- 2フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grswnd( t2(1), t2(xgr), ymin, ymax )
c      call grswnd( RUNDEF, RUNDEF, ymin, ymax )
      call grsvpt( 0.15, 0.90, 0.40, 0.60 )
*--------- 書きたいデータの中から最大値と最小値をみつける
      call usspnt(xgr,t2,u2)
      call usspnt(xgr,t2,v2)
      call usspnt(xgr,t2,dudt2)
      call usspnt(xgr,t2,dvdt2)
      call usspnt(xgr,t2,betayu2)
      call usspnt(xgr,t2,betayv2)
      call usspnt(xgr,t2,py2)
      call usspnt(xgr,t2,px2)
c      call usspnt(xgr,t2,checkx2)
c      call usspnt(xgr,t2,checky2)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 上下の座標軸を描く
      call usxaxs('B')
      call usxaxs('T')
*--------- 左右の座標軸を描く
      call usyaxs('l')
      call usyaxs('r')

*--------- X軸のタイトル指定
      call uzfact( 0.7 )
c      call uxpttl('t',1,'LAT '//csgi(210)//'3',0.0)
c      call uxsttl('b','Time',0.0)
*--------- Y軸のタイトル指定
      call uysttl('l','Amp.',0.0)
      call uzfact( 1/0.7 )

*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .TRUE. )
*--------- 折れ線のラベルの設定
      call sgsplc( 'u' )
*--------- 折れ線を描く
      call uulinz(xgr,t2,u2,1,3)
*--------- 折れ線のラベルの設定
      call sgsplc( 'v' )
      call uulinz(xgr,t2,v2,3,3)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'dudt' )
*--------- 折れ線を描く
      call uulinz(xgr,t2,dudt2,1,43)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.2) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'dvdt' )
      call uulinz(xgr,t2,dvdt2,3,43)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.9) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'yu' )
*--------- 折れ線を描く
      call uulinz(xgr,t2,betayu2,3,33)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'yv' )
      call uulinz(xgr,t2,betayv2,1,33)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.7) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'py' )
      call uulinz(xgr,t2,py2,3,23)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'px' )
      call uulinz(xgr,t2,px2,1,23)
*--------- 折れ線のラベルの設定
c      call sgrset( 'FFCT', 0.1) !" 線分の始まり部分(全長に対して)
c      call sgsplc( 'checkx' )
c      call uulinz(xgr,t2,checkx2,3,73)
*--------- 折れ線のラベルの設定
c      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
c      call sgsplc( 'checky' )
c      call uulinz(xgr,t2,checky2,1,73)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .FALSE. )

*--------- 2.5フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grsvpt( 0.15, 0.90, 0.57, 0.66 )
      call grswnd( t2(1), t2(xgr), vecty(1), vecty(3))
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call uwsgxb( t1(1), t1(xgr), numv )
      call uwsgya( vecty, 3 )

*--------- ベクトル図を描くルーチン
      call uglset( 'LMSG', .FALSE.)
      call uglset( 'LUMSG', .FALSE.)
      call ugiset( 'ICENT', 0) !" ベクトルの始点の位置設定
      call ugvect( uve2, numv, vve2, numv, numv, 3 )

*--------- 3フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grswnd( t1(1), t1(xgr), ymin, ymax )
c      call grswnd( RUNDEF, RUNDEF, ymin, ymax )
      call grsvpt( 0.15, 0.90, 0.10, 0.30 )
*--------- 書きたいデータの中から最大値と最小値をみつける
      call usspnt(xgr,t1,u1)
      call usspnt(xgr,t1,v1)
      call usspnt(xgr,t1,dudt1)
      call usspnt(xgr,t1,dvdt1)
      call usspnt(xgr,t1,betayu1)
      call usspnt(xgr,t1,betayv1)
      call usspnt(xgr,t1,py1)
      call usspnt(xgr,t1,px1)
c      call usspnt(xgr,t1,checkx1)
c      call usspnt(xgr,t1,checky1)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 上下の座標軸を描く
      call usxaxs('B')
      call usxaxs('T')
*--------- 左右の座標軸を描く
      call usyaxs('l')
      call usyaxs('r')

*--------- X軸のタイトル指定
*--------- トップのタイトル指定
      call uzfact( 0.7 )
c      call uxpttl('t',1,'LAT 1',0.0)
      call uxsttl('b','Time',0.0)
*--------- Y軸のタイトル指定
      call uysttl('l','Amp.',0.0)
      call uzfact( 1/0.7 )

*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .TRUE. )
*--------- 折れ線のラベルの設定
      call sgsplc( 'u' )
*--------- 折れ線を描く
      call uulinz(xgr,t1,u1,1,3)
*--------- 折れ線のラベルの設定
      call sgsplc( 'v' )
      call uulinz(xgr,t1,v1,3,3)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'dudt' )
*--------- 折れ線を描く
      call uulinz(xgr,t1,dudt1,1,43)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.2) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'dvdt' )
      call uulinz(xgr,t1,dvdt1,3,43)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.9) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'yu' )
*--------- 折れ線を描く
      call uulinz(xgr,t1,betayu1,3,33)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'yv' )
      call uulinz(xgr,t1,betayv1,1,33)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.7) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'py' )
      call uulinz(xgr,t1,py1,3,23)
*--------- 折れ線のラベルの設定
      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
      call sgsplc( 'px' )
      call uulinz(xgr,t1,px1,1,23)
*--------- 折れ線のラベルの設定
c      call sgrset( 'FFCT', 0.1) !" 線分の始まり部分(全長に対して)
c      call sgsplc( 'checkx' )
c      call uulinz(xgr,t1,checkx1,3,73)
*--------- 折れ線のラベルの設定
c      call sgrset( 'FFCT', 0.4) !" 線分の始まり部分(全長に対して)
c      call sgsplc( 'checky' )
c      call uulinz(xgr,t1,checky1,1,73)
*--------- ラベル文字列の許可
      call sglset( 'LCHAR' , .FALSE. )

*--------- 3.5フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
      call grsvpt( 0.15, 0.90, 0.27, 0.38 )
      call grswnd( t1(1), t1(xgr), vecty(1), vecty(3))
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call uwsgxb( t1(1), t1(xgr), numv )
      call uwsgya( vecty, 3 )

*--------- ベクトル図を描くルーチン
      call uglset( 'LMSG', .FALSE.)
      call uglset( 'LUMSG', .FALSE.)
      call ugiset( 'ICENT', 0) !" ベクトルの始点の位置設定
      call ugvect( uve1, numv, vve1, numv, numv, 3 )

*--------- ベクトルが 45 度の方向を向く時間
c      call sglnzv( 0.29, 0.1, 0.22, 0.9, 1) 
c      call sglnzv( 0.65, 0.1, 0.58, 0.9, 1) 
      call sgtxzv( 0.73, 0.13, 'LAT -1.5',0.02,0,-1,3) 
c      call sgtxzv( 0.75, 0.13, 'LAT 1.5',0.02,0,-1,3) 
      call sgtxzv( 0.73, 0.43, 'LAT -0.5',0.02,0,-1,3) 
c      call sgtxzv( 0.73, 0.43, 'LAT 2.0',0.02,0,-1,3) 
      call sgtxzv( 0.73, 0.73, 'LAT 0.0',0.02,0,-1,3) 
c      call sgtxzv( 0.73, 0.73, 'LAT 2.5',0.02,0,-1,3) 

*---------出力装置のクローズ
      call grcls
      end
