*表題  2D 非回転浅水重力波の水平構造の理論解を表示するプログラム
*
*履歴  2002/03/14 H.Taniguchi

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

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

      real XMIN, XMAX

      real gamma                !" 南北シアー (du/dy) (南北座標決定に使用)
      real beta                 !" df/dy              (南北座標決定に使用)
      real ys                   !" gamma/2*beta       (南北座標決定に使用)
      real yc                   !"                    (南北座標決定に使用)
      real dy                   !" 格子点間隔(m 単位)
      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 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                    !" 地球の半径 

      character aaa
      character csgi*1
      character title1*10

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

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

*---- y 座標データ

      gamma = 0.00001
      beta  = 2.29D-11

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

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

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

 1220 continue


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

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

      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

      kgrid = k

*---- x 座標データ

      dx = (dxlon*kgrid)/(xgr-1) !" kgrid 波長分描いた時の1グリッドあたり経度
      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

c      write(*,*) 'aaa = ?'
c      read(*,*) aaa

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

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

c      write(*,*) 'OK ?'
c      read(*,*) aaa

c      write(*,*) 'いまから k=',k,'モード成分を計算します. OK ?'
c      read(*,*) aaa

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

      write(*,*) '二次元非回転浅水重力波を描かせます.'
      write(*,*) '(1) ω > 0, (2) ω < 0'
      read(*,*) nwave

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

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

      title1 = 'n= , k=   '

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

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

      do 3120 iny=1,ygr    

        vrvect(iny) = sin(n*pi*y(iny)/l)
        vivect(iny) = 0.

        urvect(iny) = 0.
        uivect(iny) = (n*pi/(l*k*(1-(omega/k)**2)))*cos(n*pi*y(iny)/l)

        prvect(iny) = 0.
        pivect(iny) = (n*(omega/k)*pi/(l*k*(1-(omega/k)**2)))
     &       *cos(n*pi*y(iny)/l)

 3120 continue         

        rval2 = omega

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

        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方向にまびく場合)

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

           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

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

*---- ベクトルを描くためのデータ作成(まびかない場合)

c          do 3240 ixn=1,xgr
c              do 3240 iyn=1,ygr
c                 uve(ixn,iyn) = uur(ixn,iyn)
c                 vve(ixn,iyn) = vvr(ixn,iyn)
c 3240         continue

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

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

*--------- 1フレーム目の設定
      call glrget('RUNDEF',RUNDEF)
      call grfrm
*--------- ウインドウの設定
      call grswnd( XMIN, XMAX, RUNDEF, RUNDEF )
      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 )
      call uxsttl('t',title1,0.0)
      call uzfact( 1/0.85 )
*--------- 座標軸のタイトル・単位指定 (X-title,X単位) 
      call uxsttl('b','X',0.0)
c      call uxsttl('b','LONGITUDE',0.0)
*--------- 座標軸のタイトル・単位指定 (Y-title,Y単位) 
      call uysttl('L','Y',0.0)
c      call uysttl('L','LATITUDE',0.0)
*--------- 等高線図のマージン調整
      call udrset( 'RSIZEL', 0.015 )
      call udrset( 'RSIZET', 0.015 )

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

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

      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 udcntr( ppr, xgr, xgr, ygr )
*--------- ベクトル図のマージン調整
      call ugrset( 'RSIZET', 0.015 )
*--------- 単位ベクトルを描くルーチン
      call uglset( 'LUNIT' , .TRUE. )
      call ugsut ( 'X', 'u' )
      call ugsut ( 'Y', 'v' )

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

      call grcls
      stop
      end
