!-------------------------------------------------------------------------- ! Copyright (C) 2005 SPMODEL Development Group. All rights reserved. !-------------------------------------------------------------------------- ! !Title(表題): Two-dimensional barotropic model on a rotating sphere. ! (回転球面上の2次元順圧モデル) ! ! The time integration is performed with ! the 4th order Runge-Kutta scheme. ! The viscous linear term ant the beta term are transformed ! by the exp., cos, sin functions. ! Program is organized for decaying tubulence problem. ! !History(履歴): 2005/10/02 S. Takehiro, created ! 2005/10/05 S. Takehiro, m_Esp, m_EnsSp added. ! 2008/08/10 S. Takehiro, index changed, im -> 0:im-1 ! ! The governing equation(支配方程式): ! ! d\nabla^2\psi/dt ! = -J(\psi,\nabla^2\psi)/a^2 - 2\Omega/a^2 d\psi/d\lambda ! +(-1)^{p+1}\nu_{2p}(\nabla^2 + 2/a^2)^p \nabla^2\psi ! program spbaro_rn4expcs_zonalflow !== モジュール引用宣言 ================================================ use w_module use gt4_history, only : GT_HISTORY, HistoryCreate, HistoryPut, HistoryGet, & HistoryAddVariable, HistoryClose, HistoryAddAttr use dc_trace, only : SetDebug, BeginSub, EndSub, DbgMessage use dc_message, only : MessageNotify implicit none !== 宣言部 ============================================================ !---- 変数(格子点データ) ---- real(8), allocatable :: xy_VelLon(:,:) ! 速度経度成分 real(8), allocatable :: xy_VelLat(:,:) ! 速度緯度成分 real(8), allocatable :: xy_Vor(:,:) ! 渦度(鉛直成分) real(8), allocatable :: xy_StrFunc(:,:) ! 流線関数 !---- 変数(スペクトルデータ) ---- real(8), allocatable :: w_Vor(:) ! 渦度(鉛直成分) real(8), allocatable :: w_StrFunc(:) ! 流線関数 real(8), allocatable :: w_Xi(:) ! ξ : ζ=ξexp[-(iω+ν)t] real(8), allocatable :: w_Vortmp(:) ! 渦度(ルンゲクッタ計算作業用) real(8), allocatable :: w_k1(:) ! 拡散以外の渦度時間変化(1段目) real(8), allocatable :: w_k2(:) ! 拡散以外の渦度時間変化(2段目) real(8), allocatable :: w_k3(:) ! 拡散以外の渦度時間変化(3段目) real(8), allocatable :: w_k4(:) ! 拡散以外の渦度時間変化(4段目) real(8), allocatable :: nm_ESp(:,:) ! エネルギースペクトル real(8), allocatable :: nm_EnsSp(:,:) ! エンストロフィースペクトル real(8), allocatable :: n_ESp(:) ! エネルギースペクトル real(8), allocatable :: m_ESp(:) ! エネルギースペクトル real(8), allocatable :: n_EnsSp(:) ! エンストロフィースペクトル real(8), allocatable :: m_EnsSp(:) ! エンストロフィースペクトル !---- 固定パラメタ ----- real(8), parameter :: pi = 3.141592653589793D0 ! 円周率 character(len=20) :: DbgMessageFmt='*** DbgMESSAGE ***' real(8), parameter :: vmiss = -999.0 ! 欠損値 !---- NAMELIST 変数 ---- ! NAMELIST 入力用テンポラリファイル名(名前変更禁止) character(len=30),parameter :: nmlfile='spbaro_rn4expcs_zonalflow.nml' logical :: Verbose=.false. ! 冗長なメッセージ出力 On/Off logical :: DebugOn=.false. ! デバッグメッセージ On/Off 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=1.0D0 ! 球の半径 real(8) :: Omega=0.0D0 ! 回転角速度 integer :: HVOrder=1 ! 超粘性の次数(1 で普通の粘性, ! 水平ラプラシアンの階数) real(8) :: HVisc=1.0D0 ! 超粘性係数 namelist /physics/ Radius, Omega, HVOrder, HVisc ! -- 初期値 -- character(len=100) :: initial_file='' ! 初期値データファイル名 ! (空なら内部で初期値を計算) real :: initial_time=0.0 ! 初期時刻 namelist /initial/ initial_file, initial_time ! ! -- 時間積分 -- real(8) :: delta_t=1.0e-7 ! 時間積分刻み integer :: nstep=2000 ! 時間積分ステップ数 namelist /tint/ delta_t, nstep ! ! -- ヒストリー出力 -- character(len=100) :: hst_file= '' ! ヒストリーファイル名 character(len=100) :: title = & ! タイトル 'Stability of zonal flows of 2-dim barotropic 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 !---- 作業変数 ---- real(8) :: Beta=0.0d0 ! βパラメター(2*Omega/radius**2) real(8), allocatable :: w_HVisc(:) ! 超粘性係数 real(8), allocatable :: n_Damptime(:) ! 減衰時間(e-folding time) real(8), allocatable :: w_OmegaR(:) ! ロスビー波の振動数 ! (- Beta * abs(m)/(n*(n+1)) integer :: it=0 ! 時間ステップ real(8) :: time ! モデル内時間 integer :: n, m ! 全波数, 帯状波数 type(GT_HISTORY) :: hst_rst ! リスタート GT_HISTORY 変数 !---------------- 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=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_VelLon(0:im-1,jm),xy_VelLat(0:im-1,jm)) allocate(xy_Vor(0:im-1,jm),xy_StrFunc(0:im-1,jm)) allocate(w_Vor((nm+1)*(nm+1)),w_StrFunc((nm+1)*(nm+1))) allocate(w_Xi((nm+1)*(nm+1))) allocate(w_Vortmp((nm+1)*(nm+1))) allocate(w_k1((nm+1)*(nm+1)),w_k2((nm+1)*(nm+1))) allocate(w_k3((nm+1)*(nm+1)),w_k4((nm+1)*(nm+1))) allocate(w_HVisc((nm+1)*(nm+1))) allocate(w_OmegaR((nm+1)*(nm+1))) allocate(nm_ESp(0:nm,-nm:nm),nm_EnsSp(0:nm,-nm:nm)) allocate(m_ESp(-nm:nm),m_EnsSp(-nm:nm)) allocate(n_ESp(0:nm),n_EnsSp(0:nm)) allocate(n_Damptime(0:nm)) !------------------ 座標値の設定 ----------------------- call DbgMessage(fmt='call %c', c1='w_initial') call w_Initial(nm,im,jm) w_spectrum_VMiss = vmiss !------------------ 物理係数の設定 ----------------------- Beta = 2 *Omega/Radius**2 w_HVisc = (-1)**(HVOrder)*HVisc*((rn(:,1)+2.0D0)/Radius**2)**HVOrder w_HVisc(l_nm(0,0)) = vmiss ! rn(ln(0,0,1) は正の値なので修正しておく. n_Damptime(0)=vmiss ! 欠損値 do n=1,nm if ( w_HVisc(l_nm(n,0)) .NE. 0.0D0 ) then n_Damptime(n)=1.0D0/w_HVisc(l_nm(n,0)) else n_Damptime(n)=vmiss ! 欠損値 endif enddo w_OmegaR = vmiss do n=1,nm do m=-n,n w_OmegaR(l_nm(n,m)) = - Beta * abs(m)/(n*(n+1)) enddo enddo !------------------- 初期値設定 ---------------------- time = initial_time if ( initial_file == "") then call set_initial_values ! xy_Vor, xy_Strfunc を定める else call HistoryGet( trim(initial_file), 'w_vor', w_Vor, time ) endif w_Strfunc = w_LaplaInv_w(w_Vor) * Radius**2 !------------------- 時間積分(Euler 法) -------------------- 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.') do it=1,nstep time = initial_time + it * delta_t !---- 1 段目 [ k1 = f( x_n, t_n) ] ---- w_k1 = w_DXiDt_w_w(w_Vor,w_Strfunc,0.0D0) !---- 2 段目 [ k2 = f( x_n+k1*dt/2, t_n+dt/2 ) ] ---- w_Xi = w_Vor + w_k1*delta_t/2.0D0 w_Vortmp = w_Xi2Vor_w(w_Xi,delta_t/2.0D0) w_StrFunc = w_LaplaInv_w(w_Vortmp) * Radius**2 w_k2 = w_DXiDt_w_w(w_Vortmp,w_Strfunc,delta_t/2.0D0) !---- 3 段目 [ k3 = f( x_n+k2*dt/2, t_n+dt/2 ) ] ---- w_Xi = w_Vor + w_k2*delta_t/2.0D0 w_Vortmp = w_Xi2Vor_w(w_Xi,delta_t/2.0D0) w_StrFunc = w_LaplaInv_w(w_Vortmp) * Radius**2 w_k3 = w_DXiDt_w_w(w_Vortmp,w_Strfunc,delta_t/2.0D0) !---- 4 段目 [ k4 = f( x_n+k3*dt, t_n+dt ) ] ---- w_Xi = w_Vor + w_k3*delta_t w_Vortmp = w_Xi2Vor_w(w_Xi,delta_t) w_StrFunc = w_LaplaInv_w(w_Vortmp) * Radius**2 w_k4 = w_DXiDt_w_w(w_Vortmp,w_Strfunc,delta_t) !---- 積分 ---- w_Xi =w_Vor + delta_t * ( w_k1/6.0D0 + w_k2/3.0D0 & + w_k3/3.0D0 + w_k4/6.0D0 ) w_Vor = w_Xi2Vor_w(w_Xi,delta_t) w_StrFunc = w_LaplaInv_w(w_Vor) * Radius**2 if(mod(it,hst_intstep) .eq. 0)then ! ヒストリー出力 call output_history endif if(mod(it,rst_intstep) .eq. 0)then ! リスタート出力 call output_restart endif enddo 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 !------------------- 時間変化項(非線形項のみ) ---------------------- function w_DtDVorNonlinear_w_w(w_Vor,w_Strfunc) real(8) :: w_Vor(:) ! 渦度(鉛直成分) real(8) :: w_StrFunc(:) ! 流線関数 real(8) :: w_DtDVorNonlinear_w_w(size(w_Vor)) ! 渦度時間変化(非線形項のみ) w_DtDVorNonlinear_w_w = - w_Jacobian_w_w(w_StrFunc,w_Vor)/Radius**2 end function w_DtDVorNonlinear_w_w !------------------- ξ から ζ を求める ---------------------- function w_Xi2Vor_w(w_Xi,dt) real(8), intent(IN) :: w_Xi(:) ! ξ : ζ=ξexp[-(iω+ν)t] real(8), intent(IN) :: dt ! 時間刻 real(8) :: w_Xi2Vor_w(size(w_Xi)) ! 渦度ζ w_Xi2Vor_w(l_nm(0,0)) = 0.0D0 do n=1,nm do m=0,n w_Xi2Vor_w(l_nm(n,m)) & = w_Xi(l_nm(n,m)) & * exp(-w_HVisc(l_nm(n,m))*dt)*cos(w_OmegaR(l_nm(n,m))*dt) & + w_Xi(l_nm(n,-m)) & * exp(-w_HVisc(l_nm(n,m))*dt)*sin(w_OmegaR(l_nm(n,m))*dt) enddo do m=1,n w_Xi2Vor_w(l_nm(n,-m)) & = w_Xi(l_nm(n,-m)) & * exp(-w_HVisc(l_nm(n,-m))*dt)*cos(w_OmegaR(l_nm(n,-m))*dt) & - w_Xi(l_nm(n,m)) & * exp(-w_HVisc(l_nm(n,-m))*dt)*sin(w_OmegaR(l_nm(n,-m))*dt) enddo enddo end function w_Xi2Vor_w !---------------- ξ(ζ=ξexp[-(iω+ν)t])の時間変化 ------------------- function w_DXiDt_w_w(w_Vor,w_Strfunc,dt) real(8), intent(IN) :: w_Vor(:) ! 渦度ζ real(8), intent(IN) :: w_Strfunc(:) ! 流線関数 real(8), intent(IN) :: dt ! 時間刻 real(8) :: w_DXiDt_w_w(size(w_Vor)) ! ξ時間変化 real(8) :: w_G(size(w_Vor)) ! 作業領域 integer n,m w_G = w_DtDVorNonlinear_w_w(w_Vor,w_Strfunc) w_DXiDt_w_w(l_nm(0,0)) = 0.0D0 do n=1,nm do m=0,n w_DXiDt_w_w(l_nm(n,m)) & = exp(w_HVisc(l_nm(n,m))*dt)*cos(w_OmegaR(l_nm(n,m))*dt) & * w_G(l_nm(n,m)) & - exp(w_HVisc(l_nm(n,m))*dt) * sin(w_OmegaR(l_nm(n,m))*dt) & * w_G(l_nm(n,-m)) enddo do m=1,n w_DXiDt_w_w(l_nm(n,-m)) & = exp(w_HVisc(l_nm(n,-m))*dt) * cos(w_OmegaR(l_nm(n,-m))*dt) & * w_G(l_nm(n,-m)) & + exp(w_HVisc(l_nm(n,-m))*dt) * sin(w_OmegaR(l_nm(n,-m))*dt) & * w_G(l_nm(n,m)) enddo enddo end function w_DXiDt_w_w !=========================== 初期値設定 ============================ ! ! 初期値設定(リスタートファイルない場合のデフォルト設定) ! subroutine set_initial_values ! w_Vor を定める !---- 初期値帯状流 + ランダムホワイトノイズ integer :: n=4 ! 基本流の水平全波数 real(8) :: vor_amplitude =1.0D0 ! 渦度振幅 real(8) :: noise_amplitude=1.0D-6 ! 微小ホワイトノイズの振幅 integer :: seed=0 ! seed(1)に設定するランダム種の値 namelist /initvor/ n,vor_amplitude, noise_amplitude, seed integer :: random_seed_size ! 乱数の種の長さ integer, allocatable :: seedarray(:) ! 乱数の種 real, allocatable :: harvest(:) ! 乱数の値 if (verbose) write(6,nml=initvor) open(10,file=nmlfile,status='OLD') read(10,nml=initvor) ; write(6,nml=initvor) ; close(10) ! 乱数設定 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(harvest((nm+1)*(nm+1))) call random_number(harvest) w_Vor = 0.0d0 w_Vor(l_nm(n,0))=vor_amplitude/sqrt(2.0d0) ! 帯状流設定 w_Vor = w_Vor + harvest * noise_amplitude/sqrt(2.0d0) ! ノイズ設定 w_Vor(l_nm(0,0)) = 0.0d0 ! 意味のない n,m=(0,0) の成分は0 deallocate(seedarray,harvest) end subroutine set_initial_values !=========================== リスタート出力 ============================ ! ! リスタート出力初期化 ! subroutine output_restart_init call HistoryCreate( & file=trim(rst_file), & title=trim(title), & source='spbaro_rn4expcs_zonalflow.f90 (2008/08/10)', & 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 ','1 '/), & 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', xtype='double', history=hst_rst) !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ---- call HistoryAddAttr('lon','+Radius', Radius ,hst_rst) call HistoryAddAttr('lon','+delta_t', delta_t ,hst_rst) call HistoryAddAttr('lon','+Omega', Omega ,hst_rst) call HistoryAddAttr('lon','+HVOrder', HVOrder ,hst_rst) call HistoryAddAttr('lon','+HVisc', HVisc ,hst_rst) call HistoryAddAttr('lon','+delta_t', delta_t ,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) 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='spbaro_rn4expcs_zonalflow.f90 (2008/08/10)', & 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 ','1 '/), & 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='hvisc', dims=(/'n'/), & longname='hyper diffusivity', units='1', xtype='double') call HistoryAddAttr('hvisc','missing_value', vmiss ) call HistoryAddVariable( & ! 変数定義 varname='damptime', dims=(/'n'/), & longname='dissipative damping time (e-folding)', & units='1', xtype='double') call HistoryAddAttr('damptime','missing_value', vmiss ) call HistoryPut('hvisc',w_HVisc(l_nm((/(n,n=0,nm)/),0))) ! 変数出力 call HistoryPut('damptime',n_Damptime) ! 変数出力 !---- 物理変数定義 ---- call HistoryAddVariable( & ! 変数定義 varname='vor', dims=(/'lon','lat','t '/), & longname='Vorticity', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='strfunc', dims=(/'lon','lat','t '/), & longname='Stream function', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='vellon', dims=(/'lon','lat','t '/), & longname='lon-velocity', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='vellat', dims=(/'lon','lat','t '/), & longname='lat-velocity', units='1', xtype='double') !---- 診断量定義 ---- call HistoryAddVariable( & ! 変数定義 varname='ek', dims=(/'t '/), & longname='mean kinetic energy', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='ens', dims=(/'t '/), & longname='mean enstrophy', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='ektot', dims=(/'t '/), & longname='Total kinetic energy', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='enstot', dims=(/'t '/), & longname='Total enstrophy', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='mvlon', dims=(/'lat','t '/), & longname='Mean zonal flow', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='pv', dims=(/'lon','lat','t '/), & longname='Potential vorticity', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='pvgradlat', dims=(/'lon','lat','t '/), & longname='Latitudinal gradient of PV', units='1', xtype='double') !---- エネルギースペクトル診断量定義 ---- call HistoryAddVariable( & ! 変数定義 varname='ekspnm', dims=(/'n','m','t'/), & longname='Energy spectrum (n-m space)', units='1', xtype='double') call HistoryAddAttr('ekspnm','missing_value', vmiss ) call HistoryAddVariable( & ! 変数定義 varname='ekspm', dims=(/'m','t'/), & longname='Energy spectrum (m space)', units='1', xtype='double') call HistoryAddAttr('ekspm','missing_value', vmiss ) call HistoryAddVariable( & ! 変数定義 varname='ekspn', dims=(/'n','t'/), & longname='Energy spectrum', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='eksp', dims=(/'t'/), & longname='Mean energy (by spectrum)', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='eksptot', dims=(/'t'/), & longname='Total energy (by spectrum)', units='1', xtype='double') !---- エンストロフィースペクトル診断量定義 ---- call HistoryAddVariable( & ! 変数定義 varname='ensspnm', dims=(/'n','m','t'/), & longname='Enstrophy spectrum (n-m space)', units='1', xtype='double') call HistoryAddAttr('ensspnm','missing_value', vmiss ) call HistoryAddVariable( & ! 変数定義 varname='ensspn', dims=(/'n','t'/), & longname='Enstrophy spectrum', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='ensspm', dims=(/'m','t'/), & longname='Enstrophy spectrum (m-space)', units='1', xtype='double') call HistoryAddAttr('ensspm','missing_value', vmiss ) call HistoryAddVariable( & ! 変数定義 varname='enssp', dims=(/'t'/), & longname='Mean enstrophy (by spectrum)', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='enssptot', dims=(/'t'/), & longname='Total enstrophy (by spectrum)', units='1', xtype='double') !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ---- call HistoryAddAttr('lon','+Radius', Radius ) call HistoryAddAttr('lon','+delta_t', delta_t ) call HistoryAddAttr('lon','+Omega', Omega ) call HistoryAddAttr('lon','+HVOrder', HVOrder ) call HistoryAddAttr('lon','+HVisc', HVisc ) call HistoryAddAttr('lon','+delta_t', delta_t ) end subroutine output_history_init ! ! ヒストリー出力 ! subroutine output_history write(6,*) ' History file output at it = ',it, ' time = ', time call HistoryPut('t',real(time)) !---- 物理変数出力 ---- xy_Vor = xy_w(w_Vor) xy_StrFunc = xy_w(w_StrFunc) xy_VelLon = -xy_GradLat_w(w_StrFunc)/Radius xy_VelLat = xy_GradLon_w(w_StrFunc)/Radius call HistoryPut('vellon',xy_VelLon) call HistoryPut('vellat',xy_VelLat) call HistoryPut('vor',xy_Vor) call HistoryPut('strfunc',xy_StrFunc) !---- 診断量出力 ---- call HistoryPut('ek', & AvrLonLat_xy((xy_VelLon**2+xy_VelLat**2)/2.0d0)*Radius**2) call HistoryPut('ens',& AvrLonLat_xy(xy_Vor**2/2.0d0)*Radius**2) call HistoryPut('ektot', & IntLonLat_xy((xy_VelLon**2+xy_VelLat**2)/2.0d0)*Radius**2) call HistoryPut('enstot',& IntLonLat_xy(xy_Vor**2/2.0d0)*Radius**2) call HistoryPut('mvlon',y_AvrLon_xy(xy_VelLon)) call HistoryPut('pv',2*Omega*sin(xy_Lat)+xy_Vor) call HistoryPut('pvgradlat',& (2*Omega*cos(xy_Lat)+xy_GradLat_w(w_Vor))/Radius) !---- エネルギースペクトル診断量出力 ---- nm_ESp = nm_EnergyFromStreamfunc_w(w_StrFunc) n_ESp = n_EnergyFromStreamfunc_w(w_StrFunc) do m=-nm,nm m_ESp(m) = sum(nm_ESp(abs(m):nm,m)) end do do m=1,nm m_ESp(m) = m_ESp(m) + m_ESp(-m) ! m >= 0 に集める m_ESp(-m) = vmiss ! m < 0 は欠損値 end do call HistoryPut('ekspnm',nm_ESp) call HistoryPut('ekspn',n_ESp) call HistoryPut('ekspm',m_Esp) call HistoryPut('eksptot',(4*pi)*sum(n_Esp)) call HistoryPut('eksp',sum(n_Esp)/Radius**2) !---- エンストロフィースペクトル診断量出力 ---- nm_EnsSp = nm_EnstrophyFromStreamfunc_w(w_StrFunc) n_EnsSp = n_EnstrophyFromStreamfunc_w(w_StrFunc) do m=-nm,nm m_EnsSp(m) = sum(nm_ESp(abs(m):nm,m)) end do do m=1,nm m_EnsSp(m) = m_EnsSp(m) + m_EnsSp(-m) ! m >= 0 に集める m_EnsSp(-m) = vmiss ! m < 0 は欠損値 end do call HistoryPut('ensspn',n_EnsSp) call HistoryPut('ensspm',m_EnsSp) call HistoryPut('ensspnm',nm_EnsSp) call HistoryPut('enssptot',(4*pi/Radius**2)*sum(n_Enssp)) call HistoryPut('enssp',sum(n_Enssp)) end subroutine output_history ! ! ヒストリー出力終了 ! subroutine output_history_close call HistoryClose end subroutine output_history_close end program spbaro_rn4expcs_zonalflow