!= Dynamical core module ! ! Authors:: Yasuhiro MORIKAWA, Yukiko YAMADA ! Version:: $Id: dynamics.f90,v 1.8 2006/07/24 07:42:16 morikawa Exp $ ! Tag Name:: $Name: dcpam3-20060725 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2004-2005. All rights reserved. ! License:: See COPYRIGHT[link:../../../COPYRIGHT] ! ! This file provides dynamics_mod module. ! module dynamics_mod ! !== Overview ! ! Calculate Dynamical Core. ! ! 力学コア部分を演算するモジュール。演算している方程式系、 ! および離散化の手法は以下の通り。 ! ! 支配方程式系 : Governing Equations ! * 球座標プリミティブ方程式系 : ! Primitive Equations in Spherical Coordinate ! ! 水平離散化 : Horizontal Discretization ! * スペクトル法 : Spectral Method ! * 三角形切断 : Triangle Truncation ! * 変換法 : Transform Method ! ! 鉛直離散化 : Vertical Discretization ! * σ座標系 : Sigma Coordinate (Arakawa and Suarez(1983)) ! * Lorenz 格子 : Lorenz Grid ! ! 時間積分 : Time Integral ! * leap frog : リープフロッグ ! !== Reference ! ! * Arakawa, A., Suarez, M. J., 1983: ! Vertical differencing of the primitive equations ! in sigma coordinates. ! Mon. Wea. Rev., 111, 34--35. ! !== Error Handling ! !== Known Bugs ! ! * 粘性計算の wa_NumVis_wa と wa_NumVisScaler_wa をとりあえず SPMODEL ! から移植したため、名称などがそのままである。 ! 本来はモジュール内部のものとして名称その他の検討が必要であろう。 ! !== Note ! !== Future Plans ! !現在、力学過程の全てがこのモジュール内にある。本来ならば、 !解く方程式や項の種類によってモジュール化がなされるべきである。 ! use type_mod, only: STRING, REKIND, DBKIND, INTKIND implicit none private public :: dynamics_init, dynamics_leapfrog ! subroutines public :: dynamics_diagnostic, dynamics_end ! subroutines public :: dynamics_diffusion ! subroutines !------------------------------------------------------------------- ! dynamics_init において値を設定するもの !------------------------------------------------------------------- real(DBKIND) :: Kappa ! κ=Ｒ／Ｃｐ (気体定数/定圧比熱) real(DBKIND), pointer, save:: xy_Coli(:,:) =>null() ! 格子点データ(コリオリ力) real(DBKIND), pointer, save:: xy_UVFact(:,:) =>null() ! u→U, v→V のファクター real(DBKIND), pointer, save:: z_DelSigma(:) =>null() ! Δσ(整数レベル) real(DBKIND), pointer, save:: wz_DiffVorDiv(:,:) =>null() ! 運動量 (渦度発散) 水平拡散係数 real(DBKIND), pointer, save:: wz_DiffTherm(:,:) =>null() ! 熱、比湿 水平拡散係数 real(DBKIND), pointer, save:: z_TempAvrXY(:) =>null() ! 平均温度 (整数レベル) real(DBKIND), pointer, save:: r_TempAvrXY(:) =>null() ! 平均温度 (半整数レベル) real(DBKIND), pointer, save:: z_HydroAlpha(:) =>null() ! 静水圧の式の係数 α real(DBKIND), pointer, save:: z_HydroBeta(:) =>null() ! 静水圧の式の係数 β real(DBKIND), pointer, save:: z_TempInpolKappa(:) =>null() ! 温度鉛直補間の係数κ real(DBKIND), pointer, save:: z_TempInpolA(:) =>null() ! 温度鉛直補間の係数ａ real(DBKIND), pointer, save:: z_TempInpolB(:) =>null() ! 温度鉛直補間の係数ｂ ! semi-implicit用行列 real(DBKIND), pointer, save:: zz_DivMtx(:,:) =>null() ! Ｗ:発散の式での温度補正 ! (線形重力波項の効果) real(DBKIND), pointer, save:: zz_TempMtx(:,:) =>null() ! ｈ: ＱＳ (σ移流) - Ｒ real(DBKIND), pointer, save:: zz_TempMtxQ(:,:) =>null() ! Ｑ: dＴ/dσ (線形重力波項の効果) real(DBKIND), pointer, save:: zz_TempMtxS(:,:) =>null() ! Ｓ: dσ/dｔ (線形重力波項の効果) real(DBKIND), pointer, save:: zz_TempMtxR(:,:) =>null() ! Ｒ: ＲＤ = ! κＴ（∂π/∂ｔ＋ｄσ/ｄｔ／σ） ! 気圧変化の効果の係数 ! (線形重力波項の効果) !------------------------------------------------------------------- ! dynamics_leapfrog において利用するもの !------------------------------------------------------------------- real(DBKIND), pointer, save:: xy_lnPs(:,:) =>null() ! π＝ ln Ｐs (t) real(DBKIND), pointer, save:: xy_DlnPsDtN(:,:) =>null() ! πの時間変化 dπ/dt (t) real(DBKIND), pointer, save:: w_DlnPsDtN(:) =>null() ! πの時間変化 (スペクトル) (t) real(DBKIND), pointer, save:: w_lnPsN(:) =>null() ! π (t) real(DBKIND), pointer, save:: xy_GradLonlnPsN(:,:) =>null() ! dπ/dx (t) real(DBKIND), pointer, save:: xy_GradLatlnPsN(:,:) =>null() ! dπ/dy (t) real(DBKIND), pointer, save:: wz_TempN(:,:) => null() ! 温度 (t) real(DBKIND), pointer, save:: w_lnPsA(:) =>null() ! π (t+Δt) real(DBKIND), pointer, save:: xyz_lnPsAdv(:,:,:) =>null() ! π の移流 ｖ・∇π real(DBKIND), pointer, save:: xyz_lnPsAdvSum(:,:,:) =>null() ! π移流積下げ Σ[k〜K](ｖ・∇π)Δσ real(DBKIND), pointer, save:: xyz_DivSum(:,:,:) =>null() ! 発散の積下げ Σ[k〜K] ＤΔσ real(DBKIND), pointer, save:: xyr_VSigma(:,:,:) =>null() ! 鉛直σ速度(半整数レベル) (t) real(DBKIND), pointer, save:: xyr_VSigmaNonG(:,:,:) => null() ! 鉛直σ速度 非重力波成分 ! (半整数レベル) (t) real(DBKIND), pointer, save:: xyz_TempEdd(:,:,:) =>null() ! Ｔ' ：温度の擾乱 (整数レベル) (t) real(DBKIND), pointer, save:: xyr_TempEdd(:,:,:) =>null() ! Ｔ' ：温度の擾乱 (半整数レベル)(t) real(DBKIND), pointer, save:: xyz_TempVir(:,:,:) =>null() ! Ｔv ：仮温度 (t) real(DBKIND), pointer, save:: xyz_TempVirEdd(:,:,:) =>null() ! Ｔv'：仮温度の擾乱 (t) real(DBKIND), pointer, save:: xyz_UAN(:,:,:) =>null() ! 東西運動量変化項ＵＡ real(DBKIND), pointer, save:: xyz_VAN(:,:,:) =>null() ! 南北運動量変化項ＶＡ real(DBKIND), pointer, save:: wz_DVorDtN(:,:) =>null() ! 渦度変化のスペクトルデータ (Δt) real(DBKIND), pointer, save:: wz_VorA(:,:) =>null() ! 渦度のスペクトルデータ (t+Δt) real(DBKIND), pointer, save:: xyz_KE(:,:,:) =>null() ! 運動エネルギー + 仮温度補正 ! ｕ**2+ｖ**2 + ΣＷ(Ｔv-Ｔ) real(DBKIND), pointer, save:: xy_Phi(:,:) =>null() ! 地表ジオポテンシャルΦ real(DBKIND), pointer, save:: wz_DDivDtN(:,:) =>null() ! 発散変化のスペクトルデータ (Δt) real(DBKIND), pointer, save:: wz_DivA(:,:) =>null() ! 発散のスペクトルデータ (t+Δt) real(DBKIND), pointer, save:: wz_PresTendTemp(:,:) =>null() ! 温度による圧力傾度のスペクトルデータ real(DBKIND), pointer, save:: wz_PresTendPs(:,:) =>null() ! 地表圧力による圧力傾度のスペクトルデータ real(DBKIND), pointer, save:: xyz_DTempLocalDtN(:,:,:) => null() ! 温度の局所時間変化項 Ｈ ! (水平発散 Ｔ'Ｄ + σ移流 ! + π変化の効果) (全て非重力波項) real(DBKIND), pointer, save:: wz_DTempDtN(:,:) =>null() ! 温度の変化スペクトルデータ (Δt) real(DBKIND), pointer, save:: wz_TempA(:,:) =>null() ! 温度のスペクトルデータ (t+Δt) real(DBKIND), pointer, save:: xyz_QVapDivVSigmaAdv(:,:,:) => null() ! Ｒ＝ｑＤ＋σ移流 real(DBKIND), pointer, save:: wz_DQVapDtN(:,:) =>null() ! 比湿の変化スペクトルデータ (Δt) real(DBKIND), pointer, save:: wz_QVapA(:,:) =>null() ! 比湿のスペクトルデータ (t+Δt) !------------------------------------------------------------------- ! dynamics_diagnostic で計算する変数 !------------------------------------------------------------------- real(DBKIND), pointer:: wz_PsiA(:,:) =>null() ! スペクトル(流線関数) (t+Δt) real(DBKIND), pointer:: wz_ChiA(:,:) =>null() ! スペクトル(ポテンシャル) (t+Δt) logical, save :: dynamics_initialized = .false. character(STRING),parameter:: version = & & '$Id: dynamics.f90,v 1.8 2006/07/24 07:42:16 morikawa Exp $' character(STRING),parameter:: tagname = '$Name: dcpam3-20060725 $' contains subroutine dynamics_init( & & x_Lon, y_Lat, z_Sigma, r_Sigma & !(in) & ) ! ! Initialize module and Calculate Non-Predictional Value. ! ! 以降のサブルーチンで用いる変数の allocate 、および ! 時間発展しない量の演算を行なう。 ! また、変数データ出力のための初期設定も行なう。 ! use type_mod, only: STRING, REKIND, DBKIND, INTKIND use grid_3d_mod, only: im, jm, km use grid_wavenumber_mod, only: nm use constants_mod, only: constants_init, R0, Omega, Cp, RAir, & & TempAve, VisOrder, EFoldTime use time_mod, only: DelTime use spml_mod, only: spml_init, xy_Lat, rn use io_gt4_out_mod,only: io_gt4_out_init, io_gt4_out_SetVars use dc_trace, only: DbgMessage, BeginSub, EndSub, DataDump use dc_string, only: toChar implicit none real(DBKIND), intent(in) :: x_Lon(:) ! 経度座標 real(DBKIND), intent(in) :: y_Lat(:) ! 緯度座標 real(DBKIND), intent(in) :: z_Sigma(:) ! σレベル(整数)座標 real(DBKIND), intent(in) :: r_Sigma(:) ! σレベル(半整数)座標 !----- 拡散係数計算用変数 ----- real(DBKIND) :: VisCoef ! 超粘性係数 real(DBKIND), pointer :: wz_rn(:,:) =>null() ! ラプラシアンの係数 -n*(n+1) !----- 作業用内部変数 ----- ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用) integer(INTKIND) :: i, j, k, kk, l character(STRING), parameter:: subname = "dynamics_init" continue !---------------------------------------------------------------- ! Check Initialization !---------------------------------------------------------------- call BeginSub(subname) if (dynamics_initialized) then call EndSub( subname, '%c is already called.', c1=trim(subname) ) return else dynamics_initialized = .true. endif !---------------------------------------------------------------- ! Version identifier !---------------------------------------------------------------- call DbgMessage('%c :: %c', c1=trim(version), c2=trim(tagname)) !------------------------------------------------------------------- ! 物理定数の取得 !------------------------------------------------------------------- call constants_init !------------------------------------------------------------------- ! SPMODEL の初期化 !------------------------------------------------------------------- call spml_init !------------------------------------------------------------------- ! io_gt4_out の初期化 !------------------------------------------------------------------- call io_gt4_out_init !------------------------------------------------------------------- ! コリオリ力の設定 !------------------------------------------------------------------- allocate( xy_Coli(im, jm) ) xy_Coli = 2.0 * Omega * sin(xy_Lat) !-- !!$ call DataDump('xy_Coli', xy_Coli, strlen=80) !!$ call DataDump('xy_Lat', xy_Lat, strlen=80) !++ !------------------------------------------------------------------- ! u→U, v→V のファクターの計算 !------------------------------------------------------------------- allocate( xy_UVFact(im, jm) ) xy_UVFact = cos( xy_Lat ) !-- !!$ call DataDump('xy_UVFact', xy_UVFact, strlen=80) !++ !------------------------------------------------------------------- ! Δσ (整数、半整数) の計算 !------------------------------------------------------------------- allocate( z_DelSigma(km) ) do k = 1, km z_DelSigma(k) = r_Sigma(k) - r_Sigma(k+1) enddo !-- !!$ call DataDump('z_Delsigma', z_DelSigma, strlen=80) !++ !------------------------------------------------------------------- ! 運動量 (渦度発散) 拡散係数 wz_DiffVorDiv, ! 熱拡散係数 wz_DiffTherm の計算 !------------------------------------------------------------------- allocate( wz_DiffVorDiv((nm+1)*(nm+1), km) ) allocate( wz_DiffTherm ((nm+1)*(nm+1), km) ) allocate( wz_rn ((nm+1)*(nm+1), km) ) do k = 1, km wz_rn(:,k) = rn(:,1) enddo ! 粘性係数の計算 (最大波数の e-folding time が EFoldTime となるように) VisCoef = (real( (nm * (nm+1)), DBKIND ) / R0**2 )**(-VisOrder / 2) & & / EFoldTime call DbgMessage('VisCoef=<%f>', d=(/VisCoef/)) wz_DiffTherm = - VisCoef * ( (-wz_rn / R0**2)**(VisOrder / 2) ) wz_DiffVorDiv = wz_DiffTherm & & - VisCoef * ( - (2.0d0 / R0**2)**(VisOrder / 2)) !-- !!$ call DataDump('wz_DiffTherm' , wz_DiffTherm, strlen=80) !!$ call DataDump('wz_DiffVorDiv', wz_DiffVorDiv, strlen=80) !++ deallocate(wz_rn) !------------------------------------------------------------------- ! 静水圧の式の係数α、βの計算 !------------------------------------------------------------------- allocate(z_HydroAlpha(km)) allocate(z_HydroBeta(km)) Kappa = RAir / Cp ! κ=Ｒ／Ｃｐ (気体定数/定圧比熱) do k = 1, km z_HydroAlpha(k) = ( r_Sigma(k) / z_Sigma(k) )**Kappa - 1.0 z_HydroBeta(k) = 1.0 - ( r_Sigma(k+1) / z_Sigma(k) )**Kappa enddo !-- !!$ call DbgMessage('Kappa=<%f>', d=(/Kappa/) ) !!$ call DataDump('z_HydroAlpha', z_HydroAlpha, strlen=80) !!$ call DataDump('z_HydroBeta', z_HydroBeta, strlen=80) !++ !------------------------------------------------------------------- ! 温度鉛直補間の係数ａ、ｂの計算 !------------------------------------------------------------------- allocate(z_TempInpolA(km)) allocate(z_TempInpolB(km)) allocate(z_TempInpolKappa(km)) do k = 1, km z_TempInpolA(k) = & & z_HydroAlpha(k) / ( 1.0 - (z_Sigma(k) / z_Sigma(k-1))**Kappa ) z_TempInpolB(k) = & & z_HydroBeta(k) / ( (z_Sigma(k) / z_Sigma(k+1))**Kappa - 1.0 ) z_TempInpolKappa(k) = & & ( r_Sigma(k) * z_HydroAlpha(k) & & + r_Sigma(k+1) * z_HydroBeta(k) ) / z_DelSigma(k) enddo z_TempInpolA(1) = 0.0 z_TempInpolB(km) = 0.0 !-- !!$ call DataDump('z_TempInpolA', z_TempInpolA, strlen=80) !!$ call DataDump('z_TempInpolB', z_TempInpolB, strlen=80) !!$ call DataDump('z_TempInpolKappa', z_TempInpolKappa, strlen=80) !++ !------------------------------------------------------------------- ! 平均温度 (整数レベル、半整数レベル) の計算 !------------------------------------------------------------------- allocate( z_TempAvrXY(km) ) allocate( r_TempAvrXY(km+1) ) z_TempAvrXY = TempAve !----- 下端の平均温度設定 ----- ! ※ 正しいかは微妙だが一応変な値が入らないように r_TempAvrXY = 0.0d0 do k = 2, km r_TempAvrXY(k) = & & z_TempInpolA(k) * z_TempAvrXY(k) & & + z_TempInpolB(k-1) * z_TempAvrXY(k-1) enddo !-- !!$ call DataDump('z_TempAvrXY', z_TempAvrXY, strlen=80) !!$ call DataDump('r_TempAvrXY', r_TempAvrXY, strlen=80) !++ !------------------------------------------------------------------- ! semi-implicit 用行列の計算 !------------------------------------------------------------------- allocate( zz_DivMtx(km,km) ) allocate( zz_TempMtx(km,km) ) allocate( zz_TempMtxQ(km,km) ) allocate( zz_TempMtxS(km,km) ) allocate( zz_TempMtxR(km,km) ) ! Ｗ:発散の式での温度補正 (線形重力波項の効果) ! 初期化 zz_DivMtx = 0.0 do k = 1, km do kk = 1, k zz_DivMtx(k, kk) = Cp * z_HydroAlpha(kk) enddo do kk = 1, k-1 zz_DivMtx(k, kk) = zz_DivMtx(k, kk) + Cp * z_HydroBeta(kk) enddo enddo ! Ｓ: dσ/dｔ (線形重力波項の効果) ! 初期化 zz_TempMtxS = 0.0 do k = 1, km do kk = 1, km zz_TempMtxS(k,kk) = r_Sigma(k) * z_DelSigma(kk) enddo do kk = k, km zz_TempMtxS(k,kk) = zz_TempMtxS(k,kk) - z_DelSigma(kk) enddo enddo !-- !!$ call DataDump('zz_TempMtxS', zz_TempMtxS, strlen=80) !++ ! Ｑ: dＴ/dσ (線形重力波項の効果) ! 初期化 zz_TempMtxQ = 0.0 do k = 1, km zz_TempMtxQ(k,k) = ( r_TempAvrXY(k) - z_TempAvrXY(k) ) / z_DelSigma(k) enddo do k = 1, km - 1 zz_TempMtxQ(k,k+1) = ( z_TempAvrXY(k) - r_TempAvrXY(k+1) ) / z_DelSigma(k) enddo !-- !!$ call DataDump('zz_TempMtxQ', zz_TempMtxQ, strlen=80) !++ ! Ｒ: ＲＤ = κＴ（∂π/∂ｔ＋ｄσ/ｄｔ／σ） ! 気圧変化の効果の係数 (線形重力波項の効果) ! 初期化 zz_TempMtxR = 0.0 do k = 1, km do kk = k, km zz_TempMtxR(k,kk) = & & - z_HydroAlpha(k) / z_DelSigma(k) * z_DelSigma(kk) * z_TempAvrXY(k) enddo do kk = k+1, km zz_TempMtxR(k,kk) = zz_TempMtxR(k,kk) & & - z_HydroBeta(k) / z_DelSigma(k) * z_DelSigma(kk) * z_TempAvrXY(k) enddo enddo !-- !!$ call DataDump('zz_TempMtxR', zz_TempMtxR, strlen=80) !++ ! ｈ: ＱＳ (σ移流) - Ｒ ! 初期化 zz_TempMtx = 0.0 zz_TempMtx = matmul(zz_TempMtxQ, zz_TempMtxS) - zz_TempMtxR !-- !!$ call DataDump('zz_TempMtx', zz_TempMtx, strlen=80) !++ !------------------------------------------------------------------- ! データ出力用の初期設定 !------------------------------------------------------------------- call io_gt4_out_SetVars('xyz_DivSum') call io_gt4_out_SetVars('xyz_lnPsAdv') call io_gt4_out_SetVars('xyz_lnPsAdvSum') call io_gt4_out_SetVars('xyz_UAN') call io_gt4_out_SetVars('xyz_VAN') call io_gt4_out_SetVars('xyz_KE') call io_gt4_out_SetVars('xyz_PresTendTemp') call io_gt4_out_SetVars('xyz_PresTendPs') call io_gt4_out_SetVars('xyr_VSigma') call io_gt4_out_SetVars('xyr_VSigmaNonG') call io_gt4_out_SetVars('xyz_DTempLocalDtN') call io_gt4_out_SetVars('xyz_TempAdv') call io_gt4_out_SetVars('xyr_TempEdd') call io_gt4_out_SetVars('xyz_Temp_T') call io_gt4_out_SetVars('xyz_DTempLocalDtN_lnPsAdvSum') call io_gt4_out_SetVars('TotalMass') call io_gt4_out_SetVars('xyz_Vor_Diffusion') call io_gt4_out_SetVars('xyz_Div_Diffusion') call io_gt4_out_SetVars('xyz_Temp_Diffusion') call io_gt4_out_SetVars('xyz_QVap_Diffusion') !------------------------------------------------------------------- ! dynamics_leapfrog で計算する変数の allocate !------------------------------------------------------------------- allocate & & ( xy_lnPs (im, jm) , & ! π＝ ln Ｐs (t) & xyz_lnPsAdv(im, jm, km) , & ! π の移流 ｖ・∇π & xyz_lnPsAdvSum(im, jm, km), & ! π移流積下げ Σ[k〜K](ｖ・∇π)Δσ & xyz_DivSum(im, jm, km) , & ! 発散の積下げ Σ[k〜K] ＤΔσ & xy_DlnPsDtN( im, jm ) , & ! πの時間変化 dπ/dt (Δt) & w_DlnPsDtN( (nm+1)*(nm+1) ) , & ! πの時間変化 (スペクトル) (Δt) & w_lnPsN( (nm+1)*(nm+1) ) , & ! π (t) & xy_GradLonlnPsN( im, jm ), & ! dπ/dx (t) & xy_GradLatlnPsN( im, jm ), & ! dπ/dy (t) & wz_TempN( (nm+1)*(nm+1), km ) , & ! 温度 (t) & w_lnPsA( (nm+1)*(nm+1) ) , & ! π (t+Δt) & xyr_VSigma(im, jm, km+1) , & ! 鉛直σ速度(半整数レベル) (t) & xyr_VSigmaNonG(im, jm, km+1) , & ! 鉛直σ速度 非重力波成分 (半整数レベル) (t) & xyz_TempEdd(im, jm, km) , & ! Ｔ' ：温度の擾乱 (整数レベル) (t) & xyz_TempVir(im, jm, km) , & ! Ｔv ：仮温度 (virtual temperature)(t) & xyz_TempVirEdd(im, jm, km ), &! Ｔv'：仮温度の擾乱 (t) & xyr_TempEdd(im, jm, km+1), &! Ｔ' ：温度の擾乱 (半整数レベル) (t) & xyz_UAN(im, jm, km) , & ! 東西運動量変化項ＵＡ & xyz_VAN(im, jm, km) , & ! 南北運動量変化項ＶＡ & wz_DVorDtN( (nm+1)*(nm+1), km ) , & ! 渦度変化のスペクトルデータ (Δt) & wz_VorA( (nm+1)*(nm+1), km ) , & ! 渦度のスペクトルデータ (t+Δt) & xyz_KE(im, jm, km) , & ! 運動エネルギー + 仮温度補正 ! ｕ**2+ｖ**2 + ΣＷ(Ｔv-Ｔ) & wz_DDivDtN( (nm+1)*(nm+1), km) , & ! 発散のスペクトルデータ (Δt) & wz_DivA( (nm+1)*(nm+1), km) , & ! 発散のスペクトルデータ (t+Δt) & wz_PresTendTemp( (nm+1)*(nm+1), km) , & ! 温度による圧力傾度のスペクトルデータ & wz_PresTendPs( (nm+1)*(nm+1), km) , & ! 地表圧力による圧力傾度のスペクトルデータ & xyz_DTempLocalDtN(im, jm, km) , & ! 温度の局所時間変化項 Ｈ ! (水平発散 Ｔ'Ｄ + σ移流 ! + π変化の効果) (全て非重力波項) & wz_DTempDtN( (nm+1)*(nm+1), km ) , & ! 温度の変化スペクトルデータ (Δt) & wz_TempA( (nm+1)*(nm+1), km ) , & ! 温度のスペクトルデータ (t+Δt) & xyz_QVapDivVSigmaAdv( im, jm, km), & ! Ｒ＝ｑＤ＋σ移流 & wz_DQVapDtN( (nm+1)*(nm+1), km) , & ! 比湿の変化スペクトルデータ (Δt) & wz_QVapA( (nm+1)*(nm+1), km) & ! 比湿のスペクトルデータ (t+Δt) & ) !------------------------------------------------------------------- ! dynamics_diagnostic で計算する変数の allocate !------------------------------------------------------------------- allocate & & ( wz_PsiA( (nm+1)*(nm+1), km) , & ! スペクトル(流線関数) & wz_ChiA( (nm+1)*(nm+1), km) ) ! スペクトル(ポテンシャル) call EndSub(subname) end subroutine dynamics_init subroutine dynamics_leapfrog & & ( x_Lon , y_Lat , z_Sigma , r_Sigma , & !(in) & xyz_VelLonB, xyz_VelLatB, xyz_VorB, xyz_DivB , & !(in) & xyz_TempB , xyz_QVapB , xy_PsB , & !(in) & xyz_VelLonN, xyz_VelLatN, xyz_VorN, xyz_DivN , & !(in) & xyz_TempN , xyz_QVapN , xy_PsN , & !(in) & xyz_VelLonA, xyz_VelLatA, xyz_VorA, xyz_DivA , & !(out) & xyz_TempA , xyz_QVapA , xy_PsA ) !(out) ! ! Calculate Predictional Values. ! ! 時間発展する量の演算を行なう。 ! 演算するデータの出力も行なう。 ! use type_mod, only: STRING, REKIND, DBKIND, INTKIND use grid_3d_mod, only: im, jm, km use grid_wavenumber_mod, only: nm use constants_mod, only: R0, Cp, EpsVT use time_mod, only: DelTime, CurrentTime use spml_mod, only: w_xy, xy_w , xy_GradLon_w, xy_GradLat_w, & & w_Div_xy_xy, w_LaplaInv_w, & & wa_xya, xya_wa, wa_Div_xya_xya, wa_Lapla_wa use io_gt4_out_mod,only: io_gt4_out_Put use dc_trace, only: DbgMessage, BeginSub, EndSub, DataDump use dc_string, only: toChar implicit none real(DBKIND), intent(in) :: x_Lon(:) ! 経度座標 real(DBKIND), intent(in) :: y_Lat(:) ! 緯度座標 real(DBKIND), intent(in) :: z_Sigma(:) ! σレベル(整数)座標 real(DBKIND), intent(in) :: r_Sigma(:) ! σレベル(半整数)座標 real(DBKIND), intent(in) :: xyz_VelLonB(:,:,:) ! 速度経度成分 (t-Δt) real(DBKIND), intent(in) :: xyz_VelLatB(:,:,:) ! 速度緯度成分 (t-Δt) real(DBKIND), intent(in) :: xyz_VorB(:,:,:) ! 渦度 (t-Δt) real(DBKIND), intent(in) :: xyz_DivB(:,:,:) ! 発散 (t-Δt) real(DBKIND), intent(in) :: xyz_TempB(:,:,:) ! 温度 (t-Δt) real(DBKIND), intent(in) :: xyz_QVapB(:,:,:) ! 比湿 (t-Δt) real(DBKIND), intent(in) :: xy_PsB(:,:) ! 地表面気圧 (t-Δt) real(DBKIND), intent(in) :: xyz_VelLonN(:,:,:) ! 速度経度成分 (t) real(DBKIND), intent(in) :: xyz_VelLatN(:,:,:) ! 速度緯度成分 (t) real(DBKIND), intent(in) :: xyz_VorN(:,:,:) ! 渦度 (t) real(DBKIND), intent(in) :: xyz_DivN(:,:,:) ! 発散 (t) real(DBKIND), intent(in) :: xyz_TempN(:,:,:) ! 温度 (t) real(DBKIND), intent(in) :: xyz_QVapN(:,:,:) ! 比湿 (t) real(DBKIND), intent(in) :: xy_PsN(:,:) ! 地表面気圧 (t) real(DBKIND), intent(out) :: xyz_VelLonA(:,:,:)! 速度経度成分 (t+Δt) real(DBKIND), intent(out) :: xyz_VelLatA(:,:,:)! 速度緯度成分 (t+Δt) real(DBKIND), intent(out) :: xyz_VorA(:,:,:) ! 渦度 (t+Δt) real(DBKIND), intent(out) :: xyz_DivA(:,:,:) ! 発散 (t+Δt) real(DBKIND), intent(out) :: xyz_TempA(:,:,:) ! 温度 (t+Δt) real(DBKIND), intent(out) :: xyz_QVapA(:,:,:) ! 比湿 (t+Δt) real(DBKIND), intent(out) :: xy_PsA(:,:) ! 地表面気圧 (t+Δt) !----- 作業用内部変数 ----- ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用) integer(INTKIND) :: i, j, k, kk, l character(STRING), parameter:: subname = "dynamics_leapfrog" continue !---------------------------------------------------------------- ! Check Initialization !---------------------------------------------------------------- call BeginSub(subname) if (.not. dynamics_initialized) then call EndSub( subname, 'Call dynamics_init before call %c', & & c1=trim(subname) ) return endif !------------------------------------------------------------------- ! 地表気圧変化・鉛直σ速度の計算 (連続の式、熱の式で利用される) !------------------------------------------------------------------- xy_lnPs = log( xy_PsN ) w_lnPsN = w_xy(xy_lnPs) xy_GradLonlnPsN = xy_GradLon_w( w_lnPsN ) xy_GradLatlnPsN = xy_GradLat_w( w_lnPsN ) ! 地表面気圧 π の移流 = ｖ・∇π do k = 1, km xyz_lnPsAdv(:,:,k) = & & ( xyz_VelLonN(:,:,k) * xy_GradLonlnPsN & & + xyz_VelLatN(:,:,k) * xy_GradLatlnPsN & & ) / R0 enddo call io_gt4_out_Put('xyz_lnPsAdv', xyz_lnPsAdv) !-- !!$ call DataDump('xyz_lnPsAdv', xyz_lnPsAdv, strlen=80) !++ ! π移流の積み下げ Σ[k〜K] (ｖ・∇π) Δσ xyz_lnPsAdvSum(:,:,km) = xyz_lnPsAdv(:,:, km ) * z_DelSigma( km ) do k = km-1 , 1, -1 xyz_lnPsAdvSum(:,:,k) = & & xyz_lnPsAdvSum(:,:,k+1) & & + xyz_lnPsAdv(:,:,k) * z_DelSigma(k) enddo call io_gt4_out_Put('xyz_lnPsAdvSum', xyz_lnPsAdvSum) !-- !!$ call DataDump('xyz_lnPsAdvSum', xyz_lnPsAdvSum, strlen=80) !++ ! 発散の積下げ Σ[k〜K] ＤΔσ xyz_DivSum(:,:,km) = xyz_DivN(:,:, km ) * z_DelSigma( km ) do k = km -1 , 1, -1 xyz_DivSum(:,:,k) = & & xyz_DivSum(:,:,k+1) + xyz_DivN(:,:,k) * z_DelSigma(k) enddo call io_gt4_out_Put('xyz_DivSum', xyz_DivSum) !-- !!$ call DataDump('xyz_DivSum', xyz_DivSum, strlen=80) !++ !------------------------------------------------------------------- ! Ps の計算。連続の式を解く !------------------------------------------------------------------- !----- π = lnＰs の時間変化 dπ/dt の計算 ----- ! ! π移流の積み下げ Σ[1〜K] (ｖ・∇π) Δσ を格子点上で解く。 ! xy_DlnPsDtN = - xyz_lnPsAdvSum(:,:,1) ! - Σ[1〜K](ｖ・∇π)Δσ !-- !!$ call DataDump('xy_DlnPsDtN', xy_DlnPsDtN, strlen=80) !++ !----- 発散の積み下げ -Σ[1〜K] ＤΔσ をスペクトル空間で加える。----- w_DlnPsDtN = w_xy( xy_DlnPsDtN ) - w_xy( xyz_DivSum(:,:,1) ) !----- A = B + 2Δt * T----- w_lnPsA = w_xy( log( xy_PsB ) ) + 2.0*DelTime * ( w_DlnPsDtN ) !-- !!$ call DataDump('xy_PsB', xy_PsB, strlen=80) !!$ call DataDump('w_lnPsB', w_xy( log( xy_PsB ) ), strlen=80) !!$ call DataDump('w_lnPsA', w_lnPsA, strlen=80) !++ !----- xy_PsA に代入----- xy_PsA = exp( xy_w( w_lnPsA ) ) !-- !!$ call DataDump('xy_PsA', xy_PsA, strlen=80) !++ !------------------------------------------------------------------- ! 鉛直σ速度 (整数、半整数レベル) の上昇流の計算。 !------------------------------------------------------------------- ! σ速度 (3/2 〜 K-1/2) ! - σ[k-1/2] × ( Σ[1〜K](Ｄ + ｖ・∇π)Δσ ) ! - Σ[k〜K] (Ｄ + ｖ・∇π) Δσ do k = 2, km+1 - 1 xyr_VSigma(:,:,k) = & & r_Sigma(k) * ( xyz_lnPsAdvSum(:,:,1) + xyz_DivSum(:,:,1) ) & & - ( xyz_lnPsAdvSum(:,:,k) + xyz_DivSum(:,:,k) ) ! - σ[k-1/2] × ( Σ[1〜K](ｖ・∇π)Δσ ) ! - Σ[k〜K] (ｖ・∇π) Δσ xyr_VSigmaNonG(:,:,k) = & & r_Sigma(k) * xyz_lnPsAdvSum(:,:,1) - xyz_lnPsAdvSum(:,:,k) enddo ! σ速度 (1/2, K+1/2) xyr_VSigma(:,:,1) = 0.0 xyr_VSigma(:,:,km+1) = 0.0 xyr_VSigmaNonG(:,:,1) = 0.0 xyr_VSigmaNonG(:,:,km+1) = 0.0 call io_gt4_out_Put('xyr_VSigma', xyr_VSigma) call io_gt4_out_Put('xyr_VSigmaNonG', xyr_VSigmaNonG) !-- !!$ call DataDump('xyr_VSigma', xyr_VSigma, strlen=80) !++ !------------------------------------------------------------------- ! 温度・仮温度の基本場からのずれ (以降の ＵＡ、ＵＶ、Ｈなどの全てで利用) !------------------------------------------------------------------- do k = 1, km xyz_TempVir(:,:,k) = & & xyz_TempN(:,:,k) * ( 1.0 + (EpsVT * xyz_QVapN(:,:,k)) ) xyz_TempEdd(:,:,k) = xyz_TempN(:,:,k) - z_TempAvrXY(k) xyz_TempVirEdd(:,:,k) = xyz_TempVir(:,:,k) - z_TempAvrXY(k) enddo !-- !!$ call DataDump('xyz_Temp', xyz_Temp, strlen=80) !!$ call DataDump('xyz_TempEdd', xyz_TempEdd, strlen=80) !!$ call DataDump('xyz_TempVir', xyz_TempVir, strlen=80) !!$ call DataDump('xyz_TempVirEdd', xyz_TempVirEdd, strlen=75) !++ !------------------------------------------------------------------- ! 半整数レベルの温度の擾乱 (以降の ＵＡ、ＵＶ、Ｈなどの全てで利用) !------------------------------------------------------------------- do k = 2, km+1 - 1 xyr_TempEdd(:,:,k) = & & z_TempInpolA(k) * xyz_TempN(:,:,k) & & + z_TempInpolB(k-1) * xyz_TempN(:,:,k-1) & & - r_TempAvrXY(k) enddo call io_gt4_out_Put('xyr_TempEdd', xyr_TempEdd) !-- !!$ call DataDump('xyr_TempEdd', xyr_TempEdd, strlen=80) !++ !------------------------------------------------------------------- ! ＵＡ＝Ｖζ＋σ移流、ＶＡ＝Ｕζ＋σ移流 (渦度と発散の式で利用) ! ※ AGCM5 のマニュアルの ＵＡ と異なり、UVFact = cosφは ! 掛かっていないので注意。 !------------------------------------------------------------------- do k = 1, km ! (ζ + ｆ) ｖ ! - (Ｃp κ (1/ａ) Ｔ'v (dπ/dλ) ) / cosφ xyz_UAN(:,:,k) = & & ( xyz_VorN(:,:,k) + xy_Coli ) * xyz_VelLatN(:,:,k) & & - Cp * z_TempInpolKappa(k) & & * xyz_TempVirEdd(:,:,k) * xy_GradLonlnPsN / R0 ! - (ζ + ｆ) ｕ ! - (Ｃp κ (1/ａ) Ｔ'v (dπ/dμ) ) / cosφ xyz_VAN(:,:,k) = & & - ( xyz_VorN(:,:,k) + xy_Coli ) * xyz_VelLonN(:,:,k) & & - Cp * z_TempInpolKappa(k) & & * xyz_TempVirEdd(:,:,k) * xy_GradLatlnPsN / R0 enddo !-- !!$ do k = 1, km !!$ call DataDump('xyz_Div(pi)', & !!$ & xy_w( & !!$ & w_Div_xy_xy( & !!$ & - Cp * z_TempInpolKappa(k) & !!$ & * xyz_TempVirEdd(:,:,k) & !!$ & * xy_GradLon_w( w_xy( xy_lnPs ) ) & !!$ & / R0 & !!$ & , & !!$ & !!$ & - Cp * z_TempInpolKappa(k) & !!$ & * xyz_TempVirEdd(:,:,k) & !!$ & * xy_GradLat_w( w_xy( xy_lnPs ) ) & !!$ & / R0 & !!$ & !!$ & ) & !!$ & ) / R0 & !!$ & !!$ & , strlen=80, multi=(/k/)) !!$ enddo !++ !-- !!$ call DataDump('xyz_VelLon', xyz_VelLon, strlen=80) !!$ call DataDump('xyz_VelLat', xyz_VelLat, strlen=80) !!$ call DataDump('xyz_VAN', xyz_VAN(1,1,:), strlen=80, multi=(/-1/)) !!$ call DataDump('xyz_UAN', xyz_UAN(1,1,:), strlen=80, multi=(/-1/)) !!$ call DataDump('xyz_VAN', xyz_VAN(1,1,:), strlen=80, multi=(/-1/)) !++ do k = 2, km ! - (1/2Δσ[k]) * (dσ/dt)[k-1/2] * (ｕ[k-1] - ｕ[k]) xyz_UAN(:,:,k) = xyz_UAN(:,:,k) & & - 1.0 / (2.0 * z_DelSigma(k)) & & * xyr_VSigma(:,:,k) * ( xyz_VelLonN(:,:,k-1) - xyz_VelLonN(:,:,k) ) ! - (1/2Δσ[k]) * (dσ/dt)[k-1/2] * (ｖ[k-1] - ｖ[k]) xyz_VAN(:,:,k) = xyz_VAN(:,:,k) & & - 1.0 / (2.0 * z_DelSigma(k)) & & * xyr_VSigma(:,:,k) * ( xyz_VelLatN(:,:,k-1) - xyz_VelLatN(:,:,k) ) enddo !-- !!$ call DataDump('xyz_UAN', xyz_UAN(1,1,:), strlen=80, multi=(/-2/)) !!$ call DataDump('xyz_VAN', xyz_VAN(1,1,:), strlen=80, multi=(/-2/)) !++ do k = 1, km - 1 ! - (1/2Δσ[k]) * (dσ/dt)[k+1/2] * (ｕ[k] - ｕ[k+1]) xyz_UAN(:,:,k) = xyz_UAN(:,:,k) & & - 1.0 / (2.0 * z_DelSigma(k)) & & * xyr_VSigma(:,:,k+1) * ( xyz_VelLonN(:,:,k) - xyz_VelLonN(:,:,k+1) ) ! - (1/2Δσ[k]) * (dσ/dt)[k+1/2] * (ｖ[k] - ｖ[k+1]) xyz_VAN(:,:,k) = xyz_VAN(:,:,k) & & - 1.0 / (2.0 * z_DelSigma(k)) & & * xyr_VSigma(:,:,k+1) * ( xyz_VelLatN(:,:,k) - xyz_VelLatN(:,:,k+1) ) enddo !-- !!$ call DataDump('xyz_UAN', xyz_UAN(1,1,:), strlen=80, multi=(/-3/)) !!$ call DataDump('xyz_VAN', xyz_VAN(1,1,:), strlen=80, multi=(/-3/)) !++ !-- !!$ call DataDump('xyz_UAN', xyz_UAN, strlen=80) !!$ call DataDump('xyz_VAN', xyz_VAN, strlen=80) !++ call io_gt4_out_Put('xyz_UAN', xyz_UAN) call io_gt4_out_Put('xyz_VAN', xyz_VAN) !------------------------------------------------------------------- ! 渦度 Vor の計算。渦度方程式を解く !------------------------------------------------------------------- ! dζ/dt の計算 wz_DVorDtN = wa_Div_xya_xya( xyz_VAN, - xyz_UAN ) / R0 !& + w_DiffV(k) * wa_xya(xyz_Vor(:,:,k)) ! 拡散係数の導出が厄介なので後回し !----- A = B + 2Δt * T----- wz_VorA = wa_xya( xyz_VorB ) + 2.0d0*DelTime * ( wz_DVorDtN ) !!! & + 2.0d0*DelTime * wa_NumVis_wa( wa_xya( xyz_VorB ) ) !----- xyz_VorA に代入----- xyz_VorA = xya_wa( wz_VorA ) !-- !!$ call DataDump('wz_DVorDtN', wz_DVorDtN, strlen=70) !!$ call DataDump('xyz_VorB', xyz_VorB, strlen=80) !!$ call DataDump('xyz_Vor', xyz_Vor, strlen=80) !!$ call DataDump('xyz_Vor_T', xya_wa(wz_DVorDtN), strlen=80) !!$ call DataDump('xyz_VorA', xyz_VorA, strlen=80) !++ !------------------------------------------------------------------- ! Ｅ=ｕ**2+ｖ**2 + 仮温度補正 ( ΣＷ(Ｔv-Ｔ) ) (発散の式で利用) !------------------------------------------------------------------- ! 仮温度補正 ( ΣＷ(Ｔv-Ｔ) ) の計算。格子点上で積み上げ。 ! 初期化 xyz_KE = 0.0 xyz_KE(:,:,1) = Cp * z_HydroAlpha(1) & & * ( xyz_TempVir(:,:,1) - xyz_TempN(:,:,1) ) do k = 2, km xyz_KE(:,:,k) = xyz_KE(:,:,k-1) & & + Cp * z_HydroAlpha(k) & & * ( xyz_TempVir(:,:,k) - xyz_TempN(:,:,k) ) & & + Cp * z_HydroBeta(k-1) & & * ( xyz_TempVir(:,:,k-1) - xyz_TempN(:,:,k-1) ) enddo !-- !!$ call DataDump('xyz_Temp', xyz_Temp, strlen=80) !!$ call DataDump('xyz_TempVir', xyz_TempVir, strlen=80) !!$ call DataDump('xyz_KE-w(Tv-T)', xyz_KE, strlen=80) !++ ! 運動エネルギー ｕ**2+ｖ**2 の加算 xyz_KE = xyz_KE + ( xyz_VelLonN**2 + xyz_VelLatN**2 ) / 2.0 !-- !!$ call DataDump('xyz_KE-u**2+v**2', xyz_KE, strlen=80) !++ call io_gt4_out_Put('xyz_KE', xyz_KE) !------------------------------------------------------------------- ! 発散 Div の計算。発散方程式を解く !------------------------------------------------------------------- ! ＷＴ = Ｃp α Ｔ[k] + Ｃp β Ｔ[k-1] をスペクトル空間で積み上げ ! 初期化 wz_PresTendTemp = 0.0 wz_TempN = wa_xya(xyz_TempN) wz_PresTendTemp(:,1) = Cp * z_HydroAlpha(1) * wz_TempN(:,1) do k = 2, km wz_PresTendTemp(:,k) = wz_PresTendTemp(:,k-1) & & + Cp * z_HydroAlpha(k) * wz_TempN(:,k) & & + Cp * z_HydroBeta(k-1) * wz_TempN(:,k-1) enddo !-- !!$ call io_gt4_out_Put('xyz_PresTendTemp', xya_wa(wa_Lapla_wa(wz_PresTendTemp))/R0**2) !++ !-- !!$ call DataDump('wz_PresTendTemp', wz_PresTendTemp, strlen=70) !!$ call DataDump('xyz_Lapla(wz_PresTendTemp)/R0**2', & !!$ & xya_wa( wa_Lapla_wa( wz_PresTendTemp ) )/R0**2, strlen=60) !++ ! Ｇπ = (Ｃp κ Ｔ) π do k = 1, km wz_PresTendPs(:,k) = & & Cp * z_TempInpolKappa(k) * z_TempAvrXY(k) * w_xy( log( xy_PsN ) ) enddo !-- !!$ call io_gt4_out_Put('xyz_PresTendPs', xya_wa(wa_Lapla_wa(wz_PresTendPs))/R0**2) !++ !-- !!$ call DataDump('wz_PresTendPs', wz_PresTendPs, strlen=70) !!$ call DataDump('xyz_Lapla(wz_PresTendPs)/R0**2', & !!$ & xya_wa( wa_Lapla_wa( wz_PresTendPs ) )/R0**2, strlen=60) !++ !-- !!$ call DataDump('Div(wz_UAN - wz_VAN)', & !!$ & wa_Div_xya_xya( xyz_UAN, - xyz_VAN ), & !!$ & strlen=70) !!$ call DataDump('Div(wz_UAN - wz_VAN)/R0', & !!$ & wa_Div_xya_xya( xyz_UAN, - xyz_VAN ) /R0, & !!$ & strlen=70) !!$ call DataDump('xyz_Div(wz_UAN - wz_VAN)/R0', & !!$ & xya_wa( wa_Div_xya_xya( xyz_UAN, - xyz_VAN ) ) /R0, & !!$ & strlen=70) !!$ call DataDump('Lapla(wz_KE)', & !!$ & wa_Lapla_wa( wa_xya(xyz_KE) ), strlen=70) !++ ! dＤ/dt の計算 ! - ∇**2 (Φ + ＷＴ + Ｇπ) wz_DDivDtN = & & wa_Div_xya_xya( xyz_UAN, xyz_VAN ) / R0 & & - wa_Lapla_wa( wa_xya(xyz_KE) ) / R0**2 & ! !& + w_DiffV(k) * wa_xya(xyz_Div(:,:,k)) !現在, 拡散の効果は !dynamics_diffusion で与えている. ! semi-implicitだとここには無い項達 & - wa_Lapla_wa( & !& w_xy( xy_Phi ) & ! 地表ジオポテンシャルΦは後回し & wz_PresTendTemp & ! ＷＴ & + wz_PresTendPs & ! Ｇπ = (Ｃp κ Ｔ)π & ) / R0**2 !----- A = B + 2Δt * T----- wz_DivA = wa_xya( xyz_DivB ) + 2.0d0*DelTime * ( wz_DDivDtN ) !!! & + 2.0d0*DelTime * wa_NumVis_wa( wa_xya( xyz_DivB ) ) !----- xyz_DivA に代入----- xyz_DivA = xya_wa( wz_DivA ) !-- !!$ call DataDump('xyz_DivB', xyz_DivB, strlen=80) !!$ call DataDump('xyz_Div', xyz_Div, strlen=80) !!$ call DataDump('xyz_Div_T', xya_wa(wz_DDivDtN), strlen=80) !!$ call DataDump('xyz_DivA', xyz_DivA, strlen=80) !++ !------------------------------------------------------------------- ! Ｈ＝水平発散 Ｔ'Ｄ + σ移流 + π変化の効果 (熱力学の式) !------------------------------------------------------------------- ! 温度の局所時間変化項 Ｈ の計算 ! Ｔ' Ｄ ! κＴv (ｖ・∇π) : πの移流の効果 ! - (α/Δσ) {Ｔv Σ[k〜K] (ｖ・∇π) + Ｔv' Σ[k〜K] (ＤΔσ)} do k = 1, km xyz_DTempLocalDtN(:,:,k) = & & xyz_TempEdd(:,:,k) * xyz_DivN(:,:,k) & & + z_TempInpolKappa(k) * xyz_TempVir(:,:,k) & & * xyz_lnPsAdv(:,:,k) & & - z_HydroAlpha(k) / z_DelSigma(k) & & * ( & & xyz_TempVir(:,:,k) * xyz_lnPsAdvSum(:,:,k) & & + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k) & & ) enddo ! - (dσ/dt) * (１/Δσ) * (Ｔ'[k-1/2] − Ｔ'[k]) ! - (dσ/dt)^NG * (１/Δσ) * (Ｔ￣[k-1/2] − Ｔ￣[k]) do k = 2, km xyz_DTempLocalDtN(:,:,k) = xyz_DTempLocalDtN(:,:,k) & & - xyr_VSigma(:,:,k) / z_DelSigma(k) & & * ( xyr_TempEdd(:,:,k) - xyz_TempEdd(:,:,k) ) & & - xyr_VSigmaNonG(:,:,k) / z_DelSigma(k) & & * ( r_TempAvrXY(k) - z_TempAvrXY(k) ) enddo ! - (dσ/dt) * (１/Δσ) * (Ｔ'[k] − Ｔ'[k+1/2]) ! - (dσ/dt)^NG * (１/Δσ) * (Ｔ￣[k] − Ｔ￣[k+1/2]) ! - (β/Δσ) { Ｔv Σ[k+1〜K] (ｖ・∇π) ! + Ｔv' Σ[k+1〜K] (ＤΔσ)} do k = 1, km - 1 xyz_DTempLocalDtN(:,:,k) = xyz_DTempLocalDtN(:,:,k) & & - xyr_VSigma(:,:,k+1) / z_DelSigma(k) & & * ( xyz_TempEdd(:,:,k) - xyr_TempEdd(:,:,k+1) ) & & - xyr_VSigmaNonG(:,:,k+1) / z_DelSigma(k) & & * ( z_TempAvrXY(k) - r_TempAvrXY(k+1) ) & & - z_HydroBeta(k) / z_DelSigma(k) & & * ( & & xyz_TempVir(:,:,k) * xyz_lnPsAdvSum(:,:,k+1) & & + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k+1) & & ) enddo !-- !!$ call io_gt4_out_Put('xyz_DTempLocalDtN_lnPsAdvSum', & !!$ & - z_HydroAlpha(2) / z_DelSigma(2) & !!$ & * ( & !!$ & xyz_TempVir(:,:,2) * xyz_lnPsAdvSum(:,:,2) & !!$ & ) & !!$ & - z_HydroBeta(2) / z_DelSigma(2) & !!$ & * ( & !!$ & xyz_TempVir(:,:,2) * xyz_lnPsAdvSum(:,:,2+1) & !!$ & ) & !!$ & ) !++ !-- !!$ call DataDump('xyz_DivSum*TempVirEdd*(Alpha+Beta)', & !!$ & - z_HydroAlpha(k) / z_DelSigma(k) & !!$ & * ( xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k) ) & !!$ & - z_HydroBeta(k) / z_DelSigma(k) & !!$ & * ( xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k+1) ) & !!$ & , strlen=70) !++ !-- !!$ call DataDump('xyz_DTempLocalDtN', xyz_DTempLocalDtN, strlen=80) !++ call io_gt4_out_Put('xyz_DTempLocalDtN', xyz_DTempLocalDtN) !-- !------------------------------------------------------------------- ! ＵＴ' 、ＶＴ' (熱力学の式) !------------------------------------------------------------------- !!$ DO 5100 K = 1, KMAX !!$ DO 5100 IJ = 1, IDIM*JDIM !!$ GTUT( IJ,K ) = GAU ( IJ,K ) * GATED ( IJ,K ) !!$ & * UVFACT ( IJ ) !!$ GTVT( IJ,K ) = GAV ( IJ,K ) * GATED ( IJ,K ) !!$ & * UVFACT ( IJ ) !!$ 5100 CONTINUE !++ !------------------------------------------------------------------- ! 温度 Ｔ の計算。熱力学の式を解く !------------------------------------------------------------------- ! 非重力波項成分 (非線形項) の部分に関する温度変化を計算 ! ! dＵＴ'/dλ + dＶＴ'/dμ ! 温度の局所時間変化項 Ｈ do k = 1, km wz_DTempDtN(:,k) = & & - w_Div_xy_xy( & & xyz_VelLonN(:,:,k) * xyz_TempEdd(:,:,k), & & xyz_VelLatN(:,:,k) * xyz_TempEdd(:,:,k) & & ) / R0 & & + w_xy( xyz_DTempLocalDtN(:,:,k) ) !& ! 非断熱加熱 Ｑ もまだ入ってないよ !& ! 摩擦熱および熱拡散は後回し enddo !-- !!$ call DataDump('xyz_Temp_T_NonGrav', xya_wa(wz_DTempDtN), strlen=70) !++ ! 重力波項成分 (線形項) の部分に関する温度変化を計算 ! ! ｈＤ (重力波項成分によるπ変化の効果) do k = 1, km do kk = 1, km wz_DTempDtN(:,k) = wz_DTempDtN(:,k) & & - zz_TempMtx(k,kk) * w_xy( xyz_DivN(:,:,kk) ) !& ! 拡散項は後回し enddo enddo !----- A = B + 2Δt * T----- wz_TempA = wa_xya( xyz_TempB ) + 2.0d0*DelTime * ( wz_DTempDtN ) !!! & + 2.0d0*DelTime * wa_NumVisScaler_wa( wa_xya( xyz_TempB ) ) !----- xyz_TempA に代入----- xyz_TempA = xya_wa( wz_TempA ) call io_gt4_out_Put( 'xyz_Temp_T', xya_wa( wz_DTempDtN ) ) !-- !!$ call DataDump('xyz_TempB', xyz_TempB, strlen=80) !!$ call DataDump('xyz_Temp', xyz_Temp, strlen=80) !!$ call DataDump('xyz_Temp_T', xya_wa(wz_DTempDtN), strlen=80) !!$ call DataDump('xyz_TempA', xyz_TempA, strlen=80) !++ !------------------------------------------------------------------- ! 比湿 QVap の計算。水蒸気の式を解く。 !------------------------------------------------------------------- ! Ｒ＝ｑＤ＋σ移流 (非重力波項) を格子点上で計算 xyz_QVapDivVSigmaAdv = xyz_QVapN * xyz_DivN do k = 2, km xyz_QVapDivVSigmaAdv(:,:,k) = xyz_QVapDivVSigmaAdv(:,:,k) & & - xyr_VSigma(:,:,k) / ( 2.0 * z_DelSigma(k) ) & & * ( xyz_QVapN(:,:,k-1) - xyz_QVapN(:,:,k) ) enddo do k = 1, km - 1 xyz_QVapDivVSigmaAdv(:,:,k) = xyz_QVapDivVSigmaAdv(:,:,k) & & - xyr_VSigma(:,:,k+1) / ( 2.0 * z_DelSigma(k) ) & & * ( xyz_QVapN(:,:,k) - xyz_QVapN(:,:,k+1) ) enddo ! 水平移流 dＵｑ/dλ + dＶｑ/dμ と水平拡散項、水蒸気ソースを計算 ! dＵｑ/dλ + dＶｑ/dμ ! Ｒ＝ｑＤ＋σ移流 do k = 1, km wz_DQVapDtN(:,k) = & & - w_Div_xy_xy( & & xyz_VelLonN(:,:,k) * xyz_QVapN(:,:,k), & & xyz_VelLatN(:,:,k) * xyz_QVapN(:,:,k) & & ) / R0 & & + w_xy( xyz_QVapDivVSigmaAdv(:,:,k) ) !& ! 水平拡散項、水蒸気ソースは後回し enddo !----- A = B + 2Δt * T----- wz_QVapA = wa_xya( xyz_QVapB ) + 2.0d0*DelTime * ( wz_DQVapDtN ) !!! & + 2.0d0*DelTime * wa_NumVisScaler_wa( wa_xya( xyz_QVapB ) ) !----- xyz_QVapA に代入----- xyz_QVapA = xya_wa( wz_QVapA ) !-- !!$ call DataDump('xyz_QVapB', xyz_QVapB, strlen=80) !!$ call DataDump('xyz_QVap', xyz_QVap, strlen=80) !!$ call DataDump('xyz_QVap_T', xya_wa(wz_DQVapDtN), strlen=80) !!$ call DataDump('xyz_QVapA', xyz_QVapA, strlen=80) !++ call EndSub(subname) end subroutine dynamics_leapfrog subroutine dynamics_diffusion( & & xyz_VorB , xyz_DivB , xyz_TempB , xyz_QVapB , & & xyz_VorA , xyz_DivA , xyz_TempA , xyz_QVapA ) ! Calculate Diffusion Term ! ! t-Δt の値から水平拡散項を求め、それを t+Δt の値に加える。 ! use type_mod, only: STRING, REKIND, DBKIND, INTKIND use time_mod, only: DelTime, CurrentTime use grid_3d_mod, only: km use grid_wavenumber_mod, only: nm use spml_mod, only: wa_xya, xya_wa, l_nm use io_gt4_out_mod,only: io_gt4_out_Put use dc_trace, only: DbgMessage, BeginSub, EndSub, DataDump implicit none real(DBKIND), intent(in) :: xyz_VorB(:,:,:) ! 渦度 (t-Δt) real(DBKIND), intent(in) :: xyz_DivB(:,:,:) ! 発散 (t-Δt) real(DBKIND), intent(in) :: xyz_TempB(:,:,:) ! 温度 (t-Δt) real(DBKIND), intent(in) :: xyz_QVapB(:,:,:) ! 比湿 (t-Δt) real(DBKIND), intent(inout) :: xyz_VorA(:,:,:) ! 渦度 (t+Δt) real(DBKIND), intent(inout) :: xyz_DivA(:,:,:) ! 発散 (t+Δt) real(DBKIND), intent(inout) :: xyz_TempA(:,:,:) ! 温度 (t+Δt) real(DBKIND), intent(inout) :: xyz_QVapA(:,:,:) ! 比湿 (t+Δt) real(DBKIND) :: wz_VorA((nm+1)**2, km) ! 渦度 (t+Δt) real(DBKIND) :: wz_DivA((nm+1)**2, km) ! 発散 (t+Δt) real(DBKIND) :: wz_TempA((nm+1)**2, km) ! 温度 (t+Δt) real(DBKIND) :: wz_QVapA((nm+1)**2, km) ! 比湿 (t+Δt) integer(INTKIND) :: i character(STRING), parameter:: subname = "dynamics_diffusion" continue !---------------------------------------------------------------- ! Check Initialization !---------------------------------------------------------------- call BeginSub(subname) if (.not. dynamics_initialized) then call EndSub( subname, 'Call dynamics_init before call %c', & & c1=trim(subname) ) return endif !------------------------------------------------------------------- ! それぞれの量の拡散項を計算 !------------------------------------------------------------------- xyz_VorA = & & xya_wa( wa_xya(xyz_VorA) & & + 2.0d0*DelTime * wz_DiffVorDiv * wa_xya( xyz_VorB ) ) !-- !!$ call io_gt4_out_Put('xyz_Vor_Diffusion', & !!$ & xya_wa( 2.0d0*DelTime & !!$ & * wz_DiffVorDiv * wa_xya( xyz_VorB ) ) ) !++ xyz_DivA = & & xya_wa( wa_xya(xyz_DivA) & & + 2.0d0*DelTime * wz_DiffVorDiv * wa_xya( xyz_DivB ) ) !-- !!$ call io_gt4_out_Put('xyz_Div_Diffusion', & !!$ & xya_wa( 2.0d0*DelTime & !!$ & * wz_DiffVorDiv * wa_xya( xyz_DivB ) ) ) !++ xyz_TempA = & & xya_wa( wa_xya(xyz_TempA) & & + 2.0d0*DelTime * wz_DiffTherm * wa_xya( xyz_TempB ) ) !-- !!$ call io_gt4_out_Put('xyz_Temp_Diffusion', & !!$ & xya_wa( 2.0d0*DelTime & !!$ & * wz_DiffTherm * wa_xya( xyz_TempB ) ) ) !++ xyz_QVapA = & & xya_wa( wa_xya(xyz_QVapA) & & + 2.0d0*DelTime * wz_DiffTherm * wa_xya( xyz_QVapB ) ) !-- !!$ call io_gt4_out_Put('xyz_QVap_Diffusion', & !!$ & xya_wa( 2.0d0*DelTime & !!$ & * wz_DiffTherm * wa_xya( xyz_QVapB ) ) ) !++ ! スペクトルデータとして見て拡散が効いているか確認するための ! データ出力を記述したソース。本計算にはまったく必要でない。 !!wz_VorA = wa_xya(xyz_VorA ) !!wz_DivA = wa_xya(xyz_DivA ) !!wz_TempA = wa_xya(xyz_TempA) !!wz_QVapA = wa_xya(xyz_QVapA) !! !!write(80, *) CurrentTime, & !! & wz_VorA( lNm(21,0), 1) , wz_DivA( lNm(21,0), 1) , & !! & wz_TempA( lNm(21,0), 1) , wz_QVapA( lNm(21,0), 1) !!call flush(80) call EndSub(subname) end subroutine dynamics_diffusion subroutine dynamics_diagnostic & & ( x_Lon , y_Lat , z_Sigma , r_Sigma , & !(in) & xyz_VelLonA, xyz_VelLatA, & !(out) & xyz_VorA , xyz_DivA , xyz_TempA, xyz_QVapA, xy_PsA ) !(in) ! ! Calculate Diagnostic Values. ! ! 診断的に得られる量の演算を行なう。 ! 現在は渦度発散から速度成分 (経度方向、緯度方向) の演算を行なうのみである。 ! use type_mod, only: STRING, REKIND, DBKIND, INTKIND use grid_3d_mod, only: im, jm, km use grid_wavenumber_mod, only: nm use constants_mod, only: R0, Grav use spml_mod, only: wa_xya, xya_GradLon_wa, xya_GradLat_wa, & & wa_LaplaInv_wa, IntLonLat_xy use io_gt4_out_mod,only: io_gt4_out_Put use dc_trace, only: DbgMessage, BeginSub, EndSub, DataDump use dc_string, only: toChar implicit none real(DBKIND), intent(in):: x_Lon(:) ! 経度座標 real(DBKIND), intent(in):: y_Lat(:) ! 緯度座標 real(DBKIND), intent(in):: z_Sigma(:) ! σレベル(整数)座標 real(DBKIND), intent(in):: r_Sigma(:) ! σレベル(半整数)座標 real(DBKIND), intent(in):: xyz_VorA(:,:,:) ! 渦度 (t+Δt) real(DBKIND), intent(in):: xyz_DivA(:,:,:) ! 発散 (t+Δt) real(DBKIND), intent(in):: xyz_TempA(:,:,:)! 温度 (t+Δt) real(DBKIND), intent(in):: xyz_QVapA(:,:,:)! 比湿 (t+Δt) real(DBKIND), intent(in):: xy_PsA(:,:) ! 地表面気圧 (t+Δt) real(DBKIND), intent(out) :: xyz_VelLonA(:,:,:) ! 速度経度成分 (t+Δt) real(DBKIND), intent(out) :: xyz_VelLatA(:,:,:) ! 速度緯度成分 (t+Δt) real(DBKIND) :: TotalMass ! 全質量 !----- 作業用内部変数 ----- ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用) integer(INTKIND) :: i, j, k, l character(STRING), parameter:: subname = "dynamics_diagnostic" !---------------------------------------------------------------- ! Check Initialization !---------------------------------------------------------------- call BeginSub(subname) if (.not. dynamics_initialized) then call EndSub( subname, 'Call dynamics_init before call %c', & & c1=trim(subname) ) return endif !------------------------------------------------------------------- ! ｕ と ｖ を出力 !------------------------------------------------------------------- wz_PsiA = wa_LaplaInv_wa( wa_xya( xyz_VorA ) ) * R0**2 wz_ChiA = wa_LaplaInv_wa( wa_xya( xyz_DivA ) ) * R0**2 xyz_VelLonA = & & ( xya_GradLon_wa( wz_ChiA ) - xya_GradLat_wa( wz_PsiA ) ) / R0 xyz_VelLatA = & & ( xya_GradLon_wa( wz_PsiA ) + xya_GradLat_wa( wz_ChiA ) ) / R0 !-- !!$ call DataDump('xyz_VelLonA', xyz_VelLonA, strlen=80) !!$ call DataDump('xyz_VelLatA', xyz_VelLatA, strlen=80) !++ !------------------------------------------------------------------- ! xy_Ps から全質量を求める。 !------------------------------------------------------------------- TotalMass = IntLonLat_xy( xy_PsA / Grav ) call io_gt4_out_Put('TotalMass', TotalMass) call EndSub(subname) end subroutine dynamics_diagnostic subroutine dynamics_end ! ! Terminate module ! ! dynamics_init で allocate した変数を deallocate し、 ! 演算した値も全て破棄する。 ! use type_mod, only: STRING, REKIND, DBKIND, INTKIND use dc_trace, only: BeginSub, EndSub, DbgMessage implicit none !----------------------------------------------------------------- ! 変数定義 !----------------------------------------------------------------- !----- 作業用内部変数 ----- character(STRING), parameter:: subname = "dynamics_end" continue !----------------------------------------------------------------- ! Check Initialization !----------------------------------------------------------------- call BeginSub(subname) if ( .not. dynamics_initialized) then call EndSub( subname, 'dynamics_init was not called', & & c1=trim(subname) ) return else dynamics_initialized = .false. endif !----------------------------------------------------------------- ! allocate した変数の nullify とか !----------------------------------------------------------------- deallocate( xy_Coli , & & xy_UVFact , & & z_DelSigma , & & wz_DiffVorDiv , & & wz_DiffTherm , & & z_HydroAlpha , & & z_HydroBeta , & & z_TempInpolA , & & z_TempInpolB , & & z_TempInpolKappa , & & z_TempAvrXY , & & r_TempAvrXY , & & zz_DivMtx , & & zz_TempMtx , & & zz_TempMtxQ , & & zz_TempMtxS , & & zz_TempMtxR ) deallocate & & ( xy_lnPs , & ! π＝ ln Ｐs (t) & xyz_lnPsAdv , & ! π の移流 ｖ・∇π & xyz_lnPsAdvSum, & ! π移流積下げ Σ[k〜K](ｖ・∇π)Δσ & xyz_DivSum , & ! 発散の積下げ Σ[k〜K] ＤΔσ & xy_DlnPsDtN , & ! πの時間変化 dπ/dt (Δt) & w_DlnPsDtN , & ! πの時間変化 (スペクトル) (Δt) & w_lnPsA , & ! π (t+Δt) & xyr_VSigma , & ! 鉛直σ速度(半整数レベル) (t) & xyr_VSigmaNonG , & ! 鉛直σ速度 非重力波成分 (半整数レベル) (t) & xyz_TempEdd , &! Ｔ' ：温度の擾乱 (整数レベル) (t) & xyz_TempVir , &! Ｔv ：仮温度 (virtual temperature)(t) & xyz_TempVirEdd , &! Ｔv'：仮温度の擾乱 (t) & xyr_TempEdd, &! Ｔ' ：温度の擾乱 (半整数レベル) (t) & xyz_UAN , & ! 東西運動量変化項ＵＡ & xyz_VAN , & ! 南北運動量変化項ＶＡ & wz_DVorDtN , & ! 渦度変化のスペクトルデータ (Δt) & wz_VorA , & ! 渦度のスペクトルデータ (t+Δt) & xyz_KE , & ! 運動エネルギー + 仮温度補正 ! ｕ**2+ｖ**2 + ΣＷ(Ｔv-Ｔ) & wz_DDivDtN , & ! 発散のスペクトルデータ (Δt) & wz_DivA , & ! 発散のスペクトルデータ (t+Δt) & wz_PresTendTemp , & ! 温度による圧力傾度のスペクトルデータ & wz_PresTendPs , & ! 地表圧力による圧力傾度のスペクトルデータ & xyz_DTempLocalDtN , & ! 温度の局所時間変化項 Ｈ ! (水平発散 Ｔ'Ｄ + σ移流 ! + π変化の効果) (全て非重力波項) & wz_DTempDtN , & ! 温度の変化スペクトルデータ (Δt) & wz_TempA , & ! 温度のスペクトルデータ (t+Δt) & xyz_QVapDivVSigmaAdv, & ! Ｒ＝ｑＤ＋σ移流 & wz_DQVapDtN , & ! 比湿の変化スペクトルデータ (Δt) & wz_QVapA & ! 比湿のスペクトルデータ (t+Δt) & ) deallocate & & ( wz_PsiA , & ! スペクトル(流線関数) & wz_ChiA & ! スペクトル(ポテンシャル) & ) call EndSub(subname) end subroutine dynamics_end !-- !! !=== Function wa_NumVis_wa : 数値粘性項の計算 !! ! !! !数値粘性項を計算する内部関数である。 !! ! !! function wa_NumVis_wa(wa_data) !! !==== Dependency !! use type_mod , only: STRING, REKIND, DBKIND, INTKIND !! use constants_mod , only: R0, VisOrder, EFoldTime !! use grid_3d_mod , only: km !! use grid_wavenumber_mod, only: nm !! use dc_trace , only: BeginSub, EndSub, DbgMessage, DataDump !! implicit none !! !==== Input !! ! !! real(DBKIND), intent(in) :: wa_data((nm+1)*(nm+1), km) !! ! !! !==== Output !! ! !! real(DBKIND) :: wa_NumVis_wa((nm+1)*(nm+1), km) !! !! real(DBKIND) :: VisCoef ! 超粘性係数 !! character(STRING), parameter:: subname = "wa_NumVis_wa" !! continue !! call BeginSub(subname) !! !! VisCoef = ( real( ( nm*(nm+1) ), DBKIND ) & !! & / R0**2 )**(-VisOrder/2) & !! & / EFoldTime !! !!!!$ call DbgMessage('VisCoef=<%f>', d=(/VisCoef/)) !! !! wa_NumVis_wa = & !! & wa_NumVisScaler_wa(wa_data) & !! & - VisCoef & !! & * ( - (2.0d0/R0**2)**(VisOrder/2) ) * wa_data !! !!!!$ call DbgMessage('-VisCoef*(2/R0**2)**(VisOrder/2)=<%f>', & !!!!$ & d=(/ - VisCoef * ( - (2.0d0/R0**2)**(VisOrder/2) ) /) ) !! !! !!!!!$ hdc = ( dble( (hdmaxn*(hdmaxn+1) ) ) / pradisq )**(-hdord/2) & !!!!!$ / hdts !!!!!$ !!!!!$ hdord : 水平拡散の次数（ナブラの hdord 乗、ラプラシアンの hdord/2 乗） !!!!!$ hdmaxn : 打ち切り波数（ T21 なら 21、T42 なら 42、... ） !!!!!$ pradisq : 惑星半径の二乗 [m2] !!!!!$ hdts : 全波数 hdmaxn における減衰の時定数 [s] !! !! call EndSub(subname) !! end function wa_NumVis_wa !! !! !=== Function wa_NumVisScaler_wa : スカラーに対する数値粘性項の計算 !! ! !! !スカラーに対する数値粘性項を計算する内部関数である。 !! ! !! function wa_NumVisScaler_wa(wa_data) !! !==== Dependency !! use type_mod , only: STRING, REKIND, DBKIND, INTKIND !! use constants_mod , only: R0, VisOrder, EFoldTime !! use grid_3d_mod , only: km !! use grid_wavenumber_mod, only: nm !! use spml_mod , only: rn !! use dc_trace , only: BeginSub, EndSub, DbgMessage, DataDump !! implicit none !! !==== Input !! ! !! real(DBKIND), intent(in) :: wa_data((nm+1)*(nm+1), km) !! ! !! !==== Output !! ! !! real(DBKIND) :: wa_NumVisScaler_wa((nm+1)*(nm+1), km) !! !! real(DBKIND) :: VisCoef ! 超粘性係数 !! real(DBKIND) :: wa_rn((nm+1)*(nm+1), km) !! integer(INTKIND) :: k !! character(STRING), parameter:: subname = "wa_NumVisScaler_wa" !! continue !! call BeginSub(subname) !! !! do k = 1, km !! wa_rn(:,k) = rn(:,1) !! enddo !! !! VisCoef = ( real( ( nm*(nm+1) ), DBKIND ) & !! & / R0**2 )**(-VisOrder/2) & !! & / EFoldTime !! !!!!$ call DbgMessage('VisCoef=<%f>', d=(/VisCoef/)) !! !! wa_NumVisScaler_wa = & !! & - VisCoef * ( ( - wa_rn / R0**2 )**(VisOrder/2) ) * wa_data !! !! call EndSub(subname) !! end function wa_NumVisScaler_wa !++ end module dynamics_mod