!---------------------------------------------------------------------- ! Copyright (C) 2005--2010 SPMODEL Development Group. All rights reserved. !---------------------------------------------------------------------- ! Sample program for gt4f90io and ISPACK 2005/09/28 S.Takehiro ! ! Solving a non-linear 2-D shallow water system on a sphere. ! ! d\zeta/dt = -\Dinv{a(1-\mu^2)}\DP{[(f+\zeta)U)]}{\lamdba} ! -\Dinv{a}\DP{[(f+\zeta)V)]}{\mu} - F_\zeta^{Diff}, ! ! dD/dt = \Dinv{a(1-\mu^2)}\DP{[(f+\zeta)V)]}{\lamdba} ! -\Dinv{a}\DP{[(f+\zeta)U)]}{\mu} ! -\Dlapla[g(h+h_s)+E] - F_D^{Diff}, ! ! dh/dt = -\Dinv{a(1-\mu^2)}\DP{(hU)]}{\lamdba} ! -\Dinv{a}\DP{(hV)}{\mu} - F_h^{Diff}, ! ! \Dlapla\psi = \zeta, \Dlapla\Xi = D, ! ! U = -\frac{(1-\mu^2)}{a}\DP{\psi}{\mu} + \Dinv{a}\DP{\Xi}{\lambda} ! V = \Dinv{a}\DP{\psi}{\lambda} + \frac{(1-\mu^2)}{a}\DP{\Xi}{\mu} ! ! where hs is hight of topography, E is Kinetic energy, (U^2+V^2)/2(1-\mu^2), ! F_*^{Diff} are hyper viscosity ! ! The time integration is performed by using 4th-order Runge-Kutta scheme ! except for the dissipation terms (Crank-Nicolson scheme). ! ! The program is organized for freedecaying turbulence problem. ! ! History: 2005/10/09 S. Takehiro, created ! 2005/10/10 S. Takehiro, longname of div fixed. ! 2005/10/31 S. Nishizawa, define work variables for subroutines as global variable ! 2005/11/02 S. Nishizawa, change arguments of subroutine ! 2008/08/10 S. Takehiro, index changed, im -> 0:im-1 ! 2010/02/07 S. Takehiro, initial perturbation can bu given ! as hsfc, div as well as vor. ! 2015/07/17 T. Kawahara, 渦度強制及びレイリー摩擦を追加 ! program spshallow_zd_rn4cn_freedecay !== モジュール引用宣言 ================================================ use w_module_sjpack use gt4_history, only : HistoryCreate, HistoryPut, HistoryGet, & HistoryAddVariable, HistoryAddAttr, & HistoryClose, GT_HISTORY use dc_trace, only : SetDebug, BeginSub, EndSub, DbgMessage use dc_message, only : MessageNotify use dc_types use dc_string, only : StoA use dc_args implicit none !== 宣言部 ============================================================ !---- 変数(格子点データ)--- real(8), allocatable :: xy_Hsfc(:,:) ! 表面変位擾乱 (t) real(8), allocatable :: xy_Vor(:,:) ! 渦度 (t) real(8), allocatable :: xy_Div(:,:) ! 発散 (t) real(8), allocatable :: xy_Coli(:,:) ! 惑星渦度 real(8), allocatable :: xy_Htopo(:,:) ! 表面地形 real(8), allocatable :: xy_VelLonCosLat(:,:) ! U = xy_VelLon*Cos(xy_Lat) real(8), allocatable :: xy_VelLatCosLat(:,:) ! V = xy_VelLat*Cos(xy_Lat) real(8), allocatable :: xy_VelLon(:,:) ! 速度経度成分(出力用変数) real(8), allocatable :: xy_VelLat(:,:) ! 速度緯度成分(出力用変数) real(8), allocatable :: xy_StrFunc(:,:) ! 流線関数 real(8), allocatable :: xy_VelPot(:,:) ! 速度ポテンシャル real(8), allocatable :: xy_FluxVorLon(:,:) ! (ζ+f)U (副プログラム作業用) real(8), allocatable :: xy_FluxVorLat(:,:) ! (ζ+f)V (副プログラム作業用) ! 20150714 渦度強制用に追加 real(8), allocatable :: xy_DVorDtF(:,:) ! 渦度時間変化(渦度強制) !---- 変数(スペクトルデータ) ---- real(8), allocatable :: w_Hsfc(:) ! 表面変位擾乱 (t) real(8), allocatable :: w_Vor(:) ! 渦度 (t) real(8), allocatable :: w_Div(:) ! 発散 (t) real(8), allocatable :: w_HsfcTmp(:) ! 表面変位擾乱 (Runge-Kutta 用) real(8), allocatable :: w_VorTmp(:) ! 渦度 (Runge-Kutta 用) real(8), allocatable :: w_DivTmp(:) ! 発散 (Runge-Kutta 用) real(8), allocatable :: w_StrFunc(:) ! 流線関数 real(8), allocatable :: w_VelPot(:) ! 速度ポテンシャル real(8), allocatable :: w_Htopo(:) ! 表面地形 real(8), allocatable :: w_KEnergy(:) ! 運動エネルギー real(8), allocatable :: w_DHsfcDtK1(:) ! 変位時間変化(RK 第1段) real(8), allocatable :: w_DHsfcDtK2(:) ! 変位時間変化(RK 第2段) real(8), allocatable :: w_DHsfcDtK3(:) ! 変位時間変化(RK 第3段) real(8), allocatable :: w_DHsfcDtK4(:) ! 変位時間変化(RK 第4段) real(8), allocatable :: w_DVorDtK1(:) ! 渦度時間変化(RK 第1段) real(8), allocatable :: w_DVorDtK2(:) ! 渦度時間変化(RK 第2段) real(8), allocatable :: w_DVorDtK3(:) ! 渦度時間変化(RK 第3段) real(8), allocatable :: w_DVorDtK4(:) ! 渦度時間変化(RK 第4段) ! 20150714 渦度強制用に追加 real(8), allocatable :: w_DVorDtF(:) ! 渦度時間変化(渦度強制) real(8), allocatable :: w_DDivDtK1(:) ! 発散時間変化(RK 第1段) real(8), allocatable :: w_DDivDtK2(:) ! 発散時間変化(RK 第2段) real(8), allocatable :: w_DDivDtK3(:) ! 発散時間変化(RK 第3段) real(8), allocatable :: w_DDivDtK4(:) ! 発散時間変化(RK 第4段) real(8), allocatable :: nm_EvSp(:,:) ! 渦エネルギースペクトル real(8), allocatable :: nm_EdSp(:,:) ! 発散エネルギースペクトル real(8), allocatable :: nm_EkSp(:,:) ! 運動エネルギースペクトル real(8), allocatable :: nm_EpSp(:,:) ! 位置エネルギースペクトル real(8), allocatable :: nm_ESp(:,:) ! 全エネルギースペクトル real(8), allocatable :: nm_EnsSp(:,:) ! エンストロフィースペクトル real(8), allocatable :: n_EvSp(:) ! 渦エネルギースペクトル real(8), allocatable :: n_EdSp(:) ! 発散エネルギースペクトル real(8), allocatable :: n_EkSp(:) ! 運動エネルギースペクトル real(8), allocatable :: n_EpSp(:) ! 位置エネルギースペクトル real(8), allocatable :: n_ESp(:) ! 全エネルギースペクトル real(8), allocatable :: n_EnsSp(:) ! エンストロフィースペクトル !---- 固定パラメタ ----- real(8), parameter :: pi = 3.141592653589793D0 ! 円周率 character(len=20) :: DbgMessageFmt='*** DbgMESSAGE ***' real(8), parameter :: vmiss = -999.0 ! 欠損値 !---- 作業変数 ---- real(8), allocatable :: w_HVisc(:) ! 超粘性係数(運動方程式) real(8), allocatable :: w_HDiff(:) ! 超拡散係数(質量保存式) integer :: it=0 ! 時間ステップ real(8) :: time ! モデル内時間 integer :: n, m ! 全波数, 帯状波数 type(GT_HISTORY) :: hst_rst ! リスタート GT_HISTORY 変数 type(ARGS) :: arg ! コマンド引数解析用 character(STRING), pointer :: argv(:) => null() character(STRING) :: nmlfile ! NAMELIST 設定ファイル名 !--- 乱数生成用 integer :: random_seed_size ! 乱数の種の長さ integer, allocatable :: seedarray(:) ! 乱数の種 real(8), allocatable :: harvest(:) ! 乱数の値 real(8), allocatable :: theta(:) integer :: loopsum ! ループ回数の合計 !---- NAMELIST 変数 ---- ! (メッセージ出力制御) logical :: Verbose=.false. ! 出力メッセージレベル logical :: DebugOn=.false. ! デバッグ出力コントロール namelist /message/ Verbose, DebugOn ! ! (空間解像度) integer :: nm=21 ! 切断全波数 integer :: im=64 ! 経度方向格子点数 ( > 3*nm + 1) integer :: jm=32 ! 緯度方向格子点数 ( > 3*nm/2 ) namelist /gridset/ nm, im, jm ! (物理パラメタ) real(8) :: Radius = 6.37122D6 ! 惑星半径 real(8) :: Omega = 7.292D-5 ! 回転角速度 real(8) :: Alpha = 0.0 ! 回転軸方向(極との角度,deg.) real(8) :: Grav = 9.80616D0 ! 重力加速度 real(8) :: Hbar = 4.0D3 ! 平均水深 integer :: HVOrder = 2 ! 超粘性の次数 ! (水平ラプラシアンの階数) real(8) :: HVisc = 0.5D16 ! 超粘性係数 integer :: HDOrder = 2 ! 超拡散の次数 ! (水平ラプラシアンの階数) real(8) :: HDiff = 0.5D16 ! 超拡散係数 ! 20150715 散逸過程用に追加 real(8) :: TauRD = 1.0D30 ! レイリー摩擦の時定数 real(8) :: TauNC = 0.25D0 ! ニュートン冷却の時定数 ! 20150714 渦度強制用に追加 real(8) :: epsilon_Forcing = 5.0d7 ! 強制のエネルギー注入率 integer :: n_Forcing = 42 ! 強制の中心波数 integer :: delta_n = 4 ! 強制の波数幅 namelist /physics/ Radius, Omega, Alpha, Grav, Hbar, & HVOrder, HVisc, HDOrder, HDiff, & TauRD, TauNC, & epsilon_Forcing, n_Forcing, delta_n ! (初期値) character(len=100) :: initial_file='' ! 初期値データファイル名 real :: initial_time=0.0 ! 初期時刻 namelist /initial/ initial_file, initial_time ! ! (実験パラメター, 内部設定初期値) integer :: seed=0 ! seed(1)に設定する種の値 integer :: nzero=10 ! 初期エネルギースペクトル分布のパラメタ integer :: gamma=100 ! 初期エネルギースペクトル分布のパラメタ real(8) :: Etotal=1.0D0 ! 初期全エネルギーの値 character(len=10) :: disttype='YIHY2002' ! 初期エネルギー分布のタイプ ! (YIHY2002/YY1993) character(len=4) :: fieldtype='Vor' ! 初期場の種類(VOR/DIV/HSFC) namelist /initvalue/ seed, nzero, gamma, Etotal, disttype, fieldtype ! (時間積分) real(8) :: delta_t=1.0e-7 ! 時間積分刻み integer :: nstep=2000 ! 時間積分ステップ数 namelist /tint/ delta_t, nstep ! ! (ヒストリー出力) character(len=100) :: hst_file= '' ! ヒストリーファイル名 character(len=100) :: source = & ! ソースファイル名 'spshallow_zd_rn4cn_freedecay.f90 (2010/02/07)' character(len=100) :: title = & ! タイトル '2-dim shallow water fluid on a rotating sphere' integer :: hst_intstep=200 ! ヒストリー出力間隔ステップ数 namelist /history/ hst_file, title, hst_intstep ! (リスタート出力) character(len=100) :: rst_file='' ! リスタート出力ファイル名 integer :: rst_intstep=200 ! リスタート出力間隔ステップ数 namelist /restart/ rst_file, rst_intstep !== メインルーチン ========================================================== !---------------- コマンドライン解析 --------------------- call MessageNotify('M','main', trim(source)) call Open(arg) call Debug(arg) ; call Help(arg) ; call Strict(arg) call Get(arg, argv) if ( size(argv) .le. 0 ) then call MessageNotify('M','main',& 'Usage: spshallow_zd_rn4cn_freedecay.out [namelist file]') call MessageNotify('E','main','There is no argument. & & Please set the namelist file name after the execution.') else nmlfile=argv(1) call MessageNotify('M','main','Namelist file is '//trim(nmlfile)) endif deallocate(argv) call Close(arg) !---------------- NAMELIST 読み込み --------------------- write(6,nml=message) open(10,file=nmlfile,status='OLD') read(10,nml=message) ; write(6,nml=message) ; close(10) if (verbose) write(6,nml=gridset) open(10,file=nmlfile,status='OLD') read(10,nml=gridset) ; write(6,nml=gridset) ; close(10) if (verbose) write(6,nml=physics) open(10,file=nmlfile,status='OLD') read(10,nml=physics) ; write(6,nml=physics) ; close(10) if (verbose) write(6,nml=initial) open(10,file=nmlfile,status='OLD') read(10,nml=initial) ; write(6,nml=initial) ; close(10) if(verbose)write(6,nml=initvalue) open(10,file=nmlfile,status='OLD') read(10,nml=initvalue) ; write(6,nml=initvalue) ; close(10) if (verbose) write(6,nml=tint) open(10,file=nmlfile,status='OLD') read(10,nml=tint) ; write(6,nml=tint) ; close(10) if (verbose) write(6,nml=history) open(10,file=nmlfile,status='OLD') read(10,nml=history) ; write(6,nml=history) ; close(10) if (verbose) write(6,nml=restart) open(10,file=nmlfile,status='OLD') read(10,nml=restart) ; write(6,nml=restart) ; close(10) !---------------- デバッグ出力制御設定 ----------------- if (DebugOn) then call SetDebug end if !------------------ 変数の割り付け --------------------- allocate(xy_Vor(0:im-1,jm),xy_Div(0:im-1,jm)) allocate(xy_Hsfc(0:im-1,jm),xy_Htopo(0:im-1,jm)) ! 20150714 渦度強制用に追加 allocate(xy_DVorDtF(0:im-1,jm)) allocate(w_Vor((nm+1)*(nm+1)),w_Div((nm+1)*(nm+1))) allocate(w_Hsfc((nm+1)*(nm+1)),w_Htopo((nm+1)*(nm+1))) allocate(w_StrFunc((nm+1)*(nm+1)),w_VelPot((nm+1)*(nm+1))) allocate(w_VorTmp((nm+1)*(nm+1)),w_DivTmp((nm+1)*(nm+1))) allocate(w_HsfcTmp((nm+1)*(nm+1))) allocate(w_KEnergy((nm+1)*(nm+1))) allocate(xy_VelLonCosLat(0:im-1,jm),xy_VelLatCosLat(0:im-1,jm)) allocate(w_DVorDtK1((nm+1)*(nm+1)),w_DVorDtK2((nm+1)*(nm+1)), & w_DVorDtK3((nm+1)*(nm+1)),w_DVorDtK4((nm+1)*(nm+1)) ) allocate(w_DDivDtK1((nm+1)*(nm+1)),w_DDivDtK2((nm+1)*(nm+1)), & w_DDivDtK3((nm+1)*(nm+1)),w_DDivDtK4((nm+1)*(nm+1)) ) allocate(w_DHsfcDtK1((nm+1)*(nm+1)),w_DHsfcDtK2((nm+1)*(nm+1)), & w_DHsfcDtK3((nm+1)*(nm+1)),w_DHsfcDtK4((nm+1)*(nm+1)) ) ! 20150714 渦度強制用に追加 allocate(w_DVorDtF((nm+1)*(nm+1))) allocate(xy_Coli(0:im-1,jm)) allocate(w_HVisc((nm+1)*(nm+1)),w_HDiff((nm+1)*(nm+1))) allocate(nm_EvSp(0:nm,-nm:nm),nm_EdSp(0:nm,-nm:nm)) allocate(nm_EkSp(0:nm,-nm:nm),nm_EpSp(0:nm,-nm:nm)) allocate(nm_ESp(0:nm,-nm:nm),nm_EnsSp(0:nm,-nm:nm)) allocate(n_EvSp(0:nm),n_EdSp(0:nm),n_EkSp(0:nm),n_EpSp(0:nm)) allocate(n_ESp(0:nm),n_EnsSp(0:nm)) allocate(xy_VelLon(0:im-1,jm),xy_VelLat(0:im-1,jm)) allocate(xy_StrFunc(0:im-1,jm),xy_VelPot(0:im-1,jm)) allocate(xy_FluxVorLon(0:im-1,jm),xy_FluxVorLat(0:im-1,jm)) !------------------ 座標値の設定 ----------------------- call DbgMessage(fmt='call %c', c1='w_initial') call w_Initial(nm,im,jm) w_spectrum_VMiss = vmiss !------------------ 物理係数の設定 ----------------------- ! 20150715 散逸過程用に編集 xy_Coli = 2 * Omega * ( -cos(xy_Lon)*cos(xy_Lat)*sin(Alpha*pi/180.0) & + sin(xy_Lat)*cos(Alpha*pi/180.0) ) if (TauRD > 0) then w_HVisc = 1 / TauRD + HVisc & *( (-rn(:,1)/Radius**2)**HVOrder & -(2.0D0/Radius**2)**HVOrder ) else w_HVisc = HVisc & *( (-rn(:,1)/Radius**2)**HVOrder & -(2.0D0/Radius**2)**HVOrder ) endif w_HDiff = 1 / TauNC + HDiff * (-rn(:,1)/Radius**2)**HDOrder ! rn(ln(0,0,1) は正の値なので修正しておく. w_HVisc(l_nm(0,0)) = 0.0D0 !------------------- 初期値設定 ---------------------- time = initial_time if ( initial_file == "") then ! リスタートファイルを指定しない場合, ! 内部で w_Vor, w_Div, w_Hsfc, w_Htopo を与える. if ( fieldtype(1:3) == 'VOR' ) then call set_initial_values_vor else if ( fieldtype(1:3) == 'DIV' ) then call set_initial_values_div else if ( fieldtype(1:4) == 'HSFC' ) then call set_initial_values_hsfc else call MessageNotify('E','main','fieldtype in NML initivaule invalid') endif else ! 初期値設定(リスタートファイルからの読みこみ) call HistoryGet( trim(initial_file), 'w_vor', w_Vor, time ) call HistoryGet( trim(initial_file), 'w_div', w_div, time ) call HistoryGet( trim(initial_file), 'w_hsfc', w_Hsfc, time ) call HistoryGet( trim(initial_file), 'w_htopo', w_Htopo ) endif ! 渦度・発散から流線関数・速度ポテンシャルならびに速度成分を計算 call VorDiv2Velocity( & w_Vor, w_Div ) !(in) 渦度, 発散 xy_Vor = xy_w(w_Vor) ; xy_Div = xy_w(w_Div) ; xy_HSfc = xy_w(w_Hsfc) !------------- 時間積分(Runge-Kutta + Crank-Nicolson 法) -------------- call output_restart_init call output_history_init if ( initial_file == '' ) call output_history ! 内部で与えた初期値は出力 call DbgMessage(fmt='%c %c', & & c1=DbgMessageFmt, & & c2='Time integration starts.') ! 乱数生成 call random_seed(size=random_seed_size) allocate(seedarray(random_seed_size)) call random_seed(get=seedarray) seedarray(1) = seed call random_seed(put=seedarray) loopsum = n_Forcing*(2*delta_n+1) allocate(harvest(1:loopsum), theta(1:loopsum)) do it=1,nstep time = initial_time + it * delta_t !---- 渦度強制 (20150714 追加) ---- call TimeDerivetives_Forcing( & w_DVorDtF) xy_DVorDtF = xy_w(w_DVorDtF) w_Vor = w_Vor + w_DVorDtF * delta_t xy_Vor = xy_w(w_Vor) !---- 1 段目 [ k1 = f( x_n, t_n) ] ---- call TimeDerivetives_noDisppation( & w_Vor, w_Div, w_Hsfc, &!(in) 渦度, 発散, 変位 w_DVorDtK1, w_DDivDtK1, w_DHsfcDtK1) !(out) 時間変化 !---- 2 段目 [ k2 = f( x_n+k1*dt/2, t_n+dt/2 ) ] ---- w_VorTmp = w_Vor + w_DVorDtK1 * delta_t/2.0D0 w_DivTmp = w_Div + w_DDivDtK1 * delta_t/2.0D0 w_HsfcTmp = w_Hsfc + w_DHsfcDtK1 * delta_t/2.0D0 xy_Vor = xy_w(w_VorTmp) xy_Div = xy_w(w_DivTmp) xy_HSfc = xy_w(w_HsfcTmp) call TimeDerivetives_noDisppation( & w_VorTmp, w_DivTmp, w_HsfcTmp, &!(in) 渦度, 発散, 変位 w_DVorDtK2, w_DDivDtK2, w_DHsfcDtK2) !(out) 時間変化 !---- 3 段目 [ k3 = f( x_n+k2*dt/2, t_n+dt/2 ) ] ---- w_VorTmp = w_Vor + w_DVorDtK2 * delta_t/2.0D0 w_DivTmp = w_Div + w_DDivDtK2 * delta_t/2.0D0 w_HsfcTmp = w_Hsfc + w_DHsfcDtK2 * delta_t/2.0D0 xy_Vor = xy_w(w_VorTmp) xy_Div = xy_w(w_DivTmp) xy_HSfc = xy_w(w_HsfcTmp) call TimeDerivetives_noDisppation( & w_VorTmp, w_DivTmp, w_HsfcTmp, &!(in) 渦度, 発散, 変位 w_DVorDtK3, w_DDivDtK3, w_DHsfcDtK3) !(out) 時間変化 !---- 4 段目 [ k4 = f( x_n+k3*dt, t_n+dt ) ] ---- w_VorTmp = w_Vor + w_DVorDtK3 * delta_t w_DivTmp = w_Div + w_DDivDtK3 * delta_t w_HsfcTmp = w_Hsfc + w_DHsfcDtK3 * delta_t xy_Vor = xy_w(w_VorTmp) xy_Div = xy_w(w_DivTmp) xy_HSfc = xy_w(w_HsfcTmp) call TimeDerivetives_noDisppation( & w_VorTmp, w_DivTmp, w_HsfcTmp, &!(in) 渦度, 発散, 変位 w_DVorDtK4, w_DDivDtK4, w_DHsfcDtK4) !(out) 時間変化 !---- 積分(Runge-Kutta) ---- w_Vor = w_Vor + delta_t * ( w_DVorDtK1/6.0D0 + w_DVorDtK2/3.0D0 & + w_DVorDtK3/3.0D0 + w_DVorDtK4/6.0D0 ) w_Div = w_Div + delta_t * ( w_DDivDtK1/6.0D0 + w_DDivDtK2/3.0D0 & + w_DDivDtK3/3.0D0 + w_DDivDtK4/6.0D0 ) w_Hsfc = w_Hsfc + delta_t * ( w_DHsfcDtK1/6.0D0 + w_DHsfcDtK2/3.0D0 & + w_DHsfcDtK3/3.0D0 + w_DHsfcDtK4/6.0D0 ) !---- 散逸項積分(Crank-Nicolson) ---- w_Vor = (1-w_HVisc*delta_t/2.0D0)/(1+w_HVisc*delta_t/2.0D0) * w_Vor w_Div = (1-w_HVisc*delta_t/2.0D0)/(1+w_HVisc*delta_t/2.0D0) * w_Div w_Hsfc = (1-w_HDiff*delta_t/2.0D0)/(1+w_HDiff*delta_t/2.0D0) * w_Hsfc xy_Vor = xy_w(w_Vor) ; xy_Div = xy_w(w_Div) ; xy_HSfc = xy_w(w_Hsfc) if(mod(it,hst_intstep) .eq. 0)then ! ヒストリー出力 ! 渦度・発散から流線関数・速度ポテンシャルならびに速度成分を計算 call VorDiv2Velocity( & w_Vor, w_Div ) !(in) 渦度, 発散 call output_history endif if(mod(it,rst_intstep) .eq. 0)then ! リスタート出力 call output_restart endif enddo deallocate(harvest, seedarray, theta) call DbgMessage(fmt='%c %c', & & c1=DbgMessageFmt, & & c2='Time integration end.') if(.not. mod(it-1,rst_intstep) .eq. 0)then ! 最終出力 call output_restart endif call output_restart_close call output_history_close ! 以上 メインプログラム !----------------------------------------------------------------------------- ! 以下 サブルーチン contains !=========================== 初期値設定 ============================ ! ! 初期値設定(リスタートファイルない場合のデフォルト設定, 渦度擾乱) ! subroutine set_initial_values_vor ! w_Vor 擾乱を与える. !---- 初期値乱数分布 integer :: random_seed_size ! 乱数の種の長さ integer, allocatable :: seedarray(:) ! 乱数の種 real, allocatable :: harvest(:) ! 乱数の値 real(8), allocatable :: n_EspInit(:) ! 初期エネルギースペクトル分布 ! 乱数設定 call random_seed(size=random_seed_size) allocate(seedarray(random_seed_size)) call random_seed(get=seedarray) seedarray(1)=seed call random_seed(put=seedarray) allocate(n_EspInit(0:nm)) ! 初期のエネルギースペクトル分布の形状を与える if ( trim(disttype) .eq. 'YIHY2002' ) then call MessageNotify('M','set_initial_values_vor',& 'Selected initial energy spectrum is "YIHY2002" type. ') n_EspInit = (/(dble(n)**(gamma/2)/dble(n+nzero)**gamma,n=0,nm)/) elseif ( trim(disttype) .eq. 'YY1993' ) then call MessageNotify('M','set_initial_values_vor',& 'Selected initial energy spectrum is "YY1993" type. ') n_EspInit = (/(dble(n)**5/exp(-dble(n)/2),n=0,nm)/) else call MessageNotify('E','set_initial_values_vor',& 'The parameter "disttype" should be "YIHY2002" or "YY1993"') endif w_StrFunc = 0.0d0 do n=1,nm allocate(harvest(-n:n)) call random_number(harvest) w_StrFunc(l_nm(n,(/(m,m=-n,n)/)))=2.0 * harvest - 1 ! [-1,1] の一様乱数 deallocate(harvest) enddo n_ESp = n_EnergyFromStreamfunc_w(w_StrFunc) ! スペクトル分布の形を決める do n=0,nm do m=-n,n if ( n_ESp(n) .ne. 0.0d0 ) then w_StrFunc(l_nm(n,m)) = w_StrFunc(l_nm(n,m)) & * sqrt( n_EspInit(n)/n_ESp(n) ) endif enddo enddo ! エネルギースペクトル再計算 n_ESp = n_EnergyFromStreamfunc_w(w_StrFunc) ! エネルギーの大きさを定める. w_StrFunc = w_StrFunc * sqrt(Etotal/sum(n_ESp)) xy_StrFunc = xy_w(w_StrFunc) w_Vor = w_Lapla_w(w_StrFunc)/Radius**2 xy_Vor = xy_w(w_Vor) ! 渦度場以外は 0 w_Div = 0.0D0 w_Hsfc = 0.0D0 w_Htopo = 0.0D0 deallocate(n_EspInit,seedarray) end subroutine set_initial_values_vor ! ! 初期値設定(リスタートファイルない場合のデフォルト設定, 発散擾乱) ! subroutine set_initial_values_div ! w_Div 擾乱を与える. !---- 初期値乱数分布 integer :: random_seed_size ! 乱数の種の長さ integer, allocatable :: seedarray(:) ! 乱数の種 real, allocatable :: harvest(:) ! 乱数の値 real(8), allocatable :: n_EspInit(:) ! 初期エネルギースペクトル分布 ! 乱数設定 call random_seed(size=random_seed_size) allocate(seedarray(random_seed_size)) call random_seed(get=seedarray) seedarray(1)=seed call random_seed(put=seedarray) allocate(n_EspInit(0:nm)) ! 初期のエネルギースペクトル分布の形状を与える if ( trim(disttype) .eq. 'YIHY2002' ) then call MessageNotify('M','set_initial_values_div',& 'Selected initial energy spectrum is "YIHY2002" type. ') n_EspInit = (/(dble(n)**(gamma/2)/dble(n+nzero)**gamma,n=0,nm)/) elseif ( trim(disttype) .eq. 'YY1993' ) then call MessageNotify('M','set_initial_values_div',& 'Selected initial energy spectrum is "YY1993" type. ') n_EspInit = (/(dble(n)**5/exp(-dble(n)/2),n=0,nm)/) else call MessageNotify('E','set_initial_values_div',& 'The parameter "disttype" should be "YIHY2002" or "YY1993"') endif w_VelPot = 0.0d0 do n=1,nm allocate(harvest(-n:n)) call random_number(harvest) w_VelPot(l_nm(n,(/(m,m=-n,n)/)))=2.0 * harvest - 1 ! [-1,1] の一様乱数 deallocate(harvest) enddo n_ESp = n_EnergyFromStreamfunc_w(w_VelPot) ! スペクトル分布の形を決める do n=0,nm do m=-n,n if ( n_ESp(n) .ne. 0.0d0 ) then w_VelPot(l_nm(n,m)) = w_VelPot(l_nm(n,m)) & * sqrt( n_EspInit(n)/n_ESp(n) ) endif enddo enddo ! エネルギースペクトル再計算 n_ESp = n_EnergyFromStreamfunc_w(w_VelPot) ! エネルギーの大きさを定める. w_VelPot = w_VelPot * sqrt(Etotal/sum(n_ESp)) xy_VelPot = xy_w(w_VelPot) w_Div = w_Lapla_w(w_VelPot)/Radius**2 xy_Div = xy_w(w_Div) ! 渦度場以外は 0 w_Vor = 0.0D0 w_Hsfc = 0.0D0 w_Htopo = 0.0D0 deallocate(n_EspInit,seedarray) end subroutine set_initial_values_div ! ! 初期値設定(リスタートファイルない場合のデフォルト設定, 表面変位場擾乱) ! subroutine set_initial_values_hsfc ! w_Hsfc 擾乱を与える. !---- 初期値乱数分布 integer :: random_seed_size ! 乱数の種の長さ integer, allocatable :: seedarray(:) ! 乱数の種 real, allocatable :: harvest(:) ! 乱数の値 real(8), allocatable :: n_EspInit(:) ! 初期エネルギースペクトル分布 ! 乱数設定 call random_seed(size=random_seed_size) allocate(seedarray(random_seed_size)) call random_seed(get=seedarray) seedarray(1)=seed call random_seed(put=seedarray) allocate(n_EspInit(0:nm)) ! 初期のエネルギースペクトル分布の形状を与える if ( trim(disttype) .eq. 'YIHY2002' ) then call MessageNotify('M','set_initial_values_hsfc',& 'Selected initial energy spectrum is "YIHY2002" type. ') n_EspInit = (/(dble(n)**(gamma/2)/dble(n+nzero)**gamma,n=0,nm)/) elseif ( trim(disttype) .eq. 'YY1993' ) then call MessageNotify('M','set_initial_values_hsfc',& 'Selected initial energy spectrum is "YY1993" type. ') n_EspInit = (/(dble(n)**5/exp(-dble(n)/2),n=0,nm)/) else call MessageNotify('E','set_initial_values_hsfc',& 'The parameter "disttype" should be "YIHY2002" or "YY1993"') endif w_Hsfc = 0.0d0 do n=1,nm allocate(harvest(-n:n)) call random_number(harvest) w_Hsfc(l_nm(n,(/(m,m=-n,n)/)))=2.0 * harvest - 1 ! [-1,1] の一様乱数 deallocate(harvest) enddo n_ESp = n_EnergyFromSurfaceHeight_w(w_Hsfc) ! スペクトル分布の形を決める do n=0,nm do m=-n,n if ( n_ESp(n) .ne. 0.0d0 ) then w_Hsfc(l_nm(n,m)) = w_Hsfc(l_nm(n,m)) & * sqrt( n_EspInit(n)/n_ESp(n) ) endif enddo enddo ! エネルギースペクトル再計算 n_ESp = n_EnergyFromSurfaceHeight_w(w_Hsfc) ! エネルギーの大きさを定める. w_Hsfc = w_Hsfc * sqrt(Etotal/sum(n_ESp)) ! 表面変位場以外は 0 w_Vor = 0.0D0 w_Div = 0.0D0 w_Htopo = 0.0D0 deallocate(n_EspInit,seedarray) end subroutine set_initial_values_hsfc !---------------- ポテンシャルエネルギースペクトル計算 ---------------- function nm_EnergyFromSurfaceHeight_w(w_Hsfc) ! ! 表面変位のスペクトルデータからポテンシャルエネルギーの球面調和函数成分 ! (スペクトル)を計算する(1 層用). ! ! * 全波数 n, 帯状波数 m の表面変位のスペクトル成分Hsfc(n,m) から ! エネルギースペクトルは (1/2)Grav Hsfc(n,m)^2 と計算される. ! Rd は変形半径である ! ! * 全てのエネルギースペクトル成分の和に4πをかけたものが球面上での ! 全エネルギーに等しい. ! ! * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される. ! 欠損値の値はモジュール変数 w_spectrum_VMiss によって設定できる ! (初期値は -999.0) ! real(8), intent(in) :: w_Hsfc(:) !(in) 表面変位(スペクトルデータ) real(8), dimension(0:nm,-nm:nm) :: nm_EnergyFromSurfaceHeight_w !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間) integer n,m ! DO 変数 nm_EnergyFromSurfaceHeight_w = w_spectrum_VMiss do n=0,nm do m=-n,n nm_EnergyFromSurfaceHeight_w(n,m) & = 0.5 * Grav* w_Hsfc(l_nm(n,m))**2 enddo enddo end function nm_EnergyFromSurfaceHeight_w function n_EnergyFromSurfaceHeight_w(w_Hsfc) ! ! 表面変位のスペクトルデータから各全波数のポテンシャルエネルギー成分 ! (スペクトル)を計算する(1 層用). ! ! * 全波数 n の表面変位のスペクトル成分H(n,m) から ! エネルギースペクトルはΣ[m=-nm]^nm(1/2)Grav*H(n,m)^2 ! と計算される. ! ! * 全てのエネルギースペクトル成分の和に 4πをかけたものが ! 球面上での全エネルギーに等しい. ! real(8), intent(in) :: w_Hsfc(:) !(in) 表面変位(スペクトルデータ) real(8), dimension(0:nm) :: n_EnergyFromSurfaceHeight_w !(out) エネルギースペクトル (水平全波数 n 空間) integer n,m ! DO 変数 do n=0,nm n_EnergyFromSurfaceHeight_w(n) & = 0.5 * Grav * sum(w_Hsfc(l_nm(n,(/(m,m=-n,n)/)))**2,1) enddo end function n_EnergyFromSurfaceHeight_w !=========================== 時間変化項計算 ============================ ! ! 渦度矯正による時間変化計算(20150714 追加) ! subroutine TimeDerivetives_Forcing( & w_DVorDtF) real(8), intent(out) :: w_DVorDtF(:) ! 時間変化 integer :: loopnum ! カウンタ ! integer :: loopsum ! ループ回数の合計 ! integer :: random_seed_size ! 乱数の種の長さ ! integer, allocatable :: seedarray(:) ! 乱数の種 ! real(8), allocatable :: harvest(:) ! 乱数の値 ! real(8), allocatable :: theta(:) ! 一度seed を設定すれば同一の計算の中ではランダムに次から次へと出力されるのでseed の設定し直しはいらない! ! なので時間積分の前にseedを設定すれば良い! ! call random_seed(size=random_seed_size) ! allocate(seedarray(random_seed_size)) ! call random_seed(get=seedarray) ! seedarray(1) = seed ! call random_seed(put=seedarray) ! loopsum = n_Forcing*(2*delta_n+1) ! allocate(harvest(1:loopsum), theta(1:loopsum)) call random_number(harvest) theta = 2 * pi * harvest loopnum = 0 do n = n_Forcing - delta_n, n_Forcing + delta_n do m = 1, n loopnum = loopnum+1 w_DVorDtF(l_nm(n,m)) = sqrt(2*n*(n+1)*epsilon_Forcing / ((2*n+1)*delta_n*delta_t))*cos(theta(loopnum)) w_DVorDtF(l_nm(n,-m)) = sqrt(2*n*(n+1)*epsilon_Forcing / ((2*n+1)*delta_n*delta_t))*sin(theta(loopnum)) enddo enddo ! deallocate(harvest, seedarray, theta) end subroutine TimeDerivetives_Forcing ! ! 散逸項を除いた時間変化計算 ! subroutine TimeDerivetives_noDisppation( & w_Vor, w_Div, w_Hsfc, & !(in) 渦度, 発散, 変位 w_DVorDtNoDisp, w_DDivDtNoDisp, w_DHsfcDtNoDisp) !(out) 時間変化 !---- 変数(入力) ---- real(8), intent(in) :: w_Vor(:) ! 渦度 real(8), intent(in) :: w_Div(:) ! 発散 real(8), intent(in) :: w_Hsfc(:) ! 変位 !---- 変数(出力) ---- real(8), intent(out) :: w_DHsfcDtNoDisp(:) ! 変位時間変化 real(8), intent(out) :: w_DVorDtNoDisp(:) ! 渦度時間変化 real(8), intent(out) :: w_DDivDtNoDisp(:) ! 発散時間変化 ! 渦度・発散から流線関数・速度ポテンシャルならびに速度成分を計算 call VorDiv2Velocity( & w_Vor, w_Div ) !(in) 渦度, 発散 xy_FluxVorLon = (xy_Vor + xy_Coli) * xy_VelLonCosLat xy_FluxVorLat = (xy_Vor + xy_Coli) * xy_VelLatCosLat w_KEnergy = w_xy((xy_VelLonCosLat**2 + xy_VelLatCosLat**2) & /(2.0*cos(xy_Lat)**2)) ! 時間変化項の見積り w_DVorDtNoDisp = - w_DivLambda_xy(xy_FluxVorLon)/Radius & - w_DivMu_xy(xy_FluxVorLat)/Radius w_DDivDtNoDisp = w_DivLambda_xy(xy_FluxVorLat)/Radius & - w_DivMu_xy(xy_FluxVorLon)/Radius & - w_Lapla_w( Grav*(w_Hsfc+ w_Htopo)+w_KEnergy )/Radius**2 w_DHsfcDtNoDisp = - w_DivLambda_xy(xy_Hsfc*xy_VelLonCosLat)/Radius & - w_DivMu_xy(xy_Hsfc*xy_VelLatCosLat)/Radius & - Hbar * w_Div end subroutine TimeDerivetives_noDisppation ! ! 渦度・発散から流線関数・速度ポテンシャルならびに速度成分を計算 ! subroutine VorDiv2Velocity( & w_Vor, w_Div ) !(in) 渦度, 発散 !---- 変数(入力) ---- real(8), intent(in) :: w_Vor(:) ! 渦度 real(8), intent(in) :: w_Div(:) ! 発散 ! 渦度・発散から流線関数・速度ポテンシャルを計算 w_Strfunc = w_LaplaInv_w(w_Vor) * Radius**2 w_VelPot = w_LaplaInv_w(w_Div) * Radius**2 ! 流線関数・速度ポテンシャルから速度成分を計算 xy_VelLonCosLat = -xy_Gradmu_w(w_StrFunc)/Radius & + xy_GradLambda_w(w_VelPot)/Radius xy_VelLatCosLat = xy_GradLambda_w(w_StrFunc)/Radius & + xy_Gradmu_w(w_VelPot)/Radius end subroutine VorDiv2Velocity !=========================== リスタート出力 ============================ ! ! リスタート出力初期化 ! subroutine output_restart_init call HistoryCreate( & file=trim(rst_file), title=trim(title), source=trim(source), & institution='GFD_Dennou Club SPMODEL project', & dims=(/'lon','lat','nm ','t '/), & dimsizes=(/im,jm,(nm+1)**2,0/),& longnames=(/'Longitude ','Latitude ',& 'Hor.wave number index','time '/),& units=(/'radian','radian','1 ','sec '/), & origin=real(time), interval=real(rst_intstep*delta_t), & xtypes=(/'real'/), history=hst_rst) !---- 座標変数定義, 出力 ---- call HistoryPut('lon',x_Lon, hst_rst) ! 変数出力 call HistoryAddattr('lon','topology','circular', hst_rst) ! 周期属性 call HistoryAddattr('lon','modulo',2*pi, hst_rst) ! 周期属性 call HistoryPut('lat',y_Lat, hst_rst) ! 変数出力 call HistoryPut('nm',(/(dble(n),n=0,(nm+1)**2)/), hst_rst)! 変数出力 call HistoryAddVariable( & ! 変数定義 varname='lon_weight', dims=(/'lon'/), & longname='weight for integration in longitude', & units='radian', xtype='double',history=hst_rst) call HistoryAddVariable( & ! 変数定義 varname='coslat_lat_weight', dims=(/'lat'/), & longname='cos(lat) weight for integration in latitide', & units='1', xtype='double',history=hst_rst) call HistoryPut('lon_weight',x_Lon_weight,hst_rst) ! 変数出力 call HistoryPut('coslat_lat_weight',y_Lat_weight,hst_rst) ! 変数出力 !---- 物理変数定義 ---- call HistoryAddVariable( & ! 変数定義 varname='w_vor', dims=(/'nm','t '/), & longname='Vorticity', & units='1/sec', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 varname='w_div', dims=(/'nm','t '/), & longname='Divergence', & units='1/sec', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 varname='w_hsfc', dims=(/'nm','t '/), & longname='Suface height', & units='m', xtype='double', history=hst_rst) call HistoryAddVariable( & ! 変数定義 varname='w_htopo', dims=(/'nm'/), & longname='Bottom topography', & units='m', xtype='double', history=hst_rst) !---- 物理変数出力 ---- call HistoryPut('w_htopo', w_Htopo, hst_rst) !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ---- call HistoryAddAttr('lon','+Radius', Radius, hst_rst) call HistoryAddAttr('lon','+Omega', Omega, hst_rst) call HistoryAddAttr('lon','+Alpha', Alpha, hst_rst) call HistoryAddAttr('lon','+Grav', Grav, hst_rst) call HistoryAddAttr('lon','+Hbar', Hbar, hst_rst) call HistoryAddAttr('lon','+HVOrder', HVOrder, hst_rst) call HistoryAddAttr('lon','+HVisc', HVisc, hst_rst) call HistoryAddAttr('lon','+HDOrder', HDOrder, hst_rst) call HistoryAddAttr('lon','+HDiff', HDiff, hst_rst) call HistoryAddAttr('lon','+delta_t', delta_t, hst_rst) call HistoryAddAttr('lon','+TauRD', TauRD, hst_rst) call HistoryAddAttr('lon','+TauNC', TauNC, hst_rst) call HistoryAddAttr('lon','+epsilon_Forcing', epsilon_Forcing, hst_rst) call HistoryAddAttr('lon','+n_Forcing', n_Forcing, hst_rst) call HistoryAddAttr('lon','+delta_n', delta_n, hst_rst) end subroutine output_restart_init ! ! リスタート出力 ! subroutine output_restart write(6,*) ' Restart file output at it = ',it, ' time = ', time call HistoryPut('t',real(time),hst_rst) !---- 物理変数出力 ---- call HistoryPut('w_vor', w_Vor, hst_rst) call HistoryPut('w_div', w_Div, hst_rst) call HistoryPut('w_hsfc',w_Hsfc,hst_rst) end subroutine output_restart ! ! リスタート出力終了 ! subroutine output_restart_close call HistoryClose(hst_rst) end subroutine output_restart_close !=========================== ヒストリー出力 ============================ ! ! ヒストリー出力初期化 ! subroutine output_history_init !---- ヒストリーファイル作成 ---- call HistoryCreate( & file=trim(hst_file), title=trim(title), source=trim(source), & institution='GFD_Dennou Club SPMODEL project', & dims=(/'lon','lat','nm ','n ','m ','t '/), & dimsizes=(/im,jm,(nm+1)**2,nm+1,2*nm+1,0/),& longnames=(/'Longitude ','Latitude ',& 'Hor.wave number index','Hor.total wave number',& 'zonal wave number ','time '/),& units=(/'degree','degree','1 ','1 ','1 ','sec '/), & origin=real(time), interval=real(hst_intstep*delta_t), & xtypes=(/'real'/)) !---- 座標変数定義, 出力 ---- call HistoryPut('lon',x_Lon/pi*180) ! 変数出力 call HistoryAddattr('lon','topology','circular') ! 周期属性 call HistoryAddattr('lon','modulo',360.0) ! 周期属性 call HistoryPut('lat',y_Lat/pi*180) ! 変数出力 call HistoryPut('nm',(/(dble(n),n=0,(nm+1)**2)/)) ! 変数出力 call HistoryPut('n',(/(dble(n),n=0,nm)/)) ! 変数出力 call HistoryPut('m',(/(dble(m),m=-nm,nm)/)) ! 変数出力 call HistoryAddVariable( & ! 変数定義 varname='lon_weight', dims=(/'lon'/), & longname='weight for integration in longitude', & units='radian', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='coslat_lat_weight', dims=(/'lat'/), & longname='cos(lat) weight for integration in latitide', & units='1', xtype='double') call HistoryPut('lon_weight',x_Lon_weight) ! 変数出力 call HistoryPut('coslat_lat_weight',y_Lat_weight) ! 変数出力 !---- 物理変数定義 ---- call HistoryAddVariable( & ! 変数定義 varname='vor', dims=(/'lon','lat','t '/), & longname='Vorticity', units='1/sec', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='div', dims=(/'lon','lat','t '/), & longname='Divergence', units='1/sec', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='hsfc', dims=(/'lon','lat','t '/), & longname='Surface height', units='m', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='strfunc', dims=(/'lon','lat','t '/), & longname='Stream function', units='m2/sec', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='velpot', dims=(/'lon','lat','t '/), & longname='Velocity potential', units='m2/sec', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='vellon', dims=(/'lon','lat','t '/), & longname='lon-velocity', units='m/sec', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='vellat', dims=(/'lon','lat','t '/), & longname='lat-velocity', units='m/sec', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='am', dims=(/'lat','t '/), & longname='Zonal mean angular momentum (U cos\phi)', & units='m^2/sec)', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='qv', dims=(/'lon','lat','t '/), & longname='Potential vorticity', units='1/(m sec)', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='htopo', dims=(/'lon','lat'/), & longname='Bottom topography', units='m', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='ek', dims=(/'t'/), & longname='Kinetic energy', units='m^3/sec^2', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='ep', dims=(/'t'/), & longname='Potential energy', units='m^3/sec^2', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='en', dims=(/'t'/), & longname='Total energy', units='m^3/sec^2', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='F', dims=(/'lon','lat','t '/), & longname='Random Forcing', units='1/sec^2', xtype='double') !---- エネルギースペクトル診断変数定義 ---- call HistoryAddVariable( & ! 変数定義 varname='evspnm', dims=(/'n','m','t'/), & longname='Vortical kinetic energy spectrum (n-m space)', & units='1', xtype='double') call HistoryAddAttr('evspnm','missing_value', vmiss ) call HistoryAddVariable( & ! 変数定義 varname='evspn', dims=(/'n','t'/), & longname='Vortical kinetic energy spectrum', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='evsp', dims=(/'t'/), & longname='Mean vortical kinetic energy (by spectrum)', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='evsptot', dims=(/'t'/), & longname='Total vortical kinetic energy (by spectrum)', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='edspnm', dims=(/'n','m','t'/), & longname='Divergent kinetic energy spectrum (n-m space)', & units='1', xtype='double') call HistoryAddAttr('edspnm','missing_value', vmiss ) call HistoryAddVariable( & ! 変数定義 varname='edspn', dims=(/'n','t'/), & longname='Divergent kinetic energy spectrum', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='edsp', dims=(/'t'/), & longname='Mean divergent kinetic energy (by spectrum)', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='edsptot', dims=(/'t'/), & longname='Total divergent kinetic energy (by spectrum)', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='ekspnm', dims=(/'n','m','t'/), & longname='Kinetic energy spectrum (n-m space)', & units='1', xtype='double') call HistoryAddAttr('ekspnm','missing_value', vmiss ) call HistoryAddVariable( & ! 変数定義 varname='ekspn', dims=(/'n','t'/), & longname='Kinetic energy spectrum', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='eksp', dims=(/'t'/), & longname='Mean kinetic energy (by spectrum)', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='eksptot', dims=(/'t'/), & longname='Total kinetic energy (by spectrum)', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='epspnm', dims=(/'n','m','t'/), & longname='potential energy spectrum (n-m space)', & units='1', xtype='double') call HistoryAddAttr('epspnm','missing_value', vmiss ) call HistoryAddVariable( & ! 変数定義 varname='epspn', dims=(/'n','t'/), & longname='potential energy spectrum', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='epsp', dims=(/'t'/), & longname='Mean potential energy (by spectrum)', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='epsptot', dims=(/'t'/), & longname='Total potential energy (by spectrum)', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='enspnm', dims=(/'n','m','t'/), & longname='Total energy spectrum (n-m space)', & units='1', xtype='double') call HistoryAddAttr('enspnm','missing_value', vmiss ) call HistoryAddVariable( & ! 変数定義 varname='enspn', dims=(/'n','t'/), & longname='Total energy spectrum', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='ensp', dims=(/'t'/), & longname='Mean total energy (by spectrum)', & units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='ensptot', dims=(/'t'/), & longname='Total energy (by spectrum)', & units='1', xtype='double') !---- 物理変数出力 ---- call HistoryPut('htopo', xy_Htopo) !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ---- call HistoryAddAttr('lon','+Radius', Radius ) call HistoryAddAttr('lon','+Omega', Omega ) call HistoryAddAttr('lon','+Alpha', Alpha ) call HistoryAddAttr('lon','+Grav', Grav ) call HistoryAddAttr('lon','+Hbar', Hbar ) call HistoryAddAttr('lon','+HVOrder', HVOrder) call HistoryAddAttr('lon','+HVisc', HVisc ) call HistoryAddAttr('lon','+HDOrder', HDOrder) call HistoryAddAttr('lon','+HDiff', HDiff ) call HistoryAddAttr('lon','+delta_t', delta_t) call HistoryAddAttr('lon','+TauRD', TauRD ) call HistoryAddAttr('lon','+TauNC', TauNC ) call HistoryAddAttr('lon','+epsilon_Forcing', epsilon_Forcing) call HistoryAddAttr('lon','+n_Forcing', n_Forcing) call HistoryAddAttr('lon','+delta_n', delta_n) end subroutine output_history_init ! ! ヒストリー出力 ! subroutine output_history write(6,*) ' History file output at it = ',it, ' time = ', time call HistoryPut('t',real(time)) !---- 物理変数出力 ---- call HistoryPut('vor', xy_Vor) call HistoryPut('div', xy_Div) call HistoryPut('hsfc',xy_Hsfc) call HistoryPut('strfunc', xy_w(w_StrFunc)) call HistoryPut('velpot', xy_w(w_VelPot)) xy_VelLon = xy_VelLonCosLat/cos(xy_Lat) xy_VelLat = xy_VelLatCosLat/cos(xy_Lat) call HistoryPut('vellon',xy_VelLon) call HistoryPut('vellat',xy_VelLat) call HistoryPut('am', y_AvrLon_xy(xy_VelLonCosLat)*Radius) call HistoryPut('qv', (xy_Coli + xy_Vor)/(Hbar + xy_Hsfc)) call HistoryPut('ek',0.5*AvrLonLat_xy((Hbar+xy_hsfc) & *(xy_Vellon**2+xy_Vellat**2))) call HistoryPut('ep',0.5*Grav*AvrLonLat_xy(xy_Hsfc**2)) call HistoryPut('en',0.5*AvrLonLat_xy(& (Hbar+xy_hsfc)*(xy_Vellon**2+xy_Vellat**2)+Grav*xy_Hsfc**2)) ! 20150714 渦度強制用に追加 call HistoryPut('F', xy_DVorDtF) !---- エネルギースペクトル診断出力 ---- nm_EvSp = nm_EnergyFromStreamfunc_w(w_StrFunc) n_EvSp = n_EnergyFromStreamfunc_w(w_StrFunc) call HistoryPut('evspnm',nm_EvSp/Radius**2) call HistoryPut('evspn',n_EvSp/Radius**2) call HistoryPut('evsptot',(4*pi)*sum(n_Evsp)) call HistoryPut('evsp',sum(n_Evsp)/Radius**2) nm_EdSp = nm_EnergyFromStreamfunc_w(w_VelPot) n_EdSp = n_EnergyFromStreamfunc_w(w_VelPot) call HistoryPut('edspnm',nm_EdSp/Radius**2) call HistoryPut('edspn',n_EdSp/Radius**2) call HistoryPut('edsptot',(4*pi)*sum(n_Edsp)) call HistoryPut('edsp',sum(n_Edsp)/Radius**2) nm_EkSp = w_spectrum_VMiss do n=0,nm do m=-n,n nm_EkSp(n,m) = nm_EvSp(n,m) + nm_EdSp(n,m) enddo enddo n_EkSp = n_EvSp + n_EdSp call HistoryPut('ekspnm',nm_EkSp/Radius**2) call HistoryPut('ekspn',n_EkSp/Radius**2) call HistoryPut('eksptot',(4*pi)*sum(n_EkSp)) call HistoryPut('eksp',sum(n_EkSp)/Radius**2) nm_EpSp = nm_EnergyFromSurfaceHeight_w(w_Hsfc) n_EpSp = n_EnergyFromSurfaceHeight_w(w_Hsfc) call HistoryPut('epspnm',nm_EpSp/Radius**2) call HistoryPut('epspn',n_EpSp/Radius**2) call HistoryPut('epsptot',(4*pi)*sum(n_Epsp)) call HistoryPut('epsp',sum(n_Epsp)/Radius**2) nm_ESp = w_spectrum_VMiss do n=0,nm do m=-n,n nm_ESp(n,m) = nm_EkSp(n,m) + nm_EpSp(n,m) enddo enddo n_ESp = n_EkSp + n_EpSp call HistoryPut('enspnm',nm_ESp/Radius**2) call HistoryPut('enspn',n_ESp/Radius**2) call HistoryPut('ensptot',(4*pi)*sum(n_ESp)) call HistoryPut('ensp',sum(n_ESp)/Radius**2) end subroutine output_history ! ! ヒストリー出力終了 ! subroutine output_history_close call HistoryClose end subroutine output_history_close end program spshallow_zd_rn4cn_freedecay