*表題 1D 非回転浅水重力波の水平構造の理論解を表示するプログラム * *履歴 2002/03/19 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 *---- !" 固有値理論解 c write(*,*) 'モード番号 n=?' c read(*,*) n title1 = 'k= ' c title1 = 'n= , k= ' c write(title1(3:3) ,'(i1)' ) n write(title1(3:5),'(f3.1)') k if (nwave.eq.1) omega = k if (nwave.eq.2) omega = -k do 3120 iny=1,ygr vrvect(iny) = 0. vivect(iny) = 0. urvect(iny) = 1. uivect(iny) = 0. prvect(iny) = 1. pivect(iny) = 0. 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