*表題  赤道波の水平構造の理論解と各項の位相の水平構造を表示するプログラム
*
*履歴  2001/04/15 H.Taniguchi
*      2001/06/25 H.Taniguchi (指定した波長分だけ描く, コンター間隔均等化)

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

      parameter (xgr=61,ygr=61) !" グリッド数の設定
      parameter (nvxg=16, nvyg=16,pi=3.1415927)

      real XMIN, XMAX

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

      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 omega
      real rval2                !" 固有値実数部分(固有ベクトル用)
      real ival2                !" 固有値虚数部分(固有ベクトル用)

      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 ur(xgr,ygr)          !" 東西速度成分
      real ui(xgr,ygr)          !" 東西速度成分
      real vr(xgr,ygr)          !" 南北速度成分
      real vi(xgr,ygr)          !" 南北速度成分
      real pr(xgr,ygr)          !" ジオポテンシャル(phi)
      real ppr(xgr,ygr)         !" ジオポテンシャル(phi) ※一波長分データ作成用
      real uur(xgr,ygr)         !" 東西速度(u)           ※一波長分データ作成用
      real vvr(xgr,ygr)         !" 南北速度(v)           ※一波長分データ作成用
      real pim(xgr,ygr)         !" ジオポテンシャル(phi)
      real t                    !" 変動時間

      real k                    !" 波数(入力)
      real kgrid
      real a                    !" 地球の半径 

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

      character aaa
      character csgi*1
      character title1*10
      character title2*5
      character title3*3
      character title4*6
      character title5*5
      character title6*4
      character title7*6
      character title8*5
      character title9*6
      character title10*6
      character output*22       !" output ps file 名
      character output1*22      !" output ps file 名
      character output2*22      !" output ps file 名
      character output3*17      !" output ps file 名
      character output4*21      !" output ps file 名

      title2 = csgi(216)//'u/'//csgi(216)//'t'
      title3 = csgi(153)//'yv'
      title4 = '-'//csgi(216)//csgi(148)//'/'//csgi(216)//'x'
      title5 = csgi(216)//'v/'//csgi(216)//'t'
      title6 = '-'//csgi(153)//'yu'
      title7 = '-'//csgi(216)//csgi(148)//'/'//csgi(216)//'y'
      title8 =      csgi(216)//csgi(148)//'/'//csgi(216)//'t'
      title9 = '-'//csgi(216)//'u/'//csgi(216)//'x'
      title10= '-'//csgi(216)//'v/'//csgi(216)//'y'

*                   5   10   15   20   25   30   35   40   45   50
*               ----+----+----+----+----+----+----+----+----+----+'
      output1= 'n _k   _ig-w-east-term'
      output2= 'n _k   _ig-w-west-term'
      output3= 'n _k   _ro-w-term'
      output4= 'n-1k   _kelvin-w-term'

*-----------------------------------------------------------------------
*     1.0 描く波の選択と東西波数
*-----------------------------------------------------------------------

*---- !" 描く波の選択

      write(*,*) '書かせる波動の種類?'
      write(*,*) '(1) 東進慣性重力波(n＝0), (2) 東進慣性重力波(n≠0)'
      write(*,*) '(3) 混合ロスビー重力波(n＝0),(4) 西進慣性重力波(n≠0)'
      write(*,*) '(5) ケルビン波,           (6) ロスビー波'
      read(*,*) nwave

      write(*,*) 'モード番号 n=?'
      read(*,*) n

      write(*,*) '書かせる波数 k ?'
      read(*,*) k

      write(output1(5:7),'(f3.1)') k 
      write(output2(5:7),'(f3.1)') k 
      write(output3(5:7),'(f3.1)') k 
      write(output4(5:7),'(f3.1)') k 

*-----------------------------------------------------------------------
*     1.1 欠損値処理
*-----------------------------------------------------------------------

      call GLRGET( 'RMISS', RMISS )
      call GLLSET( 'LMISS', .TRUE. )

*-----------------------------------------------------------------------
*     1.2 座標データの作成
*-----------------------------------------------------------------------

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

*---- y 座標データ

      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 1210 ii=0,ygr

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

 1210 continue

      do 1220 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座標

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

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

 1220 continue


*-----------------------------------------------------------------------
*     3.2 モードの成分, それぞれの項のデータ作成(u,phi,v)
*-----------------------------------------------------------------------

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

      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
c      kgrid = k

*---- 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 座標 (ラジアン)
c         xxx(ix) = (beta*dg/gamma)*xx(ix) !" 無次元 x 座標 (経度)
         write(*,*) 'xxx(',ix,')=',xxx(ix)
 1230 continue

      XMIN  = xx(1)
      XMAX  = xx(xgr)

      write(*,*) 'XMIN=',XMIN
      write(*,*) 'XMAX=',XMAX

*---- !" 固有値理論解

      if (nwave.ne.5) then

      title1 = 'n= , k=   '

      write(title1(3:3) ,'(i1)'  ) n
      write(title1(8:10),'(f3.1)') k
      write(output1(2:2),'(i1)') n 
      write(output2(2:2),'(i1)') n 
      write(output3(2:2),'(i1)') n 

      if (nwave.eq.1) omega  = 0.5*( k + sqrt( k**2 + 4. ) )
      if (nwave.eq.1) output = output1 
      if (nwave.eq.2) omega  = sqrt(k**2 + 2.*n + 1.)
      if (nwave.eq.2) output = output1
      if (nwave.eq.3) omega  = 0.5*( k - sqrt( k**2 + 4. ) )
      if (nwave.eq.3) output = output2
      if (nwave.eq.4) omega = -sqrt(k**2 + 2.*n + 1.)
      if (nwave.eq.4) output = output2
      if (nwave.eq.6) omega = -k/(k**2 + 2.*n + 1.)
      if (nwave.eq.6) output = output3

      write(*,*) 'output=',output

*---- !" for inertio gravity & mixed rossby & rossby waves

      write(*,*) 
     & 'drowing for inertio gravity & mixed rossby & rossby waves'

      do 3120 iny=1,ygr    

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

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

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

 3120 continue         

      else 

*---- !" for kelvin

      if (nwave.eq.5) omega = real(k)
      if (nwave.eq.5) output = output4

      write(*,*) 'output=',output
      title1 = 'k=   '
      write(title1(3:5),'(f3.1)') k

      write(*,*) 
     & 'drowing for kelvin waves'

      do 3130 iny=0,ygr+1    

        vrvect(iny) = 0.
        vivect(iny) = 0.

        urvect(iny) = exp( -(1./2.)*y(iny)**2 )
        uivect(iny) = 0.

        prvect(iny) = exp( -(1./2.)*y(iny)**2 )
        pivect(iny) = 0.

 3130 continue         

      end if

        rval2 = omega
        ival2 = 0.

*---- !" モードの成分の計算(データ作成)

        write(*,*) 't=?'
        read(*,*) t
c        t = 0.

        do 3210 ix=1,xgr

           do 3220 iy = 1,ygr
              ur(ix,iy) = ( urvect(iy)*cos(k*xxx(ix)-rval2*t)
     &                    - uivect(iy)*sin(k*xxx(ix)-rval2*t) )

              ui(ix,iy) = ( urvect(iy)*sin(k*xxx(ix)-rval2*t)
     &                    + uivect(iy)*cos(k*xxx(ix)-rval2*t) )

              pr(ix,iy) = ( prvect(iy)*cos(k*xxx(ix)-rval2*t)
     &                    - pivect(iy)*sin(k*xxx(ix)-rval2*t) )

              pim(ix,iy)= ( prvect(iy)*sin(k*xxx(ix)-rval2*t)
     &                    + pivect(iy)*cos(k*xxx(ix)-rval2*t) )

              vr(ix,ivy)= ( vrvect(iy)*cos(k*xxx(ix)-rval2*t)
     &                    - vivect(iy)*sin(k*xxx(ix)-rval2*t) )

              vi(ix,ivy)= ( vrvect(iy)*sin(k*xxx(ix)-rval2*t)
     &                    + vivect(iy)*cos(k*xxx(ix)-rval2*t) )
 3220      continue

 3210   continue


*---- !" 二波長分だけのモードの成分の計算(データ作成) ※ 無次元波数で描く

        do 3270 iix = 1,xgr

           do 3250 iiy = 1,ygr 
              ppr(iix,iiy) = ( prvect(iiy)*cos(k*xxx(iix)-rval2*t)
     &                       - pivect(iiy)*sin(k*xxx(iix)-rval2*t) )

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

              vvr(iix,iiy)=  ( vrvect(iiy)*cos(k*xxx(iix)-rval2*t)
     &                       - vivect(iiy)*sin(k*xxx(iix)-rval2*t))
 3250      continue

 3270   continue

*---- ベクトルを描くためのデータ作成(y方向にまびく場合)

        do 3240 ixn=1,nvxg
           do 3240 iyn=1,nvyg
              uve(ixn,iyn)
     &             = uur(1+int((xgr-1)/(nvxg-1))*(ixn-1), 
     &               1+int((ygr-1)/(nvyg-1))*(iyn-1))
              vve(ixn,iyn)
     &             = vvr(1+int((xgr-1)/(nvxg-1))*(ixn-1), 
     &               1+int((ygr-1)/(nvyg-1))*(iyn-1))
 3240         continue

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

*---- 欠損値処理

      do inx=1,xgr          

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

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

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

      end do

      do inx=1,xgr
         do iny=2,ygr-1 

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

*---- ∂u/∂t 

      dudt(inx,iny) 
     &           = (
     &               ( 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) 
     &           = 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)
     &           = k*(
     &                 prvect(iny)*sin( k*xxx(inx) - rval2*t )
     &                +pivect(iny)*cos( k*xxx(inx) - rval2*t ) 
     &                )*exp( ival2*t )

c      write(*,*) 'dudt(',inx,iny,')=',dudt(inx,iny)
c      write(*,*) 'betayv(',inx,iny,')=',betayv(inx,iny)
c      write(*,*) 'dpdx(',inx,iny,')=',dpdx(inx,iny)
c      write(*,*) 'dudt - betayv + dpdx =',
c     &     dudt(inx,iny)-betayv(inx,iny)-dpdx(inx,iny)
c      write(*,*) ''

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

*---- ∂v/∂t 

      dvdt(inx,iny) 
     &           = (
     &               ( 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) 
     &           = -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)
     &           = -(
     &                ( ( 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 )

c      write(*,*) 'dvdt(',inx,iny,')=',dvdt(inx,iny)
c      write(*,*) 'betayu(',inx,iny,')=',betayu(inx,iny)
c      write(*,*) 'dpdy(',inx,iny,')=',dpdy(inx,iny)
c      write(*,*) 'dvdt - betayu + dpdy =',
c     &     dvdt(inx,iny)-betayu(inx,iny)-dpdy(inx,iny)
c      write(*,*) ''

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

*---- ∂Φ/∂t 

      dpdt(inx,iny) 
     &           = (
     &               ( 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) 
     &           = 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)
     &           = -( 
     &               ( ( 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 )

c      write(*,*) 'dpdt(',inx,iny,')=',dpdt(inx,iny)
c      write(*,*) 'ghdudx(',inx,iny,')=',ghdudx(inx,iny)
c      write(*,*) 'ghdvdy(',inx,iny,')=',ghdvdy(inx,iny)
c      write(*,*) 'dvdt - betayu + dpdy =',
c     &     dpdt(inx,iny)-ghdudx(inx,iny)-ghdvdy(inx,iny)
c      write(*,*) ''

      end do
      end do


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

      call swcset('FNAME',output) !" out put する ps ファイル名の指定

      call gropn(-iws)
*--------- マージン
      call slmgn(0.00, 0.00, 0.00, 0.00)

*--------- コーナマーカを消す
      call sglset('lcorner', .FALSE.)
*--------- フレーム分割
c      call sldiv ( 'S', 1, 4 )
c      call sglset( 'LFULL', .TRUE. )
c      call slrat( 1., 1.5 ) 
*--------- マージン
      call slmgn(0.00, 0.00, 0.00, 0.00)

*--------- 左 1フレーム目の設定 --------------------------------------

      call glrget('RUNDEF',RUNDEF)
      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.15, 0.40, 0.72, 0.92 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call uzlset('LABELXB', .FALSE.)
      call usxaxs('t')
      call usxaxs('b')
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELXB', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.8 )
      call uxsttl('t','U',0.0)
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
c      call uxsttl('b','LONGITUDE',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
      call uysttl('L','LAT.',0.0)
      call uzfact( 1/0.8 )
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if

*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
      call udrset('RSIZEL',0.010) 
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udcntr( uur, xgr, xgr, ygr )
c*--------- コンターインターバルの計算
c      rint = rudlev(1)
*--------- ベクトル図のマージン調整
c      call ugrset( 'RSIZET', 0.01 )
*--------- 単位ベクトルを描くルーチン
      call uglset( 'LUNIT' , .FALSE. )
c      call ugsut ( 'X', 'u' )
c      call ugsut ( 'Y', 'v' )
*--------- ベクトルのスケーリングの表示をしない
      call uglset('LMSG', .FALSE.)

*--------- ベクトル図を描くルーチン
c      call ugvect( uve, xgr, vve, xgr, nvxg, nvyg ) !" 間引く場合


*--------- 左 2フレーム目の設定

      call grfig
      call glrget('RUNDEF',RUNDEF)
c      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.15, 0.40, 0.51, 0.71 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call uzlset('LABELXB', .FALSE.)
      call usxaxs('t')
      call usxaxs('b')
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELXB', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.85 )
c      call uxsttl('t',title2,0.0)
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
c      call uxsttl('b','LONGITUDE',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
      call uysttl('L','LAT.',0.0)
      call uzfact( 1/0.85 )
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if
      
*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udrset('RSIZEL',0.01) 
*--------- コンターの描画
      call udcntr( dudt, xgr, xgr, ygr )
*--------- コンターインターバルの計算
      rint = rudlev(1)

*--------- 左 3フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
c      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.15, 0.40, 0.30, 0.50 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call uzlset('LABELXB', .FALSE.)
      call usxaxs('t')
      call usxaxs('b')
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELXB', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.85 )
c      call uxsttl('t',title3,0.0)
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
c      call uxsttl('b','LONGITUDE',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
      call uysttl('L','LAT.',0.0)
      call uzfact( 1/0.85 )
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if
      
*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udrset('RSIZEL',0.01) 
*--------- コンターラインの間隔指定(上図とコンターインターバルを揃える)
      call udgclb( betayv, xgr, xgr, ygr, rint ) 
*--------- コンターの描画
      call udcntr( betayv, xgr, xgr, ygr )

*--------- 左 4フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
c      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.15, 0.40, 0.09, 0.29 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

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

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.85 )
c      call uxsttl('t',title4,0.0)
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
      call uxsttl('b','LON.',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
      call uysttl('L','LAT.',0.0)
      call uzfact( 1/0.85 )
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if
      
*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udrset('RSIZEL',0.01) 
*--------- コンターラインの間隔指定(左上図とコンターインターバルを揃える)
      call udgclb( dpdx, xgr, xgr, ygr, rint ) 
*--------- コンターの描画
      call udcntr( dpdx, xgr, xgr, ygr )

*--------- udgcla で設定したコンターレベルを無効にする
      call udiclv

*--------- 真中 1 フレーム目の設定 --------------------------------------

      call grfig
      call glrget('RUNDEF',RUNDEF)
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.41, 0.66, 0.72, 0.92 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call uzlset('LABELXB', .FALSE.)
      call usxaxs('t')
      call usxaxs('b')
      call uzlset('LABELYL', .FALSE.)
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELXB', .TRUE.)
      call uzlset('LABELYL', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.8 )
      call uxsttl('t','V',0.0)
      call uxsttl('t',title1,0.0)
      call uzfact( 1/0.8 )
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
c      call uxsttl('b','LONGITUDE',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
c      call uysttl('L','LAT.',0.0)
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if

*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udrset('RSIZEL',0.01) 
*--------- コンターの描画
      call udcntr( vvr, xgr, xgr, ygr )
c*--------- コンターインターバルの計算
c      rint = rudlev(1)
*--------- ベクトル図のマージン調整
c      call ugrset( 'RSIZET', 0.01 )
*--------- 単位ベクトルを描くルーチン
      call uglset( 'LUNIT' , .FALSE. )
c      call ugsut ( 'X', 'u' )
c      call ugsut ( 'Y', 'v' )
*--------- ベクトルのスケーリングの表示をしない
      call uglset('LMSG', .FALSE.)

*--------- ベクトル図を描くルーチン
c      call ugvect( uve, xgr, vve, xgr, nvxg, nvyg ) !" 間引く場合

*--------- 真中 2フレーム目の設定

      call grfig
      call glrget('RUNDEF',RUNDEF)
c      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.41, 0.66, 0.51, 0.71 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call uzlset('LABELXB', .FALSE.)
      call usxaxs('t')
      call usxaxs('b')
      call uzlset('LABELYL', .FALSE.)
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELXB', .TRUE.)
      call uzlset('LABELYL', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.85 )
c      call uxsttl('t',title5,0.0)
      call uzfact( 1/0.85 )
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
c      call uxsttl('b','LONGITUDE',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
c      call uysttl('L','LAT.',0.0)
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if
      
*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udrset('RSIZEL',0.01) 
*--------- コンターの描画
      call udcntr( dvdt, xgr, xgr, ygr )
*--------- コンターインターバルの計算
      rint = rudlev(1)

*--------- 真中 3フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
c      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.41, 0.66, 0.30, 0.50 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call uzlset('LABELXB', .FALSE.)
      call usxaxs('t')
      call usxaxs('b')
      call uzlset('LABELYL', .FALSE.)
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELXB', .TRUE.)
      call uzlset('LABELYL', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.85 )
c      call uxsttl('t',title6,0.0)
      call uzfact( 1/0.85 )
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
c      call uxsttl('b','LONGITUDE',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
c      call uysttl('L','LAT.',0.0)
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if
      
*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udrset('RSIZEL',0.01) 
*--------- コンターラインの間隔指定(上図とコンターインターバルを揃える)
      call udgclb( betayu, xgr, xgr, ygr, rint ) 
*--------- コンターの描画
      call udcntr( betayu, xgr, xgr, ygr )

*--------- 真中 4フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
c      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.41, 0.66, 0.09, 0.29 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call usxaxs('t')
      call usxaxs('b')
      call uzlset('LABELYL', .FALSE.)
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELYL', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.85 )
c      call uxsttl('t',title7,0.0)
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
      call uxsttl('b','LON.',0.0)
      call uzfact( 1/0.85 )
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
c      call uysttl('L','LAT.',0.0)
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if
      
*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udrset('RSIZEL',0.01) 
*--------- コンターラインの間隔指定(上図とコンターインターバルを揃える)
      call udgclb( dpdy, xgr, xgr, ygr, rint ) 
*--------- コンターの描画
      call udcntr( dpdy, xgr, xgr, ygr )

*--------- udgcla で設定したコンターレベルを無効にする
      call udiclv

*--------- 右 1 フレーム目の設定 --------------------------------------

      call grfig
      call glrget('RUNDEF',RUNDEF)
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.67, 0.92, 0.72, 0.92 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call uzlset('LABELXB', .FALSE.)
      call usxaxs('t')
      call usxaxs('b')
      call uzlset('LABELYL', .FALSE.)
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELXB', .TRUE.)
      call uzlset('LABELYL', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.8 )
      call uxsttl('t','Phi',0.0)
      call uzfact( 1/0.8 )
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
c      call uxsttl('b','LONGITUDE',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
c      call uysttl('L','LAT.',0.0)
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if

*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
      call udrset('RSIZEL',0.01) 
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udcntr( ppr, xgr, xgr, ygr )
c*--------- コンターインターバルの計算
c      rint = rudlev(1)
*--------- ベクトル図のマージン調整
c      call ugrset( 'RSIZET', 0.01 )
*--------- 単位ベクトルを描くルーチン
      call uglset( 'LUNIT' , .FALSE. )
c      call ugsut ( 'X', 'u' )
c      call ugsut ( 'Y', 'v' )
*--------- ベクトルのスケーリングの表示をしない
      call uglset('LMSG', .FALSE.)

*--------- ベクトル図を描くルーチン
      call ugvect( uve, xgr, vve, xgr, nvxg, nvyg ) !" 間引く場合
c      call ugvect( uve, xgr, vve, xgr, xgr, ygr ) !" 間引かない場合

*--------- 右 2フレーム目の設定

      call grfig
      call glrget('RUNDEF',RUNDEF)
c      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.67, 0.92, 0.51, 0.71 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call uzlset('LABELXB', .FALSE.)
      call usxaxs('t')
      call usxaxs('b')
      call uzlset('LABELYL', .FALSE.)
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELXB', .TRUE.)
      call uzlset('LABELYL', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.85 )
c      call uxsttl('t',title8,0.0)
      call uzfact( 1/0.85 )
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
c      call uxsttl('b','LONGITUDE',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
c      call uysttl('L','LAT.',0.0)
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if
      
*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udrset('RSIZEL',0.01) 
*--------- コンターの描画
      call udcntr( dpdt, xgr, xgr, ygr )
*--------- コンターインターバルの計算
      rint = rudlev(1)

*--------- 右 3フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
c      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.67, 0.92, 0.30, 0.50 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call uzlset('LABELXB', .FALSE.)
      call usxaxs('t')
      call usxaxs('b')
      call uzlset('LABELYL', .FALSE.)
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELXB', .TRUE.)
      call uzlset('LABELYL', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.85 )
c      call uxsttl('t',title9,0.0)
      call uzfact( 1/0.85 )
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
c      call uxsttl('b','LONGITUDE',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
c      call uysttl('L','LAT.',0.0)
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if
      
*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
      call udrset('RSIZEL',0.01) 
*--------- コンターラインの間隔指定(左上図とコンターインターバルを揃える)
      call udgclb( ghdudx, xgr, xgr, ygr, rint ) 
*--------- コンターの描画
      call udcntr( ghdudx, xgr, xgr, ygr )

*--------- 右 4フレーム目の設定
      call grfig
      call glrget('RUNDEF',RUNDEF)
c      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      call grsvpt( 0.67, 0.92, 0.09, 0.29 )
      call usspnt(ygr,RUNDEF,y)
*--------- データの最大値・最小値を切りの良い数値に丸め作画範囲を決める
      call uspfit
*--------- 正規化変換を確定
      call grstrf

*--------- 座標軸(上下左右の枠)を描く
      call usxaxs('t')
      call usxaxs('b')
      call uzlset('LABELYL', .FALSE.)
      call usyaxs('l')
      call usyaxs('r')
      call uzlset('LABELYL', .TRUE.)

*--------- 添字表示の許可を与えるルーチン
      call sglset('lcntl', .TRUE.)
*--------- トップのタイトル指定
      call uzfact( 0.85 )
c      call uxsttl('t',title10,0.0)
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
      call uxsttl('b','LON.',0.0)
      call uzfact( 1/0.85 )
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
c      call uysttl('L','LAT.',0.0)
*--------- 等高線図のマージン調整
      call udrset( 'RSIZET', 0.01 )

*--------- 等高線図を描くルーチン

c      write(*,*) 'contour line をカラーで描くか?'
c      write(*,*) '(1) yes (2) no'
c      read(*,*) nnsel

      nnsel = 2

      if (nnsel.eq.1) then 
      call udiset('INDXMJ', 23) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 21) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      else
      call udiset('INDXMJ', 3) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      call udiset('INDXMN', 1) !" 主曲線 line index(1:黒,2:赤,3:緑,4:青,5:黄)
      end if
      
*--------- コンター間隔の表示をしない
      call udlset('LMSG', .FALSE.)
*--------- 格子点配列の格子点座標を最小値と最大値および格子点数で指定する
c      call udrset('RSIZEL',0.01) 
*--------- コンターラインの間隔指定(左上図とコンターインターバルを揃える)
      call udgclb( ghdvdy, xgr, xgr, ygr, rint ) 
*--------- コンターの描画
      call udcntr( ghdvdy, xgr, xgr, ygr )

*--------- udgcla で設定したコンターレベルを無効にする
      call udiclv
      call grcls

      stop
      end
