!---------------------------------------------------------------------- ! Copyright (c) 2005 Masatsugu Odaka. All rights reserved. !---------------------------------------------------------------------- ! Sample program for gt4f90io and ISPACK 2005/05/16 M. Odaka ! ! Solving a non-linear 2-D shallow water system on a sphere. ! ! du/dt - 2\Omega\sin\phi v = -g/a\cos\phi dh/d\lambda + D(u), ! dv/dt + 2\Omega\sin\phi u = -g/a dh/d\phi + D(v) ! dH/dt + ( 1/\cos\phi d([H-hs] u)/d\lambda ! + 1/\cos\phi d([H-hs] v\cos\phi)/d\phi ) = D(h). ! ! where hs is hight of topography, D(a) is hyper viscosity. ! ! The time integration is performed by using leap-frog & trapezoidal implicit ! scheme. ! program shallow_zd_bench06 !== モジュール引用宣言 ================================================ use w_module use gt4_history, only : HistoryCreate, HistoryPut, HistoryGet, & & HistoryAddVariable, HistoryClose, GT_HISTORY use dc_trace, only : SetDebug, BeginSub, EndSub, DbgMessage use dc_message, only : MessageNotify implicit none !== 宣言部 ============================================================ !---------------------------------------------------------------------- ! NAMELIST ファイル名 !---------------------------------------------------------------------- character(len=100) :: NameListFile='shallow_zd_bench06.nml' !---------------------------------------------------------------------- ! 変数(格子点データ) !---------------------------------------------------------------------- real(8), allocatable :: xy_VelLon_n(:,:) ! 速度経度成分 (t) real(8), allocatable :: xy_VelLat_n(:,:) ! 速度緯度成分 (t) real(8), allocatable :: xy_Hsfc_n(:,:) ! 表面変位 (t) real(8), allocatable :: xy_Vor_n(:,:) ! 渦度 (t) real(8), allocatable :: xy_HDiv_n(:,:) ! 発散 (t) real(8), allocatable :: xy_VelLon_b(:,:) ! 速度経度成分 (t-Δt) real(8), allocatable :: xy_VelLat_b(:,:) ! 速度緯度成分 (t-Δt) real(8), allocatable :: xy_Hsfc_b(:,:) ! 表面変位 (t-Δt) real(8), allocatable :: xy_Vor_b(:,:) ! 渦度 (t-Δt) real(8), allocatable :: xy_HDiv_b(:,:) ! 発散 (t-Δt) real(8), allocatable :: xy_Coli(:,:) ! 惑星渦度 real(8), allocatable :: xy_Topo(:,:) ! 表面地形 !---------------------------------------------------------------------- ! 変数(スペクトルデータ) !---------------------------------------------------------------------- real(8), allocatable :: w_Hsfc_a(:) ! 変位 (t+Δt) real(8), allocatable :: w_Vor_a(:) ! 渦度 (t+Δt) real(8), allocatable :: w_HDiv_a(:) ! 発散 (t+Δt) real(8), allocatable :: w_Hsfc_n(:) ! 変位 (t) real(8), allocatable :: w_Vor_n(:) ! 渦度 (t) real(8), allocatable :: w_HDiv_n(:) ! 発散 (t) real(8), allocatable :: w_Hsfc_b(:) ! 変位 (t-Δt) real(8), allocatable :: w_Vor_b(:) ! 渦度 (t-Δt) real(8), allocatable :: w_HDiv_b(:) ! 発散 (t-Δt) real(8), allocatable :: w_dHsfc(:) ! 変位増分 real(8), allocatable :: w_dVor(:) ! 渦度増分 real(8), allocatable :: w_dHDiv(:) ! 発散増分 !---------------------------------------------------------------------- ! 固定パラメタ !---------------------------------------------------------------------- real(8), parameter :: pi = 3.141592653589793D0 ! 円周率 !---------------------------------------------------------------------- ! gtool4 用変数(複数ファイルに分割出力するため) !---------------------------------------------------------------------- type(GT_HISTORY) :: hst_2d, hst_1d, hst_rst !---------------------------------------------------------------------- ! NAMELIST 変数 (実験設定等) !---------------------------------------------------------------------- character(len=100) :: ExpTitle='' ! 実験名 character(len=100) :: ExpInst ='' ! 実行者所属 character(len=100) :: ExpSrc ='' ! 実行者名 character(len=100) :: ExpModel='' ! モデル設定名 namelist /expset/ ExpTitle, ExpInst, ExpSrc, ExpModel !---------------------------------------------------------------------- ! NAMELIST 変数 (時間単位) !---------------------------------------------------------------------- integer :: UnitYear = 360 ! 1 年の長さ(日) integer :: UnitDay = 24 ! 1 日の長さ(時間) integer :: UnitHour = 60 ! 1 時間の長さ(分) integer :: UnitMinit = 60 ! 1 分の長さ(秒) namelist /dateset/ UnitYear, UnitDay, UnitHour, UnitMinit !---------------------------------------------------------------------- ! NAMELIST 変数 (時間積分パラメタ) !---------------------------------------------------------------------- real(8) :: TimeIntYear ! 積分時間(年) real(8) :: TimeIntDay ! 積分時間(日) real(8) :: TimeIntHour ! 積分時間(時) real(8) :: TimeIntMinit ! 積分時間(分) real(8) :: TimeIntSec ! 積分時間(秒) real(8) :: OutputYear ! 出力間隔(年) real(8) :: OutputDay ! 出力間隔(日) real(8) :: OutputHour ! 出力間隔(時) real(8) :: OutputMinit ! 出力間隔(分) real(8) :: OutputSec ! 出力間隔(秒) real(8) :: DelTime ! 時間刻み(秒) character(len=100) :: IntScheme='' ! 時間積分法 namelist /timeset/ & & TimeIntYear, TimeIntDay, TimeIntHour, TimeIntMinit, TimeIntSec, & & OutputYear, OutputDay, OutputHour, OutputMinit, OutputSec, & & DelTime, IntScheme !---------------------------------------------------------------------- ! NAMELIST 変数 (入出力ファイル名) !---------------------------------------------------------------------- character(len=100) :: InputRstFile='' ! 入力リスタートファイル名 character(len=100) :: OutputRstfile='' ! 出力リスタートファイル名 character(len=100) :: OutputFile1D='' ! 出力ファイル名 character(len=100) :: OutputFile2D='' ! 出力ファイル名 namelist /fileset/ InputRstFile, OutputRstfile, OutputFile1D, OutputFile2D !---------------------------------------------------------------------- ! NAMELIST 変数(空間解像度設定) !---------------------------------------------------------------------- integer :: nm=21 ! 切断全波数 integer :: im=64 ! 経度方向格子点数 ( > 3*nm + 1) integer :: jm=32 ! 緯度方向格子点数 ( > 3*nm/2 ) namelist /gridset/ nm, im, jm !---------------------------------------------------------------------- ! NAMELIST 変数(物理パラメタ) !---------------------------------------------------------------------- real(8) :: Rplanet = 6.37122D6 ! 惑星半径 real(8) :: Omega = 7.292D-5 ! 回転角速度 real(8) :: Grav = 9.80616D0 ! 重力加速度 real(8) :: TfilCoef = 0.05 ! 時間フィルター係数 integer :: VisOrder = 2 ! 超粘性の次数 real(8) :: VisCoef = 2.338D16 ! 超粘性係数 namelist /paramset/ Rplanet, Omega, Grav, & & TfilCoef, VisOrder, VisCoef !---------------------------------------------------------------------- ! NAMELIST 変数 (デバッグ出力制御) !---------------------------------------------------------------------- logical :: DebugOn = .false. ! デバッグ出力制御 namelist /debugset/ DebugOn !---------------------------------------------------------------------- ! 内部変数 !---------------------------------------------------------------------- real(8) :: delta ! 浅水/順圧切替え制御 real(8) :: alpha ! 時間差分切替え制御 real(8) :: beta ! 陰解法用パラメタ real(8) :: GeoPotAvr ! 平均ポテンシャル real(8) :: HsfcAvr ! 平均ポテンシャル real(8) :: Time ! 時間(秒) integer :: NStep ! 時間ステップ数 integer :: NStepAll ! 全時間ステップ数 integer :: NStepOut ! 出力ステップ間隔 real(8) :: InitTime ! 積分開始時刻 real(8) :: TimeIntAll ! 全積分時間(秒) real(8) :: OutputAll ! 全出力間隔(秒) real(8) :: MassAll ! 全質量 real(8) :: EnergyAll ! 全エネルギー real(8) :: AngMomentAll ! 全角運動量 real(8) :: PotEnstAll !全ポテンシャルエンストロフィー real(8) :: VorAll ! 全渦度 real(8) :: HDivAll ! 全発散 real(8) :: CpuTime ! 計算 CPU 時間(frt 用) character(len=20) :: DbgMessageFmt='*** DbgMESSAGE *** ' !---------------------------------------------------------------------- ! 内部変数 (Case5 実験用) !---------------------------------------------------------------------- real(8), allocatable :: xy_VelLonRef(:,:) ! 速度経度成分(参照解) real(8), allocatable :: xy_VelLatRef(:,:) ! 速度緯度成分(参照解) real(8), allocatable :: xy_HsfcRef(:,:) ! 表面変位(参照解) real(8) :: MassAllRef ! 全質量(参照解) real(8) :: EnergyAllRef ! 全エネルギー(参照解) real(8) :: PotEnstAllRef !全ポテンシャルエンストロフィー ! (参照解) real(8) :: MassAllError ! 全質量の誤差 real(8) :: EnergyAllError ! 全エネルギーの誤差 real(8) :: PotEnstAllError !全ポテンシャルエンストロフィー ! の誤差 !== 初期化 ============================================================ !---------------------------------------------------------------------- ! NAMELIST 読み込み !---------------------------------------------------------------------- open(10,file=trim(NameListFile)) read(10,nml=expset) read(10,nml=dateset) read(10,nml=timeset) read(10,nml=fileset) read(10,nml=gridset) read(10,nml=paramset) read(10,nml=debugset) close(10) !---------------------------------------------------------------------- ! デバッグ出力制御設定 ! DebugOn が .true. なら SetDebug を呼ぶ !---------------------------------------------------------------------- if (DebugOn) then call SetDebug end if !---------------------------------------------------------------------- ! 動的割り付け配列変数の初期化 !---------------------------------------------------------------------- allocate(xy_VelLon_n(im,jm), xy_VelLat_n(im,jm), xy_Hsfc_n(im,jm), & & xy_Vor_n(im,jm) , xy_HDiv_n(im,jm) , & & xy_VelLon_b(im,jm), xy_VelLat_b(im,jm), xy_Hsfc_b(im,jm), & & xy_Vor_b(im,jm) , xy_HDiv_b(im,jm)) allocate(xy_Coli(im,jm), xy_Topo(im,jm)) allocate(w_Hsfc_a((nm+1)*(nm+1)), & & w_Vor_a((nm+1)*(nm+1)) , & & w_HDiv_a((nm+1)*(nm+1)), & & w_Hsfc_n((nm+1)*(nm+1)), & & w_Vor_n((nm+1)*(nm+1)) , & & w_HDiv_n((nm+1)*(nm+1)), & & w_Hsfc_b((nm+1)*(nm+1)), & & w_Vor_b((nm+1)*(nm+1)) , & & w_HDiv_b((nm+1)*(nm+1))) allocate(w_dHsfc((nm+1)*(nm+1)), & & w_dVor((nm+1)*(nm+1)) , & & w_dHDiv((nm+1)*(nm+1))) allocate(xy_VelLonRef(im,jm), & & xy_VelLatRef(im,jm), & & xy_HsfcRef(im,jm) ) !---------------------------------------------------------------------- ! spml ライブラリの初期化 ! 座標変数 xy_Lon, xy_Lat x_Lon, y_Lat はここで値が確定する. !---------------------------------------------------------------------- call DbgMessage(fmt='call %c', c1='w_initial') call w_initial(nm,im,jm) ! ISPACK初期化 !---------------------------------------------------------------------- ! リスタートファイル読み込み ! Hsfc = Hsfc - HsfcAvr とする. !---------------------------------------------------------------------- xy_Coli = 2 * Omega * sin(xy_Lat) call InputRestart_gtool4(InputRstFile, InitTime, xy_Topo, & & xy_VelLon_n, xy_VelLat_n, xy_Hsfc_n, xy_Vor_n, xy_HDiv_n, & & xy_VelLon_b, xy_VelLat_b, xy_Hsfc_b, xy_Vor_b, xy_HDiv_b) HsfcAvr = AvrLonLat_xy(xy_Hsfc_n) GeoPotAvr = Grav*HsfcAvr xy_Hsfc_n = xy_Hsfc_n - HsfcAvr xy_Hsfc_b = xy_Hsfc_b - HsfcAvr w_Hsfc_n = w_xy(xy_Hsfc_n) ; w_Hsfc_b = w_xy(xy_Hsfc_b) ; w_Hsfc_a = w_Hsfc_n w_Vor_n = w_xy(xy_Vor_n) ; w_Vor_b = w_xy(xy_Vor_b) ; w_Vor_a = w_Vor_n w_HDiv_n = w_xy(xy_HDiv_n) ; w_HDiv_b = w_xy(xy_HDiv_b) ; w_HDiv_a = w_HDiv_n call Output_gtool4_init ! ヒストリー初期化 !---------------------------------------------------------------------- ! 参照解の保存: 初期値を保存しておく !---------------------------------------------------------------------- xy_HsfcRef = xy_Hsfc_n xy_VelLonRef = xy_VelLon_n xy_VelLatRef = xy_VelLat_n MassAllRef = IntLonLat_xy(xy_HsfcRef + HsfcAvr - xy_Topo) EnergyAllRef = IntLonLat_xy((xy_HsfcRef + HsfcAvr - xy_Topo)*( & & xy_VelLonRef**2 + xy_VelLatRef**2 + & & Grav*(xy_HsfcRef + HsfcAvr - xy_Topo))*0.5D0 ) PotEnstAllRef = IntLonLat_xy((xy_Vor_n + xy_Coli)**2/ & & (xy_HsfcRef + HsfcAvr - xy_Topo)*0.5D0) !---------------------------------------------------------------------- ! 浅水/順圧切替え ! 変数 ExpModel で切替え. デフォルトは浅水モデル. !---------------------------------------------------------------------- select case (ExpModel) case ('baro') delta = 0.0D0 call MessageNotify('Message', & & 'shallow_zd', & & 'delta=%f; %c', & & d=(/delta/), & & c1='barotoropic vorticity equation') case default delta = 1.0D0 call MessageNotify('Message', & & 'shallow_zd', & & 'delta=%f; %c', & & d=(/delta/), & & c1='shallow water equation') end select !---------------------------------------------------------------------- ! 時間積分スキームの切替えパラメータの設定 ! 変数 IntScheme で切替え. デフォルトは semi-implicit !---------------------------------------------------------------------- select case (IntScheme) case ('explicit') alpha = 0.0D0 beta = 0.0D0 call MessageNotify('Message', & & 'shallow_zd', & & 'alpha=%f; %c', & & d=(/alpha/), & & c1='explicit time integration') case default alpha = 1.0D0 beta = GeoPotAvr*DelTime**2/Rplanet**2 call MessageNotify('Message', & & 'shallow_zd', & & 'alpha=%f; %c', & & d=(/alpha/), & & c1='semi-implicit time integration') end select !---------------------------------------------------------------------- ! 時間積分設定 ! NAMELIST で入力した全積分時間と出力時間間隔を年月日分秒から秒へ変換し, ! DelTime を用いて全計算ステップ数 NstepAll と出力ステップ数 NstepOut ! を計算する. !---------------------------------------------------------------------- Time = InitTime if (DelTime == 0.0) then call MessageNotify('Error', & & 'shallow_zd', & & 'DelTime is 0.') end if TimeIntAll = (TimeIntYear * UnitYear * UnitDay * UnitHour * UnitMinit & & + TimeIntDay * UnitDay * UnitHour * UnitMinit & & + TimeIntHour * UnitHour * UnitMinit & & + TimeIntMinit * UnitMinit & & + TimeIntSec) if (mod(TimeIntAll, DelTime) /= 0) then call MessageNotify('Error', & & 'shallow_zd', & & 'mod(TimeIntAll, DelTime) is not 0.') else NstepAll = INT(TimeIntAll / DelTime) end if OutputAll = (OutputYear * UnitYear * UnitDay * UnitHour * UnitMinit & & + OutputDay * UnitDay * UnitHour * UnitMinit & & + OutputHour * UnitHour * UnitMinit & & + OutputMinit * UnitMinit & & + OutputSec) if (mod(OutputAll, DelTime) /= 0) then call MessageNotify('Error', & & 'shallow_zd', & & 'mod(OutputAll, DelTime) is not 0.') else NstepOut = INT(OutputAll / DelTime) end if call DbgMessage(fmt='%c TimeIntAll=%f', d=(/TimeIntALL/), c1=DbgMessageFmt) call DbgMessage(fmt='%c OutPutAll=%f', d=(/OutPutALL/), c1=DbgMessageFmt) call DbgMessage(fmt='%c NstepAll=%d', i=(/NstepALL/), c1=DbgMessageFmt) call DbgMessage(fmt='%c NstepOut=%d', i=(/NstepOut/), c1=DbgMessageFmt) !== 時間積分 ========================================================== call DbgMessage(fmt='%c %c', & & c1=DbgMessageFmt, & & c2='Time integration starts.') !---------------------------------------------------------------------- ! 最初の 1 ステップ分の時間積分開始 ! InitTime = 0.0 の場合オイラースキームで計算 ! その他の場合はリープフロッグスキームで計算 !---------------------------------------------------------------------- Time = Time + DelTime call DbgMessage(fmt='%c Time=%f', d=(/Time/), c1=DbgMessageFmt) if (InitTime == 0.0) then call DbgMessage(fmt='%c InitTime=%f : %c', & & d=(/InitTime/), & & c1=DbgMessageFmt, & & c2='Euler scheme is adapted.') !---------------------------------------------------------------------- ! 増分 w_dHDiv, w_dVor, w_dHsfc を計算 !---------------------------------------------------------------------- call get_dVor(DelTime, xy_VelLon_n, xy_VelLat_n, xy_Vor_n, w_dVor) call get_dHDiv(DelTime, xy_VelLon_n, xy_VelLat_n, & & xy_Hsfc_n, xy_Vor_n, w_Hsfc_n, w_dHDiv) call get_dHsfc(DelTime, xy_VelLon_n, xy_VelLat_n, & & xy_Hsfc_n, xy_Vor_n, w_HDiv_n, w_dHsfc) !---------------------------------------------------------------------- ! 1 ステップ計算 !---------------------------------------------------------------------- w_Vor_a = w_Vor_n + w_dVor + 2.0D0*DelTime*w_NumVis_w(w_Vor_n) w_HDiv_a = ((1.0D0 + beta*rn(:,1))*w_HDiv_n + & & w_dHDiv + DelTime*w_NumVis_w(w_HDiv_n)) & & /(1.0D0 - beta*rn(:,1)) w_Hsfc_a = ((1.0D0 + beta*rn(:,1))*w_Hsfc_n + & & w_dHsfc + DelTime*w_NumVisScaler_w(w_Hsfc_n))& & /(1.0D0 - beta*rn(:,1)) w_Vor_b = w_Vor_n ; w_Vor_n = w_Vor_a w_HDiv_b = w_HDiv_n ; w_HDiv_n = w_hDiv_a w_HSfc_b = w_HSfc_n ; w_Hsfc_n = w_Hsfc_a else !---------------------------------------------------------------------- ! 増分 w_dHDiv, w_dVor, w_dHsfc を計算 ! 変数 IntScheme で切替え. デフォルトは semi-implicit !---------------------------------------------------------------------- call get_dVor(2.0*DelTime, xy_VelLon_n, xy_VelLat_n, xy_Vor_n, w_dVor) select case (IntScheme) case ('explicit') call get_dHDiv(2.0*DelTime, xy_VelLon_n, xy_VelLat_n, & & xy_Hsfc_n, xy_Vor_n, w_Hsfc_n, w_dHDiv) call get_dHsfc(2.0*DelTime, xy_VelLon_n, xy_VelLat_n, & & xy_Hsfc_n, xy_Vor_n, w_HDiv_n, w_dHsfc) case default call get_dHDiv(2.0*DelTime, xy_VelLon_n, xy_VelLat_n, & & xy_Hsfc_n, xy_Vor_n, w_Hsfc_b, w_dHDiv) call get_dHsfc(2.0D0*DelTime, xy_VelLon_n, xy_VelLat_n, & & xy_Hsfc_n, xy_Vor_n, w_HDiv_b, w_dHsfc) end select !---------------------------------------------------------------------- ! 1 ステップ計算 !---------------------------------------------------------------------- w_Vor_a = w_Vor_b + w_dVor + 2.0D0*DelTime*w_NumVis_w(w_Vor_b) w_HDiv_a = ((1.0D0 + beta*rn(:,1))*w_HDiv_b + & & w_dHDiv + 2.0D0*DelTime*w_NumVis_w(w_HDiv_b)) & & /(1.0D0 - beta*rn(:,1)) w_Hsfc_a = ((1.0D0 + beta*rn(:,1))*w_Hsfc_b + & & w_dHsfc + 2.0D0*DelTime*w_NumVisScaler_w(w_Hsfc_b))& & /(1.0D0 - beta*rn(:,1)) !---------------------------------------------------------------------- ! 時間フィルターをかけて, 変数を 1 ステップ更新する. !---------------------------------------------------------------------- call timefilter(w_Vor_a , w_Vor_n , w_Vor_b ) call timefilter(w_HDiv_a, w_HDiv_n, w_HDiv_b) call timefilter(w_Hsfc_a, w_Hsfc_n, w_Hsfc_b) end if !---------------------------------------------------------------------- ! 格子点データの計算 !---------------------------------------------------------------------- call uv(w_Vor_n, w_HDiv_n, xy_VelLon_n, xy_VelLat_n) xy_Vor_n = xy_w(w_Vor_n) xy_HDiv_n = xy_w(w_HDiv_n) xy_Hsfc_n = xy_w(w_Hsfc_n) call uv(w_Vor_b, w_HDiv_b, xy_VelLon_b, xy_VelLat_b) xy_Vor_b = xy_w(w_Vor_b) xy_HDiv_b = xy_w(w_HDiv_b) xy_Hsfc_b = xy_w(w_Hsfc_b) !---------------------------------------------------------------------- ! 2 ステップ以降の時間積分開始 !---------------------------------------------------------------------- do Nstep = 2, NstepALL Time = Time + DelTime call DbgMessage(fmt='%c Time=%f', d=(/Time/), c1=DbgMessageFmt) !---------------------------------------------------------------------- ! 増分 w_dHDiv, w_dVor, w_dHsfc を計算 ! 変数 IntScheme で切替え. デフォルトは semi-implicit !---------------------------------------------------------------------- call get_dVor(2.0*DelTime, xy_VelLon_n, xy_VelLat_n, xy_Vor_n, w_dVor) select case (IntScheme) case ('explicit') call get_dHDiv(2.0*DelTime, xy_VelLon_n, xy_VelLat_n, & & xy_Hsfc_n, xy_Vor_n, w_Hsfc_n, w_dHDiv) call get_dHsfc(2.0*DelTime, xy_VelLon_n, xy_VelLat_n, & & xy_Hsfc_n, xy_Vor_n, w_HDiv_n, w_dHsfc) case default call get_dHDiv(2.0*DelTime, xy_VelLon_n, xy_VelLat_n, & & xy_Hsfc_n, xy_Vor_n, w_Hsfc_b, w_dHDiv) call get_dHsfc(2.0D0*DelTime, xy_VelLon_n, xy_VelLat_n, & & xy_Hsfc_n, xy_Vor_n, w_HDiv_b, w_dHsfc) end select !---------------------------------------------------------------------- ! 1 ステップ計算 !---------------------------------------------------------------------- w_Vor_a = w_Vor_b + w_dVor + 2.0D0*DelTime*w_NumVis_w(w_Vor_b) w_HDiv_a = ((1.0D0 + beta*rn(:,1))*w_HDiv_b + & & w_dHDiv + 2.0D0*DelTime*w_NumVis_w(w_HDiv_b)) & & /(1.0D0 - beta*rn(:,1)) w_Hsfc_a = ((1.0D0 + beta*rn(:,1))*w_Hsfc_b + & & w_dHsfc + 2.0D0*DelTime*w_NumVisScaler_w(w_Hsfc_b))& & /(1.0D0 - beta*rn(:,1)) !---------------------------------------------------------------------- ! 時間フィルターをかけて, 変数を 1 ステップ更新する. ! * Case 1 では w_Vor と w_Hdiv は定常 !---------------------------------------------------------------------- call timefilter(w_Vor_a , w_Vor_n , w_Vor_b ) call timefilter(w_HDiv_a, w_HDiv_n, w_HDiv_b) call timefilter(w_Hsfc_a, w_Hsfc_n, w_Hsfc_b) !---------------------------------------------------------------------- ! 格子点データの計算 !---------------------------------------------------------------------- call uv(w_Vor_n, w_HDiv_n, xy_VelLon_n, xy_VelLat_n) xy_Vor_n = xy_w(w_Vor_n) xy_HDiv_n = xy_w(w_HDiv_n) xy_Hsfc_n = xy_w(w_Hsfc_n) call uv(w_Vor_b, w_HDiv_b, xy_VelLon_b, xy_VelLat_b) xy_Vor_b = xy_w(w_Vor_b) xy_HDiv_b = xy_w(w_HDiv_b) xy_Hsfc_b = xy_w(w_Hsfc_b) !---------------------------------------------------------------------- ! 出力 ! Nstep が NstepOut の整数倍になったら出力を行う. !---------------------------------------------------------------------- if(mod(Nstep, NStepOut) == 0) then !-- 保存量の計算 MassAll = IntLonLat_xy(xy_Hsfc_n + HsfcAvr - xy_Topo) EnergyAll = IntLonLat_xy((xy_Hsfc_n + HsfcAvr - xy_Topo)*( & & xy_VelLon_n**2 + xy_VelLat_n**2 + & & Grav*(xy_Hsfc_n + HsfcAvr - xy_Topo))*0.5D0 ) AngMomentAll = IntLonLat_xy((xy_Hsfc_n+HsfcAvr-xy_Topo)* & & cos(xy_Lat)*Rplanet & & *(xy_VelLon_n + Omega*cos(xy_Lat)*Rplanet)) PotEnstAll = IntLonLat_xy((xy_Vor_n + xy_Coli)**2/ & & (xy_Hsfc_n + HsfcAvr - xy_Topo)*0.5D0) VorAll = IntLonLat_xy(xy_Vor_n) HdivAll = IntLonLat_xy(xy_Hdiv_n) !--誤差の計算 MassAllError = (MassAll - MassAllRef)/MassAll EnergyAllError = (EnergyAll - EnergyAllRef)/EnergyAll PotEnstAllError= (PotEnstAll - PotEnstAllRef)/PotEnstAll ! call MessageNotify('Message', & !& 'shallow_zd', & !& 'Time=%f', d=(/Time/)) !-- 出力 call Output_gtool4 endif enddo call DbgMessage(fmt='%c %c', & & c1=DbgMessageFmt, & & c2='Time integration end.') !== 後処理 ========================================================== !---------------------------------------------------------------------- ! ヒストリー出力後処理 !---------------------------------------------------------------------- call Output_gtool4_close !---------------------------------------------------------------------- ! リスタートファイル出力 !---------------------------------------------------------------------- call OutputRestart_gtool4_init ! リスタートファイル初期化 call OutputRestart_gtool4 ! リスタートファイル出力 call OutputRestart_gtool4_close ! リスタートファイル出力終了 !---------------------------------------------------------------------- ! 動的割り付け配列変数の解放 !---------------------------------------------------------------------- deallocate(xy_VelLon_n, xy_VelLat_n, xy_Hsfc_n, & & xy_Vor_n , xy_HDiv_n , & & xy_VelLon_b, xy_VelLat_b, xy_Hsfc_b, & & xy_Vor_b , xy_HDiv_b , & & xy_Coli , xy_Topo ) deallocate(w_Hsfc_a, w_Vor_a, w_HDiv_a, & & w_Hsfc_n, w_Vor_n, w_HDiv_n, & & w_Hsfc_b, w_Vor_b, w_HDiv_b ) deallocate(w_dHsfc, & & w_dVor , & & w_dHDiv ) deallocate(xy_VelLonRef, & & xy_VelLatRef, & & xy_HsfcRef ) !---------------------------------------------------------------------- ! 動的割り付け配列変数の解放 !---------------------------------------------------------------------- !-- CPU 時間の出力 (Fujitsu frt 用) ! call clock(CpuTime, 0, 2) ! ! call MessageNotify('Message', & !& 'shallow_zd_benc01', & !& 'CpuTime=%f', d=(/CpuTime/)) stop contains !== 内部副プログラム ================================================= !---------------------------------------------------------------------- ! get_dVor ! 渦度方程式の増分 w_dVor を計算する !---------------------------------------------------------------------- subroutine get_dVor(dt, xy_VelLon, xy_VelLat, xy_Vor, w_dVor) real(8),intent(in) :: dt real(8),intent(in) :: xy_VelLon(im,jm) real(8),intent(in) :: xy_VelLat(im,jm) real(8),intent(in) :: xy_Vor(im,jm) real(8),intent(out):: w_dVor((nm+1)*(nm+1)) call BeginSub('get_dVor') w_dVor = dt*( & & - w_Div_xy_xy((xy_Coli + xy_Vor)*xy_VelLon, & & (xy_Coli + xy_Vor)*xy_VelLat) / Rplanet) call EndSub('get_dVor') end subroutine get_dVor !---------------------------------------------------------------------- ! get_dHDiv ! 発散の式の増分 w_dHDiv を計算する !---------------------------------------------------------------------- subroutine get_dHDiv(dt, xy_VelLon, xy_VelLat, xy_Hsfc, & & xy_Vor, w_Hsfc, w_dHDiv) real(8),intent(in) :: dt real(8),intent(in) :: xy_VelLon(im,jm) real(8),intent(in) :: xy_VelLat(im,jm) real(8),intent(in) :: xy_Hsfc(im,jm) real(8),intent(in) :: xy_Vor(im,jm) real(8),intent(in) :: w_Hsfc((nm+1)*(nm+1)) real(8),intent(out):: w_dHDiv((nm+1)*(nm+1)) call BeginSub('get_dHDiv') w_dHDiv = dt*( & & w_Div_xy_xy( (xy_Coli + xy_Vor)*xy_VelLat, & & - (xy_Coli + xy_Vor)*xy_VelLon)/Rplanet & & - w_Lapla_w(Grav*w_Hsfc & & + 0.5D0*w_xy(xy_VelLon**2 + xy_VelLat**2))/Rplanet**2 & & + alpha*DelTime*Grav*w_Lapla_w( & & w_Div_xy_xy((xy_Hsfc - xy_Topo)*xy_VelLon, & & (xy_Hsfc - xy_Topo)*xy_VelLat)/Rplanet)/Rplanet**2) call EndSub('get_dHDiv') end subroutine get_dHDiv !---------------------------------------------------------------------- ! get_dHsfc ! 質量保存の式の増分 w_dHsfc を計算する !---------------------------------------------------------------------- subroutine get_dHsfc(dt, xy_VelLon, xy_VelLat, xy_Hsfc, & & xy_Vor, w_HDiv, w_dHsfc) real(8),intent(in) :: dt real(8),intent(in) :: xy_VelLon(im,jm) real(8),intent(in) :: xy_VelLat(im,jm) real(8),intent(in) :: xy_Hsfc(im,jm) real(8),intent(in) :: xy_Vor(im,jm) real(8),intent(in) :: w_HDiv((nm+1)*(nm+1)) real(8),intent(out):: w_dHsfc((nm+1)*(nm+1)) call BeginSub('get_dHsfc') w_dHsfc = dt*( & & - w_Div_xy_xy((xy_Hsfc - xy_Topo)*xy_VelLon, & & (xy_Hsfc - xy_Topo)*xy_VelLat)/Rplanet & & - HsfcAvr*w_HDiv & & - alpha*DelTime*HsfcAvr*( & & w_Div_xy_xy( (xy_Coli + xy_Vor)*xy_VelLat, & & - (xy_Coli + xy_Vor)*xy_VelLon)/Rplanet & & - w_Lapla_w(w_xy(xy_VelLon**2 + xy_VelLat**2)*0.5D0)/Rplanet**2)) call EndSub('get_dHsfc') end subroutine get_dHsfc !---------------------------------------------------------------------- ! w_NumVis_w ! 数値粘性項の計算 !---------------------------------------------------------------------- function w_NumVis_w(w_data) real(8) :: w_NumVis_w((nm+1)*(nm+1)) real(8), intent(in) :: w_data((nm+1)*(nm+1)) call BeginSub('w_NumVis_w') w_NumVis_w = w_NumVisScaler_w(w_data) & & - VisCoef*(- (2.0D0/Rplanet**2)**(VisOrder/2.0D0))*w_data call EndSub('w_NumVis_w') end function w_NumVis_w !---------------------------------------------------------------------- ! w_NumVisScaler_w ! スカラーに対する数値粘性項の計算 !---------------------------------------------------------------------- function w_NumVisScaler_w(w_data) real(8) :: w_NumVisScaler_w((nm+1)*(nm+1)) real(8), intent(in) :: w_data((nm+1)*(nm+1)) call BeginSub('w_NumVisScaler_w') w_NumVisScaler_w = & & - VisCoef*( (- rn(:,1)/Rplanet**2)**(VisOrder/2.0D0) )*w_data call EndSub('w_NumVisScaler_w') end function w_NumVisScaler_w !---------------------------------------------------------------------- ! uv ! 流線関数と速度ポテンシャル(スペクトルデータ)から, 水平風(格子データ) ! を計算する. !---------------------------------------------------------------------- subroutine uv(w_Vor, w_HDiv, xy_VelLon, xy_VelLat) real(8),intent(in) :: w_Vor((nm+1)*(nm+1)) real(8),intent(in) :: w_HDiv((nm+1)*(nm+1)) real(8),intent(out) :: xy_VelLon(im,jm) real(8),intent(out) :: xy_VelLat(im,jm) real(8) :: w_Stream((nm+1)*(nm+1)) ! 流線関数 real(8) :: w_VelPot((nm+1)*(nm+1)) ! 速度ポテンシャル call BeginSub('uv') w_Stream = w_LaplaInv_w(w_Vor) * Rplanet**2 w_VelPot = w_LaplaInv_w(w_HDiv) * Rplanet**2 xy_VelLon = (delta*xy_GradLon_w(w_VelPot) & & - xy_GradLat_w(w_Stream))/Rplanet xy_VelLat = (xy_GradLon_w(w_Stream) & & + delta*xy_GradLat_w(w_VelPot))/Rplanet call EndSub('uv') end subroutine uv !---------------------------------------------------------------------- ! timefilter ! Asselin (1972) の時間フィルターを適用し, 変数を 1 ステップ更新する. !---------------------------------------------------------------------- subroutine timefilter(w_data_a, w_data_n, w_data_b) real(8),intent(inout) :: w_data_a((nm+1)*(nm+1)) real(8),intent(inout) :: w_data_n((nm+1)*(nm+1)) real(8),intent(inout) :: w_data_b((nm+1)*(nm+1)) call BeginSub('timefilter') !----- 時間フィルター ----- w_data_n = w_data_n & & + 0.5D0*TfilCoef*(w_data_b - 2.0D0*w_data_n + w_data_a) !----- 1 ステップ進める ----- w_data_b = w_data_n w_data_n = w_data_a call EndSub('timefilter') end subroutine timefilter !---------------------------------------------------------------------- ! InputRestart_gtool4 ! リスタートファイルの読み込み !---------------------------------------------------------------------- subroutine InputRestart_gtool4(filename, Time, xy_Topo, & & xy_VelLon_n, xy_VelLat_n, xy_Hsfc_n, xy_Vor_n, xy_HDiv_n, & & xy_VelLon_b, xy_VelLat_b, xy_Hsfc_b, xy_Vor_b, xy_HDiv_b ) character(len=*), intent(in) :: filename ! 入力ファイル名 real(8),intent(out) :: Time ! 時刻 real(8),intent(out) :: xy_Topo(im,jm) ! 地形 real(8),intent(out) :: xy_VelLon_n(im,jm) ! 速度経度成分 real(8),intent(out) :: xy_VelLat_n(im,jm) ! 速度緯度成分 real(8),intent(out) :: xy_Hsfc_n(im,jm) ! 変位 real(8),intent(out) :: xy_Vor_n(im,jm) ! 渦度 real(8),intent(out) :: xy_HDiv_n(im,jm) ! 発散 real(8),intent(out) :: xy_VelLon_b(im,jm) ! 速度経度成分 real(8),intent(out) :: xy_VelLat_b(im,jm) ! 速度緯度成分 real(8),intent(out) :: xy_Hsfc_b(im,jm) ! 変位 real(8),intent(out) :: xy_Vor_b(im,jm) ! 渦度 real(8),intent(out) :: xy_HDiv_b(im,jm) ! 発散 call HistoryGet(filename, 't', Time) call HistoryGet(filename, 'topo', xy_Topo) call HistoryGet(filename, 'u_n' , xy_VelLon_n, time=Time) call HistoryGet(filename, 'v_n' , xy_VelLat_n, time=Time) call HistoryGet(filename, 'h_n' , xy_Hsfc_n , time=Time) call HistoryGet(filename, 'vor_n', xy_Vor_n , time=Time) call HistoryGet(filename, 'div_n', xy_HDiv_n , time=Time) call HistoryGet(filename, 'u_b' , xy_VelLon_b, time=Time) call HistoryGet(filename, 'v_b' , xy_VelLat_b, time=Time) call HistoryGet(filename, 'h_b' , xy_Hsfc_b , time=Time) call HistoryGet(filename, 'vor_b', xy_Vor_b , time=Time) call HistoryGet(filename, 'div_b', xy_HDiv_b , time=Time) end subroutine InputRestart_gtool4 !---------------------------------------------------------------------- ! Output_gtool4_init ! 出力ヒストリーファイルの初期化 !---------------------------------------------------------------------- subroutine Output_gtool4_init call HistoryCreate( & ! ヒストリー作成 & file=trim(OutputFile2D), & & title=trim(ExpTitle), & & source=trim(ExpSrc), & & institution=trim(ExpInst), & & dims=(/'lon','lat','t '/), dimsizes=(/im,jm,0/), & & longnames=(/'longitude','latitude ','time '/), & & units=(/'deg.','deg.','sec.'/), & & origin=real(InitTime), interval=real(NStepOut*DelTime),& & history=hst_2d) call HistoryPut('lon',x_Lon*180/pi, hst_2d) ! 変数出力 call HistoryPut('lat',y_Lat*180/pi, hst_2d) ! 変数出力 call HistoryAddVariable( & ! 変数定義 & varname='topo', dims=(/'lon','lat'/), & & longname='Topography', & & units='m', & & xtype='double', & & history=hst_2d) call HistoryPut('topo',xy_Topo, hst_2d) call HistoryAddVariable( & ! 変数定義 & varname='h', dims=(/'lon','lat','t '/), & & longname='height', & & units='m', xtype='double', history=hst_2d) call HistoryAddVariable( & ! 変数定義 & varname='u', dims=(/'lon','lat','t '/), & & longname='velocity(longitude)', & & units='m/s', xtype='double', history=hst_2d) call HistoryAddVariable( & ! 変数定義 & varname='v', dims=(/'lon','lat','t '/), & & longname='velocity(latitude)', & & units='m/s', xtype='double', history=hst_2d) call HistoryAddVariable( & ! 変数定義 & varname='vor', dims=(/'lon','lat','t '/), & & longname='vorticity', & & units='1/s', xtype='double', history=hst_2d) call HistoryAddVariable( & ! 変数定義 & varname='div', dims=(/'lon','lat','t '/), & & longname='divergence', & & units='1/s', xtype='double', history=hst_2d) call HistoryCreate( & ! ヒストリー作成 & file=trim(OutputFile1D), & & title=trim(ExpTitle), & & source=trim(ExpSrc), & & institution=trim(ExpInst), & & dims=(/'t '/), dimsizes=(/0/), & & longnames=(/'time '/), & & units=(/'sec.'/), & & origin=real(InitTime), interval=real(NStepOut*DelTime),& & history=hst_1d) call HistoryAddVariable( & ! 変数定義 & varname='mass', dims=(/'t '/), & & longname='total mass', & & units='m^3', xtype='double', history=hst_1d) call HistoryAddVariable( & ! 変数定義 & varname='energy_all', dims=(/'t '/), & & longname='total energy', & & units='m^5/s^2', xtype='double', history=hst_1d) call HistoryAddVariable( & ! 変数定義 & varname='angmoment_all', dims=(/'t '/), & & longname='total angular momentum', & & units='m^5/s', xtype='double', history=hst_1d) call HistoryAddVariable( & ! 変数定義 & varname='potenst_all', dims=(/'t '/), & & longname='total potential enstrophy', & & units='m/s^-2', xtype='double', history=hst_1d) call HistoryAddVariable( & ! 変数定義 & varname='vor_all', dims=(/'t '/), & & longname='total vorticity', & & units='m^2/s', xtype='double', history=hst_1d) call HistoryAddVariable( & ! 変数定義 & varname='hdiv_all', dims=(/'t '/), & & longname='total vorticity', & & units='m^2/s', xtype='double', history=hst_1d) call HistoryAddVariable( & ! 変数定義 & varname='mass_error', dims=(/'t '/), & & longname='normalized mass error', & & units=' ', xtype='double', history=hst_1d) call HistoryAddVariable( & ! 変数定義 & varname='energy_error', dims=(/'t '/), & & longname='normalized energy error', & & units=' ', xtype='double', history=hst_1d) call HistoryAddVariable( & ! 変数定義 & varname='potenst_error', dims=(/'t '/), & & longname='normalized potential enstrophy error', & & units=' ', xtype='double', history=hst_1d) end subroutine Output_gtool4_init !---------------------------------------------------------------------- ! Output_gtool4 ! ヒストリーファイルへの出力 ! Hsfc = Hsfc + HsfcAvr とする. !---------------------------------------------------------------------- subroutine Output_gtool4 call HistoryPut('t' ,Time ,hst_2d) call HistoryPut('u' ,xy_VelLon_n ,hst_2d) call HistoryPut('v' ,xy_VelLat_n ,hst_2d) call HistoryPut('h' ,xy_Hsfc_n+HsfcAvr,hst_2d) call HistoryPut('vor' ,xy_Vor_n ,hst_2d) call HistoryPut('div' ,xy_HDiv_n ,hst_2d) call HistoryPut('t' ,Time ,hst_1d) call HistoryPut('mass' ,MassAll ,hst_1d) call HistoryPut('energy_all' ,EnergyAll ,hst_1d) call HistoryPut('angmoment_all' ,AngMomentAll,hst_1d) call HistoryPut('potenst_all' ,PotEnstALL ,hst_1d) call HistoryPut('vor_all' ,VorAll ,hst_1d) call HistoryPut('hdiv_all' ,HdivAll ,hst_1d) call HistoryPut('mass_error' ,MassAllError ,hst_1d) call HistoryPut('energy_error',EnergyAllError ,hst_1d) call HistoryPut('potenst_error',PotEnstAllError,hst_1d) end subroutine Output_gtool4 !---------------------------------------------------------------------- ! Output_gtool4_close ! ファイル入出力後処理 !---------------------------------------------------------------------- subroutine Output_gtool4_close call HistoryClose(history=hst_1d) call HistoryClose(history=hst_2d) end subroutine Output_gtool4_close !---------------------------------------------------------------------- ! OutputRestart_gtool4_init ! 出力ファイルの初期化を行う ! 出力変数は xy_VelLon, xy_VelLat, xy_Hsfc, xy_Vor, xy_HDiv. ! リープフロッグスキームを用いているため, 2 ステップ分の変数定義 ! を行う. !---------------------------------------------------------------------- subroutine OutputRestart_gtool4_init call HistoryCreate( & ! ヒストリー作成 & file=trim(OutputRstFile), & & title=trim(ExpTitle), & & source=trim(ExpSrc), & & institution=trim(ExpInst), & & dims=(/'lon','lat','t '/), dimsizes=(/im,jm,0/), & & longnames=(/'longitude','latitude ','time '/), & & units=(/'deg.','deg.','sec.'/), & & origin=real(Time), interval=real(0.0), & & history=hst_rst) call HistoryPut('lon', x_Lon*180/pi, hst_rst) ! 変数出力 call HistoryPut('lat', y_Lat*180/pi, hst_rst) ! 変数出力 call HistoryAddVariable( & ! 変数定義 & varname='topo', dims=(/'lon','lat'/), & & longname='Topography', & & units='m', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 & varname='h_n', dims=(/'lon','lat','t '/), & & longname='height', & & units='m', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 & varname='u_n', dims=(/'lon','lat','t '/), & & longname='velocity(longitude)', & & units='m/s', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 & varname='v_n', dims=(/'lon','lat','t '/), & & longname='velocity(latitude)', & & units='m/s', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 & varname='vor_n', dims=(/'lon','lat','t '/), & & longname='vorticity', & & units='1/s', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 & varname='div_n', dims=(/'lon','lat','t '/), & & longname='divergence', & & units='1/s', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 & varname='h_b', dims=(/'lon','lat','t '/), & & longname='height', & & units='m', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 & varname='u_b', dims=(/'lon','lat','t '/), & & longname='velocity(longitude)', & & units='m/s', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 & varname='v_b', dims=(/'lon','lat','t '/), & & longname='velocity(latitude)', & units='m/s', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 & varname='vor_b', dims=(/'lon','lat','t '/), & & longname='vorticity', & & units='1/s', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 & varname='div_b', dims=(/'lon','lat','t '/), & & longname='divergence', & & units='1/s', xtype='double', history=hst_rst) end subroutine OutputRestart_gtool4_init !---------------------------------------------------------------------- ! output_gtool4 ! 初期値の出力を行う. リープフロッグスキームを用いているため, ! 2 ステップ分のデータを出力する. ! Hsfc = Hsfc + HsfcAvr とする. !---------------------------------------------------------------------- subroutine OutputRestart_gtool4 call HistoryPut('t' ,Time ,hst_rst) call HistoryPut('topo' ,xy_Topo ,hst_rst) call HistoryPut('u_n' ,xy_VelLon_n ,hst_rst) call HistoryPut('v_n' ,xy_VelLat_n ,hst_rst) call HistoryPut('h_n' ,xy_Hsfc_n + HsfcAvr,hst_rst) call HistoryPut('vor_n',xy_Vor_n ,hst_rst) call HistoryPut('div_n',xy_HDiv_n ,hst_rst) call HistoryPut('u_b' ,xy_VelLon_b ,hst_rst) call HistoryPut('v_b' ,xy_VelLat_b ,hst_rst) call HistoryPut('h_b' ,xy_Hsfc_b + HsfcAvr,hst_rst) call HistoryPut('vor_b',xy_Vor_b ,hst_rst) call HistoryPut('div_b',xy_HDiv_b ,hst_rst) end subroutine OutputRestart_gtool4 !---------------------------------------------------------------------- ! OutputRestart_gtool4_close ! ファイル入出力後処理 !---------------------------------------------------------------------- subroutine OutputRestart_gtool4_close call HistoryClose(history=hst_rst) end subroutine OutputRestart_gtool4_close end program shallow_zd_bench06