| Class | dynamics_hspl_vas83 |
| In: |
dynamics/dynamics_hspl_vas83.F90
|
Note that Japanese and English are described in parallel.
力学過程を演算するモジュールです. 水平離散化にスペクトル法 (Bourke, 1988) を, 鉛直離散化には Arakawa and Suarez (1983) を用いています.
時間積分法にはリープフロッグスキームを用いています. デフォルトでは, $ Delta t $ を大きくとるために, 重力波項に セミインプリシット法を適用しています. NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで, 重力波項をエクスプリシット法によって解くことも可能です.
This is a dynamical core module. Spectral method (Bourke, 1988) (for horizontal) and Arakawa and Suarez (1983) method (for vertical) are used.
Leap-frog scheme is used as time integration. By default, semi-implicit scheme is applied to gravitational terms in order to enlarge the value of $ Delta t $ . Explicit scheme can be applied to gravitational terms by changing "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
| DynamicsHSplVAS83 : | 力学計算 |
| DynamicsHSplVAS83Init : | 初期化 |
| DynamicsHSplVAS83Finalize : | 終了処理 (モジュール内部の変数の割り付け解除) |
| ——————— : | ———— |
| DynamicsHSplVAS83 : | Calculate dynamics |
| DynamicsHSplVAS83Init : | Initialization |
| DynamicsHSplVAS83Finalize : | Termination (deallocate variables in this module) |
NAMELIST#dynamics_hspl_vas83_nml
| Subroutine : | |||
| xyz_UB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_VB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_TempB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_PsB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_UN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_VN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_TempN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_PsN(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_DUDtPhy(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DVDtPhy(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DTempDtPhy(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyzf_DQMixDtPhy(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_SurfHeight(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_UA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_VA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_TempA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xy_PsA(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xyz_ArgOMG(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out), optional
|
力学過程の演算を行い, 与えられた $ t-\Delta t $ および $ t $ の 東西風速 (xyz_UB, xyz_UN), 南北風速 (xyz_VB, xyz_VN), 温度 (xyz_TempB, xyz_TempN), 質量混合比 (xyzf_QMixB, xyzf_QMixN), 地表面気圧 (xyz_PsB, xyz_PsN), および地表面高度 (xy_SurfHeight) から, $ t+Delta t $ の 東西風速 (xyz_UA), 南北風速 (xyz_VA), 温度 (xyz_TempA), 質量混合比 (xyzf_QMixA), 地表面気圧 (xyz_PsA) を返します.
別の物理プロセスによる渦度, 発散, 温度, 比湿の変化を, 力学過程の変化に足して次のステップを計算する場合には, それらの変化を xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy に与えてください.
時間積分法にはリープフロッグスキームを用いています. デフォルトでは, $ Delta t $ を大きくとるために, 重力波項に セミインプリシット法を適用しています. NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで, 重力波項をエクスプリシット法によって解くことも可能です.
Calculate dynamical processes. Input data are Eastward wind (xyz_UB, xyz_UN), northward wind (xyz_VB, xyz_VN), temperature (xyz_TempB, xyz_TempN), mass mixing ratio (xyzf_QMixB, xyzf_QMixN), surface pressure (xyz_PsB, xyz_PsN) at $ t-\Delta t $ and $ t $, and surface height (xy_SurfHeight). Eastward wind (xyz_UA), northward wind (xyz_VA), temperature (xyz_TempA), mass mixing ratio (xyzf_QMixA), surface pressure (xyz_PsA) are returned.
In order to add tendencies of vorticity, divergence, temperature and specific humidity calculated by physical processes other than those by this dynamical process, these tendencies should be given to "xyz_DUDtPhy", "xyz_DVDtPhy", "xyz_DTempDtPhy", "xyzf_DQMixDtPhy"
Leap-frog scheme is used as time integration. By default, semi-implicit scheme is applied to gravitational terms for extension of $ Delta t $ . Explicit scheme can be applied to gravitational terms by changing "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
subroutine DynamicsHSplVAS83( xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy, xy_SurfHeight, xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA, xyz_ArgOMG )
!
! 力学過程の演算を行い, 与えられた $ t-\Delta t $ および $ t $ の
! 東西風速 (xyz_UB, xyz_UN), 南北風速 (xyz_VB, xyz_VN),
! 温度 (xyz_TempB, xyz_TempN), 質量混合比 (xyzf_QMixB, xyzf_QMixN),
! 地表面気圧 (xyz_PsB, xyz_PsN), および地表面高度 (xy_SurfHeight) から,
! $ t+\Delta t $ の
! 東西風速 (xyz_UA), 南北風速 (xyz_VA), 温度 (xyz_TempA),
! 質量混合比 (xyzf_QMixA), 地表面気圧 (xyz_PsA) を返します.
!
! 別の物理プロセスによる渦度, 発散, 温度, 比湿の変化を,
! 力学過程の変化に足して次のステップを計算する場合には,
! それらの変化を xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy
! に与えてください.
!
! 時間積分法にはリープフロッグスキームを用いています.
! デフォルトでは, $ \Delta t $ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! Calculate dynamical processes.
! Input data are
! Eastward wind (xyz_UB, xyz_UN),
! northward wind (xyz_VB, xyz_VN),
! temperature (xyz_TempB, xyz_TempN),
! mass mixing ratio (xyzf_QMixB, xyzf_QMixN),
! surface pressure (xyz_PsB, xyz_PsN) at $ t-\Delta t $ and $ t $,
! and surface height (xy_SurfHeight).
! Eastward wind (xyz_UA), northward wind (xyz_VA),
! temperature (xyz_TempA),
! mass mixing ratio (xyzf_QMixA), surface pressure (xyz_PsA)
! are returned.
!
! In order to add tendencies of vorticity, divergence,
! temperature and specific humidity calculated by
! physical processes other than those by
! this dynamical process,
! these tendencies should be given to
! "xyz_DUDtPhy", "xyz_DVDtPhy", "xyz_DTempDtPhy", "xyzf_DQMixDtPhy"
!
! Leap-frog scheme is used as time integration.
! By default, semi-implicit scheme is applied to gravitational terms
! for extension of $ \Delta t $ .
! Explicit scheme can be applied to gravitational terms by changing
! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
!
! モジュール引用 ; USE statements
!
use constants, only: RPlanet, CpDry, Grav
! $ g $ [m s-2].
! 重力加速度.
! Gravitational acceleration
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: ncmax, IndexH2OVap, CompositionInqFlagAdv
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! 座標データ設定
! Axes data settings
!
use axesset, only: r_Sigma, z_Sigma ! $ \sigma $ レベル (整数).
! Full $ \sigma $ level
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
#ifdef LIB_MPI
#ifdef SJPACK
use wa_mpi_module_sjpack, only: wa_Lapla_wa, w_Lapla_w, wa_LaplaInv_wa, wa_DivLambda_xya => wa_DivLambda_xva, wa_DivMu_xya => wa_DivMu_xva, xy_GradLambda_w => xv_GradLambda_w, xya_GradLambda_wa => xva_GradLambda_wa, xy_GradMu_w => xv_GradMu_w, xya_GradMu_wa => xva_GradMu_wa, w_xy => w_xv, xy_w => xv_w, wa_xya => wa_xva, xya_wa => xva_wa
#else
use wa_mpi_module, only: wa_Lapla_wa, w_Lapla_w, wa_LaplaInv_wa, wa_DivLambda_xya => wa_DivLambda_xva, wa_DivMu_xya => wa_DivMu_xva, xy_GradLambda_w => xv_GradLambda_w, xya_GradLambda_wa => xva_GradLambda_wa, xy_GradMu_w => xv_GradMu_w, xya_GradMu_wa => xva_GradMu_wa, w_xy => w_xv, xy_w => xv_w, wa_xya => wa_xva, xya_wa => xva_wa
#endif
#elif AXISYMMETRY
use wa_zonal_module, only: wa_Lapla_wa, w_Lapla_w, wa_LaplaInv_wa, wa_DivLambda_xya, wa_DivMu_xya, xy_GradLambda_w, xya_GradLambda_wa, xy_GradMu_w, xya_GradMu_wa, w_xy, xy_w, wa_xya, xya_wa
#elif SJPACK
use wa_module_sjpack, only: wa_Lapla_wa, w_Lapla_w, wa_LaplaInv_wa, wa_DivLambda_xya, wa_DivMu_xya, xy_GradLambda_w, xya_GradLambda_wa, xy_GradMu_w, xya_GradMu_wa, w_xy, xy_w, wa_xya, xya_wa
#elif AXISYMMETRY_SJPACK
use wa_zonal_module_sjpack, only: wa_Lapla_wa, w_Lapla_w, wa_LaplaInv_wa, wa_DivLambda_xya, wa_DivMu_xya, xy_GradLambda_w, xya_GradLambda_wa, xy_GradMu_w, xya_GradMu_wa, w_xy, xy_w, wa_xya, xya_wa
#else
use wa_module, only: wa_Lapla_wa, w_Lapla_w, wa_LaplaInv_wa, wa_DivLambda_xya, wa_DivMu_xya, xy_GradLambda_w, xya_GradLambda_wa, xy_GradMu_w, xya_GradMu_wa, w_xy, xy_w, wa_xya, xya_wa
#endif
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 文字列操作
! Character handling
!
use dc_string, only: LChar
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: DP ! 倍精度実数型. Double precision.
! 積分と平均の操作
! Operation for integral and average
!
use intavr_operate, only: IntLonLat_xy
! 質量の補正
! Mass fixer
!
use mass_fixer, only: MassFixer, MassFixerColumn
! セミラグランジュ法による物質移流計算
! Semi-Lagrangian method for tracer transport
!
use sltt, only: SLTTMain
!!$ ! 物質移流 (有限体積法)
!!$ ! Tracer Transport (finite volume method)
!!$ !
!!$ use tt_fv, only: TTFVMain
!
! Utility module for advection test
!
use adv_test, only: AdvTestSetHorVels, AdvTestSetVerVel
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_UB (0:imax-1, 1:jmax, 1:kmax)
! $ u (t-\Delta t) $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_VB (0:imax-1, 1:jmax, 1:kmax)
! $ v (t-\Delta t) $ . 南北風速. Northward wind
real(DP), intent(in):: xyz_TempB (0:imax-1, 1:jmax, 1:kmax)
! $ T (t-\Delta t) $ . 温度. Temperature
real(DP), intent(in):: xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q (t-\Delta t) $ . 比湿. Specific humidity
real(DP), intent(in):: xy_PsB (0:imax-1, 1:jmax)
! $ p_s (t-\Delta t) $ . 地表面気圧. Surface pressure
real(DP), intent(in):: xyz_UN (0:imax-1, 1:jmax, 1:kmax)
! $ u (t) $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_VN (0:imax-1, 1:jmax, 1:kmax)
! $ v (t) $ . 南北風速. Northward wind
real(DP), intent(in):: xyz_TempN (0:imax-1, 1:jmax, 1:kmax)
! $ T (t) $ . 温度. Temperature
real(DP), intent(in):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q (t) $ . 比湿. Specific humidity
real(DP), intent(in):: xy_PsN (0:imax-1, 1:jmax)
! $ p_s (t) $ . 地表面気圧. Surface pressure
real(DP), intent(in):: xyz_DUDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{u}{t}\right)^{phy} $ .
! 外力項 (物理過程) による東西風速変化.
! Eastward wind tendency by external force terms (physical processes)
real(DP), intent(in):: xyz_DVDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{v}{t}\right)^{phy} $ .
! 外力項 (物理過程) による南北風速変化.
! Northward wind tendency by external force terms (physical processes)
real(DP), intent(in):: xyz_DTempDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{T}{t}\right)^{phy} $ .
! 外力項 (物理過程) による温度変化.
! Temperature tendency by external force terms (physical processes)
real(DP), intent(in):: xyzf_DQMixDtPhy(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \left(\DP{q}{t}\right)^{phy} $ .
! 外力項 (物理過程) による比湿変化.
! Temperature tendency by external force terms (physical processes)
real(DP), intent(in):: xy_SurfHeight (0:imax-1, 1:jmax)
! $ z_s $ . 地表面高度.
! Surface height.
real(DP), intent(out):: xyz_UA (0:imax-1, 1:jmax, 1:kmax)
! $ u (t+\Delta t) $ . 東西風速. Eastward wind
real(DP), intent(out):: xyz_VA (0:imax-1, 1:jmax, 1:kmax)
! $ v (t+\Delta t) $ . 南北風速. Northward wind
real(DP), intent(out):: xyz_TempA (0:imax-1, 1:jmax, 1:kmax)
! $ T (t+\Delta t) $ . 温度. Temperature
real(DP), intent(out):: xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q (t+\Delta t) $ . 比湿. Specific humidity
real(DP), intent(out):: xy_PsA (0:imax-1, 1:jmax)
! $ p_s (t+\Delta t) $ . 地表面気圧. Surface pressure
real(DP), intent(out), optional :: xyz_ArgOMG(0:imax-1, 1:jmax, 1:kmax)
!
! OMEGA, DP/Dt
! 作業変数
! Work variables
!
real(DP):: xyz_UCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ U (t) = u (t) \cos \varphi $ .
real(DP):: xyz_VCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ V (t) = v (t) \cos \varphi $ .
real(DP):: xyz_VorN (0:imax-1, 1:jmax, 1:kmax)
! $ \zeta (t) $ . 渦度. Vorticity
real(DP):: xyz_DivN (0:imax-1, 1:jmax, 1:kmax)
! $ D (t) $ . 発散. Divergence
real(DP):: xy_PiN (0:imax-1, 1:jmax)
! $ \pi = \ln p_s $
real(DP):: wz_DivN (lmax, 1:kmax)
! $ D (t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP):: wz_TempN (lmax, 1:kmax)
! $ T (t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP):: w_PiN (lmax)
! $ \pi $ スペクトル
real(DP):: xy_GradLambdaPiN (0:imax-1, 1:jmax)
! $ \DP{\pi}{\lambda} $
real(DP):: xy_GradMuPiN (0:imax-1, 1:jmax)
! $ (1-\mu^2) \DP{\pi}{\mu} $
real(DP):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
! $ \Dvect{v} \cdot \nabla \pi $ .
! $ \pi $ の移流. Advection of $ \pi $
real(DP):: xyz_UAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ U_A (t) $ . 東西運動量移流項.
! Eastward advection of momentum
real(DP):: xyz_VAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ V_A (t) $ . 南北運動量移流項.
! Northward advection of momentum
real(DP):: xyz_TempNonLinearN (0:imax-1, 1:jmax, 1:kmax)
! $ H (t) $ . 温度時間変化項.
! Temperature tendency
real(DP):: xyzf_QMixNonLinearN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ R (t) $ . 比湿時間変化項.
! Specific humidity tendency
real(DP):: xyz_KinEngyN (0:imax-1, 1:jmax, 1:kmax)
! $ KE (t) $ . 運動エネルギー項.
! Kinetic energy
real(DP):: xyz_TempUAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ UT' (t) $ . 温度東西移流項.
! Eastward advection of temperature
real(DP):: xyz_TempVAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ VT' (t) $ . 温度南北移流項.
! Northward advection of temperature
real(DP):: xyr_SigDotN (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} $ .
! 鉛直流. Vertical flow
real(DP):: xy_DPiDtNG (0:imax-1, 1:jmax)
! $ Z $ . 地表面気圧時間変化非重力波項.
! Non-gravity wave component of surface pressure tendency
real(DP):: xyzf_QMixUAdvN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ Uq (t) $ . 比湿東西移流項.
! Eastward advection of specific humidity
real(DP):: xyzf_QMixVAdvN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ Vq (t) $ . 比湿南北移流項.
! Northward advection of specific humidity
real(DP):: wz_DVorDtNG (lmax, 1:kmax)
! $ \left( \DD{\zeta}{t} (t) \right)^{NG}$ . 渦度変化の非重力波成分 (スペクトル).
! Non-gravity wave component of vorticity tendency (spectral)
real(DP):: wz_DDivDtNG (lmax, 1:kmax)
! $ \left( \DD{D}{t} (t) \right)^{NG}$ . 発散変化の非重力波成分 (スペクトル).
! Non-gravity wave component of divergence tendency (spectral)
real(DP):: wz_DTempDtNG (lmax, 1:kmax)
! $ \left( \DD{T}{t} (t) \right)^{NG}$ . 温度変化の非重力波成分 (スペクトル).
! Non-gravity wave component of temperature tendency (spectral)
real(DP):: wzf_DQMixDtN(lmax, 1:kmax, 1:ncmax)
! $ \DD{q}{t} (t) $ . 比湿変化 (スペクトル).
! Specific humidity tendency (spectral)
real(DP):: w_DPiDtNG(lmax)
! $ \left( \DD{p_s}{t} (t) \right)^{NG}$ . 地表面気圧変化費重力波項 (スペクトル).
! Non-gravity wave component of surface pressure tendency (spectral)
real(DP):: xyz_UCosLatB (0:imax-1, 1:jmax, 1:kmax)
! $ U (t-\Delta t) = u (t-\Delta t) \cos \varphi $ .
real(DP):: xyz_VCosLatB (0:imax-1, 1:jmax, 1:kmax)
! $ V (t-\Delta t) = v (t-\Delta t) \cos \varphi $ .
real(DP):: wz_VorB (lmax, 1:kmax)
! $ \zeta (t-\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP):: wz_DivB (lmax, 1:kmax)
! $ D (t-\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP):: wz_TempB (lmax, 1:kmax)
! $ T (t-\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP):: w_PiB (lmax)
! $ \pi = \ln p_s (t-\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP):: wzf_QMixB(lmax, 1:kmax, 1:ncmax)
! $ q (t-\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP):: xyz_DUDtPhyCosLat (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{u}{t}\right)^{phy} \cos \varphi $ .
real(DP):: xyz_DVDtPhyCosLat (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{v}{t}\right)^{phy} \cos \varphi $ .
real(DP):: wz_VorA (lmax, 1:kmax)
! $ \zeta (t+\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP):: wz_DivA (lmax, 1:kmax)
! $ D (t+\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP):: wz_TempA (lmax, 1:kmax)
! $ T (t+\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP):: w_PiA (lmax)
! $ \pi = \ln p_s (t+\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP):: wzf_QMixA(lmax, 1:kmax, 1:ncmax)
! $ q (t+\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP):: wz_Psi (lmax, 1:kmax)
! $ \psi $ . 流線関数. Streamline function
real(DP):: wz_Chi (lmax, 1:kmax)
! $ \chi $ . ポテンシャル. Potential
real(DP):: xyz_UCosLatA (0:imax-1, 1:jmax, 1:kmax)
! $ U (t+\Delta t) = u (t+\Delta t) \cos \varphi $ .
real(DP):: xyz_VCosLatA (0:imax-1, 1:jmax, 1:kmax)
! $ V (t+\Delta t) = v (t+\Delta t) \cos \varphi $ .
real(DP):: wz_VorDiffA (lmax, 1:kmax)
! $ \mathscr{D}(\zeta) $ .
! 運動量水平拡散による渦度変化 (スペクトル).
! Vorticity tendency by
! horizontal momentum diffusion (spectral)
real(DP):: wz_DivDiffA (lmax, 1:kmax)
! $ \mathscr{D}(D) $ .
! 運動量水平拡散による発散変化 (スペクトル).
! Divergence tendency by
! horizontal momentum diffusion (spectral)
real(DP):: wz_PsiDiff (lmax, 1:kmax)
! 運動量水平拡散による流線関数 $ \psi $ 変化
! Streamline function tendency by
! horizontal momentum diffusion
real(DP):: wz_ChiDiff (lmax, 1:kmax)
! 運動量水平拡散によるポテンシャル $ \chi $ 変化
! Potential tendency by
! horizontal momentum diffusion
real(DP):: xyz_UDiff (0:imax-1, 1:jmax, 1:kmax)
! 運動量水平拡散による東西風変化.
! Eastward wind tendency by
! horizontal momentum diffusion
real(DP):: xyz_VDiff (0:imax-1, 1:jmax, 1:kmax)
! 運動量水平拡散による南北風変化.
! Northward wind tendency by
! horizontal momentum diffusion
! 地表面高度関連
! Surface height etc.
!
real(DP):: w_SurfGeoPot (lmax)
! $ \Phi_s $ . 地表ジオポテンシャル.
! Surface geo-potential
! Variables for mass fixing
!
real(DP):: MassB
real(DP):: MassA
real(DP):: xyr_PressA(0:imax-1, 1:jmax, 0:kmax)
real(DP):: xyr_PressB(0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyz_PressA (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_HDifPsA (0:imax-1, 1:jmax)
real(DP) :: xyz_HDifPressA (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyzf_DQMixDtHDCor(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
integer :: kp, kn
real(DP) :: xyz_OMG (0:imax-1, 1:jmax, 1:kmax)
!
! OMEGA, DP/Dt
real(DP) :: xyzf_QMixTentative(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
logical :: FlagExistNonAdvComp
real(DP) :: xyz_UTracerAdv(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_VTracerAdv(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyr_AdvTestStreamFunc(0:imax-1, 1:jmax, 0:kmax)
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. dynamics_hspl_vas83_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! TimeIntegration で使用する係数の設定
! Configure coefficients for "TimeIntegration"
!
call SemiImplMatrix
! 格子点値をスペクトル値へ ( $ t-\Delta t$ )
! Exchange grid values to spectral values ( $ t-\Delta t$ )
!
! 渦度発散の計算
! Calculate vorticity and divergence
do k = 1, kmax
xyz_UCosLatB(:,:,k) = xyz_UB(:,:,k) * xy_CosLat
xyz_VCosLatB(:,:,k) = xyz_VB(:,:,k) * xy_CosLat
end do
wz_VorB = ( wa_DivLambda_xya( xyz_VCosLatB ) - wa_DivMu_xya ( xyz_UCosLatB ) ) / RPlanet
wz_DivB = ( wa_DivLambda_xya( xyz_UCosLatB ) + wa_DivMu_xya ( xyz_VCosLatB ) ) / RPlanet
!
wz_TempB = wa_xya( xyz_TempB )
!
w_PiB = w_xy( log( xy_PsB ) )
! 格子点値をスペクトル値へ ( $ t $ )
! Exchange grid values to spectral values ( $ t $ )
!
! 渦度発散の計算
! Calculate vorticity and divergence
!
do k = 1, kmax
xyz_UCosLatN(:,:,k) = xyz_UN(:,:,k) * xy_CosLat
xyz_VCosLatN(:,:,k) = xyz_VN(:,:,k) * xy_CosLat
end do
xyz_VorN = xya_wa( wa_DivLambda_xya( xyz_VCosLatN ) - wa_DivMu_xya ( xyz_UCosLatN ) ) / RPlanet
wz_DivN = ( wa_DivLambda_xya( xyz_UCosLatN ) + wa_DivMu_xya ( xyz_VCosLatN ) ) / RPlanet
xyz_DivN = xya_wa( wz_DivN )
!
select case ( IDTimeIntegScheme )
case ( IDTimeIntegSchemeExplicit )
wz_TempN = wa_xya( xyz_TempN )
end select
!
xy_PiN = log( xy_PsN )
w_PiN = w_xy( xy_PiN )
xy_GradLambdaPiN = xy_GradLambda_w( w_PiN ) / RPlanet
xy_GradMuPiN = xy_GradMu_w ( w_PiN ) / RPlanet
! 地表ジオポテンシャルの計算
! Calculate surface geo-potential
!
w_SurfGeoPot = w_xy( Grav * xy_SurfHeight )
! 格子点上での非線形力学項の計算
! Calculate non-linear dynamical terms on grid points
!
call NonLinearOnGrid( xyz_UCosLatN, xyz_VCosLatN, xyz_VorN, xyz_DivN, xyz_TempN, xyzf_QMixN(:,:,:,IndexH2OVap), xy_GradLambdaPiN, xy_GradMuPiN, xyz_PiAdv, xyz_UAdvN, xyz_VAdvN, xyz_TempNonLinearN, xyz_KinEngyN, xyz_TempUAdvN, xyz_TempVAdvN, xyr_SigDotN, xy_DPiDtNG )
! 非線形項と物理過程による時間変化項の和をスペクトル変換
! Spectral transformation of sum of non-linear term and tendencies by physical
! processes
!
do k = 1, kmax
xyz_DUDtPhyCosLat(:,:,k) = xyz_DUDtPhy(:,:,k) * xy_CosLat
xyz_DVDtPhyCosLat(:,:,k) = xyz_DVDtPhy(:,:,k) * xy_CosLat
end do
wz_DVorDtNG = ( wa_DivLambda_xya( xyz_VAdvN + xyz_DVDtPhyCosLat ) - wa_DivMu_xya ( xyz_UAdvN + xyz_DUDtPhyCosLat ) ) / RPlanet
wz_DDivDtNG = ( wa_DivLambda_xya( xyz_UAdvN + xyz_DUDtPhyCosLat ) + wa_DivMu_xya ( xyz_VAdvN + xyz_DVDtPhyCosLat ) ) / RPlanet - wa_Lapla_wa( wa_xya(xyz_KinEngyN) ) / RPlanet**2
wz_DTempDtNG = - ( wa_DivLambda_xya( xyz_TempUAdvN ) + wa_DivMu_xya ( xyz_TempVAdvN ) ) / RPlanet + wa_xya( xyz_TempNonLinearN + xyz_DTempDtPhy )
w_DPiDtNG = w_xy( xy_DPiDtNG )
! 時間積分
! Time integration
!
call TimeIntegration( w_SurfGeoPot, wz_DVorDtNG, wz_DDivDtNG, wz_DTempDtNG, w_DPiDtNG, wz_DivN, wz_TempN, w_PiN, wz_VorB, wz_DivB, wz_TempB, w_PiB, wz_VorA, wz_DivA, wz_TempA, w_PiA )
! Divergence damping
!
call DivergenceDamping( wz_DivA )
! スペクトル値を格子点値へ ( $ t+\Delta t$ )
! Exchange spectral values to grid values ( $ t+\Delta t$ )
!
wz_Psi = wa_LaplaInv_wa( wz_VorA ) * RPlanet**2
wz_Chi = wa_LaplaInv_wa( wz_DivA ) * RPlanet**2
xyz_UCosLatA = ( xya_GradLambda_wa( wz_Chi ) - xya_GradMu_wa ( wz_Psi ) ) / RPlanet
xyz_VCosLatA = ( xya_GradLambda_wa( wz_Psi ) + xya_GradMu_wa ( wz_Chi ) ) / RPlanet
do k = 1, kmax
xyz_UA(:,:,k) = xyz_UCosLatA(:,:,k) / xy_CosLat
xyz_VA(:,:,k) = xyz_VCosLatA(:,:,k) / xy_CosLat
end do
xyz_TempA = xya_wa( wz_TempA )
xy_PsA = exp( xy_w( w_PiA ) )
! 水平拡散とスポンジ層による散逸の補正
! Correction for dissipation of kinetic enery by horizontal diffusion and sponge
! layer
!
! 水平拡散とスポンジ層による渦度発散の時間変化
! Vorticity and divergence tendency by horizontal diffusion and damping in a sponge
! layer
!
wz_VorDiffA = wz_VorA * wz_DisCoefM
wz_DivDiffA = wz_DivA * wz_DisCoefM
! 水平拡散とスポンジ層による散逸の摩擦熱補正
! Correction for internal energy by horizontal diffusion and damping in a sponge
! layer
!
wz_PsiDiff = wa_LaplaInv_wa( wz_VorDiffA ) * RPlanet**2
wz_ChiDiff = wa_LaplaInv_wa( wz_DivDiffA ) * RPlanet**2
xyz_UDiff = ( xya_GradLambda_wa( wz_ChiDiff ) - xya_GradMu_wa ( wz_PsiDiff ) ) / RPlanet
xyz_VDiff = ( xya_GradLambda_wa( wz_PsiDiff ) + xya_GradMu_wa ( wz_ChiDiff ) ) / RPlanet
do k = 1, kmax
xyz_TempA(:,:,k) = xyz_TempA(:,:,k) - ( xyz_UA(:,:,k) * xyz_UDiff(:,:,k) + xyz_VA(:,:,k) * xyz_VDiff(:,:,k) ) / xy_CosLat / CpDry * 2.0_DP * DelTime
end do
! Mass fixer
! Total Mass
! NOTICE: It is assumed that the total atmospheric mass does not change.
!
if ( FlagTotMassFix ) then
MassB = IntLonLat_xy( xy_PsB )
MassA = IntLonLat_xy( xy_PsA )
if ( MassA /= 0.0_DP ) then
xy_PsA = MassB / MassA * xy_PsA
end if
end if
! Set velocities for tracer advection
!
xyz_UTracerAdv = xyz_UN
xyz_VTracerAdv = xyz_VN
! 主に移流テストのために使う if 文
!
if ( FlagAdvTest ) then
xyz_UA = xyz_UB
xyz_VA = xyz_VB
xyz_TempA = xyz_TempB
xy_PsA = xy_PsB
xyr_SigDotN = 0.0_DP
!
! Utility module for advection test
!
call AdvTestSetHorVels( xyz_UTracerAdv, xyz_VTracerAdv )
! calculation of variables required in spectral advection calculation
do k = 1, kmax
xyz_UCosLatN(:,:,k) = xyz_UTracerAdv(:,:,k) * xy_CosLat
xyz_VCosLatN(:,:,k) = xyz_VTracerAdv(:,:,k) * xy_CosLat
end do
wz_DivN = ( wa_DivLambda_xya( xyz_UCosLatN ) + wa_DivMu_xya ( xyz_VCosLatN ) ) / RPlanet
xyz_DivN = xya_wa( wz_DivN )
call AdvTestSetVerVel( xyr_SigDotN, xyr_AdvTestStreamFunc )
call HistoryAutoPut( TimeN, 'AdvTestStreamFunc', xyr_AdvTestStreamFunc )
end if
if ( FlagSLTT ) then
!!$ select case ( IDTTMethod )
!!$ case ( IDTTMethodSL )
! Semi-Lagrangian advection
do k = 0, kmax
xyr_PressA(:,:,k) = xy_PsA * r_Sigma(k)
xyr_PressB(:,:,k) = xy_PsB * r_Sigma(k)
end do
! Mass fixer is included in SLTTMain
call SLTTMain( xyr_PressB, xyr_PressA, xyz_UTracerAdv, xyz_VTracerAdv, xyr_SigDotN, xyzf_DQMixDtPhy, xyzf_QMixB, xyzf_QMixA )
!!$ case ( IDTTMethodFV )
!!$ ! Finite volume method
!!$
!!$ do k = 0, kmax
!!$ xyr_PressA(:,:,k) = xy_PsA * r_Sigma(k)
!!$ xyr_PressB(:,:,k) = xy_PsB * r_Sigma(k)
!!$ end do
!!$ call TTFVMain( &
!!$ & xyr_PressB, xyr_PressA, &
!!$ & xyz_UTracerAdv, xyz_VTracerAdv, xyr_SigDotN, & ! (in)
!!$ & xyzf_DQMixDtPhy, & ! (in)
!!$ & xyzf_QMixB, & ! (in)
!!$ & xyzf_QMixA & ! (out)
!!$ & )
else
!!$ case ( IDTTMethodSp )
! Spectral advection
do n = 1, ncmax
wzf_QMixB(:,:,n) = wa_xya( xyzf_QMixB(:,:,:,n) )
end do
!
call NonLinearOnGridQMix( xyz_UCosLatN, xyz_VCosLatN, xyz_DivN, xyr_SigDotN, xyzf_QMixN, xyzf_QMixNonLinearN, xyzf_QMixUAdvN, xyzf_QMixVAdvN )
do n = 1, ncmax
wzf_DQMixDtN(:,:,n) = - ( wa_DivLambda_xya( xyzf_QMixUAdvN(:,:,:,n) ) + wa_DivMu_xya ( xyzf_QMixVAdvN(:,:,:,n) ) ) / RPlanet + wa_xya( xyzf_QMixNonLinearN(:,:,:,n) + xyzf_DQMixDtPhy (:,:,:,n) )
end do
call TimeIntegrationQMix( wzf_DQMixDtN, wzf_QMixB, wzf_QMixA )
do n = 1, ncmax
xyzf_QMixA(:,:,:,n) = xya_wa( wzf_QMixA(:,:,n) )
end do
if ( FlagMassHorDifCor ) then
! Correction for horizontal diffusion of constituents.
! This is performed for the aim of correcting the diffusion in the viscinity of
! steep terrain slope.
!
do k = 1, kmax
xyz_PressA(:,:,k) = xy_PsA * z_Sigma(k)
end do
xy_HDifPsA = xy_w( w_DisCoefQ * w_xy( xy_PsA ) )
do k = 1, kmax
xyz_HDifPressA(:,:,k) = z_Sigma(k) * xy_HDifPsA(:,:)
end do
do n = 1, ncmax
do k = 1, kmax
kp = k-1
kn = k+1
if( k == 1 ) kp = 1
if( k == kmax ) kn = kmax
xyzf_DQMixDtHDCor(:,:,k,n) = - ( xyzf_QMixA(:,:,kn,n) - xyzf_QMixA(:,:,kp,n) ) / ( xyz_PressA(:,:,kn) - xyz_PressA(:,:,kp) ) * xyz_HDifPressA(:,:,k)
end do
end do
xyzf_QMixA = xyzf_QMixA + xyzf_DQMixDtHDCor * 2.0_DP * DelTime
endif
! Mass fixer
! Constituents
!
do k = 0, kmax
xyr_PressA(:,:,k) = xy_PsA * r_Sigma(k)
xyr_PressB(:,:,k) = xy_PsB * r_Sigma(k)
end do
call MassFixer( xyr_PressA, xyzf_QMixA, xyr_PressRef = xyr_PressB, xyzf_QMixRef = ( xyzf_QMixB+xyzf_DQMixDtPhy*2.0_DP*DelTime ) )
!!$ end select
end if
! 主に移流テストのために使う if 文
!
!!$ if ( FlagAdvTest ) then
!!$ xyr_SigDotN = 0.0_DP
!!$ end if
! Calculation for non-advection case
! オイラー移流計算のとき、ピーキーな分布の物質を移流させない。計算安定の
! ため。
!
FlagExistNonAdvComp = .false.
do n = 1, ncmax
if ( .not. CompositionInqFlagAdv( n ) ) then
xyzf_QMixA(:,:,:,n) = xyzf_QMixB(:,:,:,n) + xyzf_DQMixDtPhy(:,:,:,n) * 2.0_DP * DelTime
xyzf_QMixTentative(:,:,:,n) = xyzf_QMixA(:,:,:,n)
FlagExistNonAdvComp = .true.
else
xyzf_QMixTentative(:,:,:,n) = 0.0_DP
end if
end do
if ( FlagExistNonAdvComp ) then
! Mass fixer
! Constituents
call MassFixerColumn( xyr_PressB, xyzf_QMixTentative, xyr_PressRef = xyr_PressB, xyzf_QMixRef = ( xyzf_QMixB+xyzf_DQMixDtPhy*2.0_DP*DelTime ) )
do n = 1, ncmax
if ( .not. CompositionInqFlagAdv( n ) ) then
xyzf_QMixA(:,:,:,n) = xyzf_QMixTentative(:,:,:,n)
end if
end do
end if
! ヒストリデータ出力
! History data output
!
call OutputDiagnosedVariables( xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA, xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy, xyr_SigDotN, xy_DPiDtNG, xyz_PiAdv, xyz_OMG )
if ( present( xyz_ArgOMG ) ) then
xyz_ArgOMG = xyz_OMG
end if
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine DynamicsHSplVAS83
| Subroutine : |
モジュール内部の変数の割り付け解除を行います.
Deallocate variables in this module.
subroutine DynamicsHSplVAS83Finalize
!
! モジュール内部の変数の割り付け解除を行います.
!
! Deallocate variables in this module.
!
! 宣言文 ; Declaration statements
!
implicit none
! 実行文 ; Executable statement
!
if ( .not. dynamics_hspl_vas83_inited ) return
! デフォルト値へ戻す
! Return to default values
!
DelTimeSave = - 1.0_DP
! 割り付け解除
! Deallocation
!
if ( allocated( xy_SinLat ) ) deallocate( xy_SinLat )
if ( allocated( xy_CosLat ) ) deallocate( xy_CosLat )
if ( allocated( xy_Cori ) ) deallocate( xy_Cori )
if ( allocated( z_HydroAlpha ) ) deallocate( z_HydroAlpha )
if ( allocated( z_HydroBeta ) ) deallocate( z_HydroBeta )
if ( allocated( z_TInpCoefA ) ) deallocate( z_TInpCoefA )
if ( allocated( z_TInpCoefB ) ) deallocate( z_TInpCoefB )
if ( allocated( z_TInpCoefK ) ) deallocate( z_TInpCoefK )
if ( allocated( z_RefTemp ) ) deallocate( z_RefTemp )
if ( allocated( r_RefTemp ) ) deallocate( r_RefTemp )
if ( allocated( w_LaplaEigVal ) ) deallocate( w_LaplaEigVal )
if ( allocated( wz_DisCoefM ) ) deallocate( wz_DisCoefM )
if ( allocated( wz_DisCoefH ) ) deallocate( wz_DisCoefH )
if ( allocated( w_DisCoefQ ) ) deallocate( w_DisCoefQ )
if ( allocated( z_siMtxC ) ) deallocate( z_siMtxC )
if ( allocated( z_siMtxG ) ) deallocate( z_siMtxG )
if ( allocated( zz_siMtxH ) ) deallocate( zz_siMtxH )
if ( allocated( zz_siMtxDiH ) ) deallocate( zz_siMtxDiH )
if ( allocated( wzz_siMtxWDiH ) ) deallocate( wzz_siMtxWDiH )
if ( allocated( zz_siMtxGCt ) ) deallocate( zz_siMtxGCt )
if ( allocated( zz_siMtxW ) ) deallocate( zz_siMtxW )
if ( allocated( zz_siMtxQ ) ) deallocate( zz_siMtxQ )
if ( allocated( zz_siMtxS ) ) deallocate( zz_siMtxS )
if ( allocated( zz_siMtxR ) ) deallocate( zz_siMtxR )
!!$ if ( allocated(nmo ) ) deallocate( nmo )
!!$ if ( allocated(wzz_siMtxM ) ) deallocate( wzz_siMtxM )
!!$ if ( allocated(z_siMtxPivWork) ) deallocate( z_siMtxPivWork )
if ( allocated(wzz_siMtxLU ) ) deallocate( wzz_siMtxLU )
if ( allocated(wz_siMtxPiv ) ) deallocate( wz_siMtxPiv )
dynamics_hspl_vas83_inited = .false.
end subroutine DynamicsHSplVAS83Finalize
| Subroutine : |
計算に必要なパラメタの設定や NAMELIST#dynamics_hspl_vas83_nml の読み込みを行います.
Configure parameters for calculation, and load "NAMELIST#dynamics_hspl_vas83_nml"
This procedure input/output NAMELIST#dynamics_hspl_vas83_nml .
subroutine DynamicsHSplVAS83Init
!
! 計算に必要なパラメタの設定や NAMELIST#dynamics_hspl_vas83_nml
! の読み込みを行います.
!
! Configure parameters for calculation,
! and load "NAMELIST#dynamics_hspl_vas83_nml"
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: RPlanet, Omega, GasRDry, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 格子点設定
! Grid points settings
!
use gridset, only: nmax, lmax, imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: ncmax, a_QMixName
! 座標データ設定
! Axes data settings
!
use axesset, only: z_Sigma, r_Sigma, z_DelSigma, AxNameX, AxNameY, AxNameZ, AxNameR, AxNameT
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
#ifdef LIB_MPI
#ifdef SJPACK
use wa_mpi_module_sjpack, only: xy_Lat => xv_Lat, w_xy => w_xv, rn, nm_l
#else
use wa_mpi_module, only: xy_Lat => xv_Lat, w_xy => w_xv, rn, nm_l
#endif
#elif AXISYMMETRY
use wa_zonal_module, only: xy_Lat, w_xy, rn, nm_l
#elif SJPACK
use wa_module_sjpack, only: xy_Lat, w_xy, rn, nm_l
#elif AXISYMMETRY_SJPACK
use wa_zonal_module_sjpack, only: xy_Lat, w_xy, rn, nm_l
#else
use wa_module, only: xy_Lat, w_xy, rn, nm_l
#endif
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoAddVariable
! gtool4 データ入力
! Gtool4 data input
!
use gtool_history, only: HistoryGet
! ファイル入出力補助
! File I/O support
!
use dc_iounit, only: FileOpen
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: STDOUT, STRING ! 文字列. Strings.
! 日付および時刻の取り扱い
! Date and time handler
!
use dc_calendar, only: DCCalConvertByUnit
! 組み込み関数 PRESENT の拡張版関数
! Extended functions of intrinsic function "PRESENT"
!
use dc_present, only: present_and_not_empty
! 文字列操作
! Character handling
!
use dc_string, only: LChar
! 時刻管理
! Time control
!
use timeset, only: TimeN ! ステップ $ t $ の時刻. Time of step $ t $.
! 質量の補正
! Mass fixer
!
use mass_fixer, only : MassFixerInit
! セミラグランジュ法による物質移流計算
! Semi-Lagrangian method for tracer transport
use sltt, only: SLTTInit
!!$ ! 物質移流 (有限体積法)
!!$ ! Tracer Transport (finite volume method)
!!$ !
!!$ use tt_fv, only: TTFVInit
!
! Utility module for advection test
!
use adv_test, only: AdvTestInit
! 宣言文 ; Declaration statements
!
implicit none
! 基準温度の設定のための作業変数
! Work variable for reference temperature
!
character(TOKEN) :: TimeIntegScheme
! 時間積分法.
! 以下の方法を選択可能.
!
! Time integration scheme.
! Available schemes are as follows.
!
! * "Semi-implicit"
! * "Explicit"
real(DP):: RefTemp
! $ \overline{T} $ . 基準温度.
! Reference temperature
! 水平拡散, スポンジ層のための作業変数
! Work variable for horizontal diffusion and sponge layer
!
real(DP) :: w_HDifCoefM (lmax)
! $ {\cal D}_M = -K_{HD} [(-1)^{N_D/2}\nabla^{N_D}- (2/a^2)^{N_D/2}] $ .
! 運動量水平拡散係数.
! Coefficient of horizontal momentum diffusion
real(DP) :: w_HDifCoefH (lmax)
! $ {\cal D}_H = -(-1)^{N_D/2}K_{HD}\nabla^{N_D} $ .
! 熱, 水水平拡散係数.
! Coefficient of horizontal thermal and water diffusion
real(DP) :: wz_SpongeLayerCoefM(lmax, 1:kmax)
! スポンジ層における運動量の減衰係数
! Damping coefficient for momentum in a sponge layer
real(DP) :: wz_SpongeLayerCoefH(lmax, 1:kmax)
! スポンジ層における熱の減衰係数
! Damping coefficient for temperature zonal perturbation
! in a sponge layer
integer:: HDOrder ! 超粘性の次数. Order of hyper-viscosity
real(DP):: VisCoef ! 超粘性係数. Hyper-viscosity coefficient
real(DP):: HDEFoldTimeValue ! 最大波数に対する e-folding time.
! 負の値を与えると, 水平拡散係数をゼロにします.
!
! E-folding time for maximum wavenumber.
! If negative value is given,
! coefficients of horizontal diffusion become zero.
character(TOKEN):: HDEFoldTimeUnit
! 最大波数に対する e-folding time の単位.
! Unit of e-folding time for maximum wavenumber
real(DP):: HDEFoldTime ! 最大波数に対する e-folding time [単位: 秒].
! E-folding time for maximum wavenumber [Unit: sec]
logical :: FlagSpongeLayer
logical :: FlagSpongeLayerforZonalMean
logical :: FlagSpongeLayerforHeat
real(DP) :: SLEFoldTimeValue
! スポンジ層の最上層における減衰時定数
! Damping time constant at the top model layer
! in a sponge layer
character(TOKEN) :: SLEFoldTimeUnit
! SLEFoldTimeValue の単位
! Unit of SLEFoldTimeValue
real(DP) :: SLEFoldTime
! スポンジ層の最上層における減衰時定数 (秒)
! Damping time constant at the top model layer
! in a sponge layer (sec)
integer :: SLOrder
! スポンジ層の減衰係数の sigma 依存性のオーダ
! Sigma dependence of damping coefficient
! in a sponge layer
integer :: SLNumLayer
! スポンジ層が適応されるモデル上端からの層の数
! Number of layers which the sponge layer is applied.
integer :: a_DegOrd(2)
real(DP) :: DivDampPeriodValue
! Period for divergence damping application
character(TOKEN) :: DivDampPeriodUnit
! Unit of DivDampPeriodValue
character(TOKEN):: TracerTransportMethod
! NonLinearOnGrid 等で使用する係数の設定のための作業変数
! Work variable for coefficients for "NonLinearOnGrid", etc.
!
real(DP):: Kappa ! $ \kappa = R/C_p $ .
! 気体定数の定圧比熱に対する比.
! Ratio of gas constant to specific heat
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n
integer:: l ! 波数方向に回る DO ループ用作業変数
! Work variables for DO loop in wavenumber
integer:: unit_nml ! NAMELIST ファイルオープン用装置番号.
! Unit number for NAMELIST file open
integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT.
! IOSTAT of NAMELIST read
! NAMELIST 変数群
! NAMELIST group name
!
namelist /dynamics_hspl_vas83_nml/ TimeIntegScheme, HDOrder, HDEFoldTimeValue, HDEFoldTimeUnit, FlagSpongeLayer, FlagSpongeLayerforZonalMean, FlagSpongeLayerforHeat, SLEFoldTimeValue, SLEFoldTimeUnit, SLOrder, SLNumLayer, RefTemp, FlagDivDamp, DivDampPeriodValue, DivDampPeriodUnit, DivDampTime0, FlagTotMassFix, FlagMassHorDifCor, FlagSLTT, FlagAdvTest
!
! デフォルト値については初期化手続 "dynamics_hspl_vas83#DynamicsInit"
! のソースコードを参照のこと.
!
! Refer to source codes in the initialization procedure
! "dynamics_hspl_vas83#DynamicsInit" for the default values.
!
! 実行文 ; Executable statement
!
if ( dynamics_hspl_vas83_inited ) return
! デフォルト値の設定
! Default values settings
!
TimeIntegScheme = 'Semi-implicit'
HDOrder = 8
HDEFoldTimeValue = 8640.0_DP
HDEFoldTimeUnit = 'sec'
FlagSpongeLayer = .false.
FlagSpongeLayerforZonalMean = .false.
FlagSpongeLayerforHeat = .false.
SLEFoldTimeValue = 86400.0_DP
SLEFoldTimeUnit = 'sec'
SLOrder = 1
SLNumLayer = kmax
RefTemp = 300.
FlagDivDamp = .false.
DivDampPeriodValue = 2.0_DP
DivDampPeriodUnit = 'day'
DivDampTime0 = 0.0_DP
FlagTotMassFix = .true.
FlagMassHorDifCor = .false.
!!$ TracerTransportMethod = "spectrum"
FlagSLTT = .false.
FlagAdvTest = .false.
! NAMELIST の読み込み
! NAMELIST is input
!
if ( trim(namelist_filename) /= '' ) then
call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)
rewind( unit_nml )
read( unit_nml, nml = dynamics_hspl_vas83_nml, iostat = iostat_nml ) ! (out)
close( unit_nml )
call NmlutilMsg( iostat_nml, module_name ) ! (in)
if ( iostat_nml == 0 ) write( STDOUT, nml = dynamics_hspl_vas83_nml )
end if
! 時間積分法のチェック
! Check time integration scheme
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')
IDTimeIntegScheme = IDTimeIntegSchemeSemiImplicit
case ('explicit')
IDTimeIntegScheme = IDTimeIntegSchemeExplicit
case default
call MessageNotify( 'E', module_name, '"TimeIntegScheme" must be "Semi-Implicit" or "Explicit".' )
end select
!
! Check tracer transport method
!
!!$ select case ( LChar( trim( TracerTransportMethod ) ) )
!!$ case ('spectrum')
!!$ IDTTMethod = IDTTMethodSp
!!$ case ('semi-lagrange')
!!$ IDTTMethod = IDTTMethodSL
!!$ case ('finite volume')
!!$ IDTTMethod = IDTTMethodFV
!!$ case default
!!$ call MessageNotify( 'E', module_name, &
!!$ & '"TracerTransportMethod" must be "spectrum", "semi-lagrange", or "finite volume".' )
!!$ end select
! SemiImplMatrix サブルーチン用に $ \Delta t $ の保存値に初期値を設定.
! Configure initial value to saved value of $ \Delta t $ for a subroutine "SemiImplMatrix"
!
DelTimeSave = - 1.0_DP
! NonLinearOnGrid 等で使用する係数の設定
! Configure coefficients for "NonLinearOnGrid", etc.
!
! $ \sin \varphi $ と $ \cos \varphi $ の計算
! Calculate $ \sin \varphi $ and $ \cos \varphi $
!
allocate( xy_SinLat (0:imax-1, 1:jmax) )
allocate( xy_CosLat (0:imax-1, 1:jmax) )
xy_SinLat = sin( xy_Lat )
xy_CosLat = cos( xy_Lat )
! コリオリパラメータの計算
! Calculate Coriolis parameter
!
allocate( xy_Cori (0:imax-1, 1:jmax) )
xy_Cori = 2.0_DP * Omega * xy_SinLat
! 静水圧の式の係数 $ \alpha $ , $ \beta $ の計算
! Calculate coefficients $ \alpha $ and $ \beta $ in hydrostatic equation.
!
allocate( z_HydroAlpha(1:kmax) )
allocate( z_HydroBeta (1:kmax) )
Kappa = GasRDry / CpDry
do k = 1, kmax
z_HydroAlpha(k) = ( r_Sigma(k-1) / z_Sigma(k) ) ** Kappa - 1.0_DP
z_HydroBeta(k) = 1.0_DP - ( r_Sigma(k) / z_Sigma(k) ) ** Kappa
enddo
! 温度鉛直補間の係数 $ \kappa $, $ a $ , $ b $ の計算
! Calculate coefficients $ \kappa $, $ a $ , $ b $
! for interpolation of temperature
!
allocate( z_TInpCoefA (1:kmax) )
allocate( z_TInpCoefB (1:kmax) )
allocate( z_TInpCoefK (1:kmax) )
do k = 1, kmax
z_TInpCoefK(k) = ( r_Sigma(k-1) * z_HydroAlpha(k) + r_Sigma(k ) * z_HydroBeta (k) ) / z_DelSigma(k)
enddo
z_TInpCoefA = 0.
do k = 2, kmax
z_TInpCoefA(k) = z_HydroAlpha(k) * ( 1.0_DP - (z_Sigma(k) / z_Sigma(k-1)) ** Kappa ) ** ( -1 )
end do
z_TInpCoefB = 0.
do k = 1, kmax - 1
z_TInpCoefB(k) = z_HydroBeta(k) * ( (z_Sigma(k) / z_Sigma(k+1) ) ** Kappa - 1.0_DP ) ** ( -1 )
end do
! 基準温度 (半整数レベル) の計算
! Calculate reference temperature on half levels
!
allocate( z_RefTemp(1:kmax) )
allocate( r_RefTemp(0:kmax) )
z_RefTemp = RefTemp
r_RefTemp(0) = 0.
r_RefTemp(kmax) = 0.
do k = 1, kmax - 1
r_RefTemp(k) = z_TInpCoefA(k+1) * z_RefTemp(k+1) + z_TInpCoefB( k ) * z_RefTemp( k )
enddo
! ルジャンドル陪関数の固有値の設定
! Set eigen value of Associated Legendre Function
!
allocate( w_LaplaEigVal(lmax) )
w_LaplaEigVal(:) = rn(:,1) / RPlanet**2
! 水平拡散係数の設定
! Configure coefficient of horizontal diffusion
!
! 粘性係数の計算 (最大波数の e-folding time が HDEFoldTime となるように)
! Calculate viscosity coefficient
!
HDEFoldTime = DCCalConvertByUnit( HDEFoldTimeValue, HDEFoldTimeUnit, 'sec' ) ! (in)
if ( HDEFoldTimeValue > 0.0_DP ) then
VisCoef = ( (nmax*(nmax+1)) / RPlanet**2 )**(-HDOrder / 2) / HDEFoldTime
else
VisCoef = 0.0_DP
end if
w_HDifCoefH = - VisCoef * ( ( - w_LaplaEigVal )**( HDOrder / 2 ) )
w_HDifCoefM = w_HDifCoefH - VisCoef * ( - ( 2.0_DP / RPlanet**2 )**( HDOrder / 2 ) )
! スポンジ層の減衰係数の設定
! Set damping coefficient in a sponge layer
!
if ( SLEFoldTimeValue <= 0.0_DP ) then
call MessageNotify( 'E', module_name, 'SLEFoldTimeValue must be greater than zero, SLEFoldTimeValue=<%f>.', d = (/ SLEFoldTimeValue /) )
end if
if ( ( SLNumLayer <= 0 ) .or. ( SLNumLayer > kmax ) ) then
call MessageNotify( 'E', module_name, 'SLNumLayer must be greater than zero and less than/equal to kmax, SLNumLayer=<%d>.', i = (/ SLNumLayer /) )
end if
! スポンジ層の減衰係数の計算
! Calculate damping coefficient for momentum in a sponge layer
!
SLEFoldTime = DCCalConvertByUnit( SLEFoldTimeValue, SLEFoldTimeUnit, 'sec' ) ! (in)
if ( FlagSpongeLayer ) then
! Sponge layer is not applied for model layers, k = 1, kmax-SLNumLayers.
!
do k = 1, kmax-SLNumLayer
wz_SpongeLayerCoefM(:,k) = 0.0_DP
end do
do k = kmax-SLNumLayer+1, kmax
wz_SpongeLayerCoefM(:,k) = 1.0_DP / ( SLEFoldTime * ( z_Sigma(k) / z_Sigma(kmax) )**SLOrder )
end do
if ( .not. FlagSpongeLayerforZonalMean ) then
! Sponge layer is not applied for zonal mean component.
! i.e., sponge layer is applied for zonal wave component only.
!
do k = kmax-SLNumLayer+1, kmax
do l = 1, lmax
a_DegOrd = nm_l( l )
if ( a_DegOrd(2) == 0 ) then
wz_SpongeLayerCoefM(l,k) = 0.0_DP
end if
end do
end do
end if
if ( FlagSpongeLayerforHeat ) then
wz_SpongeLayerCoefH = wz_SpongeLayerCoefM
! Sponge layer is not applied for zonal mean component.
! i.e., sponge layer is applied for zonal wave component only.
!
do k = kmax-SLNumLayer+1, kmax
do l = 1, lmax
a_DegOrd = nm_l( l )
if ( a_DegOrd(2) == 0 ) then
wz_SpongeLayerCoefH(l,k) = 0.0_DP
end if
end do
end do
else
wz_SpongeLayerCoefH = 0.0_DP
end if
else
! Sponge layer is not applied.
!
wz_SpongeLayerCoefM = 0.0_DP
wz_SpongeLayerCoefH = 0.0_DP
end if
! 運動量の水平拡散係数とスポンジ層減衰係数の和
! Damping coefficients by horizonatl diffusion and a sponge layer
!
allocate( wz_DisCoefM(lmax, 1:kmax) )
allocate( wz_DisCoefH(lmax, 1:kmax) )
allocate( w_DisCoefQ (lmax) )
do k = 1, kmax
wz_DisCoefM(:,k) = w_HDifCoefM - wz_SpongeLayerCoefM(:,k)
wz_DisCoefH(:,k) = w_HDifCoefH - wz_SpongeLayerCoefH(:,k)
end do
w_DisCoefQ = w_HDifCoefH
! Calculation of divergence damping period
!
DivDampPeriod = DCCalConvertByUnit( DivDampPeriodValue, DivDampPeriodUnit, 'sec' )
if ( DivDampTime0 > TimeN ) then
call MessageNotify( 'E', module_name, 'DivDampTime0, %f, must be less than or equal to TimeN, %f.', d = (/ DivDampTime0, TimeN /) )
end if
! ヒストリデータ出力のためのへの変数登録
! Register of variables for history data output
!
call HistoryAutoAddVariable( 'DUDtDyn', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'dynamical tendency of zonal wind', 'm s-2' )
call HistoryAutoAddVariable( 'DVDtDyn', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'dynamical tendency of meridional wind', 'm s-2' )
call HistoryAutoAddVariable( 'DTempDtDyn', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'dynamical tendency of temperature', 'K s-1' )
do n = 1, ncmax
call HistoryAutoAddVariable( 'D'//trim(a_QMixName(n))//'DtDyn', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'dynamical tendency of water vapor', 's-1' )
end do
call HistoryAutoAddVariable( 'DPsDtDyn', (/ AxNameX, AxNameY, AxNameT /), 'dynamical tendency of surface pressure', 'Ps s-1' )
call HistoryAutoAddVariable( 'OMG', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'vertical velocity in pressure coordinate (omega, DP/Dt)', 'Pa s-1' )
call HistoryAutoAddVariable( 'SigDot', (/ AxNameX, AxNameY, AxNameR, AxNameT /), 'sigma-vertical velocity', '1 s-1' )
call HistoryAutoAddVariable( 'DPiDt', (/ AxNameX, AxNameY, AxNameT /), 'Pi (log Ps) tendency)', 'Pa s-1' )
call HistoryAutoAddVariable( 'Vor', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'vorticity', 's-1' )
call HistoryAutoAddVariable( 'Div', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'divergence', 's-1' )
call HistoryAutoAddVariable( 'Mass', (/ AxNameX, AxNameY, AxNameT /), 'mass', 'kg' )
call HistoryAutoAddVariable( 'KinEngy', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'kinetic energy', 'J m-2' )
call HistoryAutoAddVariable( 'IntEngy', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'internal energy', 'J m-2' )
call HistoryAutoAddVariable( 'PotEngy', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'potential energy', 'J m-2' )
call HistoryAutoAddVariable( 'LatEngy', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'latent energy', 'J m-2' )
call HistoryAutoAddVariable( 'TotEngy', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'total energy', 'J m-2' )
call HistoryAutoAddVariable( 'Enstro', (/ AxNameX, AxNameY, AxNameZ, AxNameT /), 'enstrophy', 'kg m-2 s-2' )
call HistoryAutoAddVariable( 'AdvTestStreamFunc', (/ AxNameX, AxNameY, AxNameR, AxNameT /), 'mass stream function', 'kg s-1' )
! Initialization of modules used in this module
!
! 質量の補正
! Mass fixer
!
call MassFixerInit
if ( FlagAdvTest ) then
!
! Utility module for advection test
!
call AdvTestInit
end if
if ( FlagSLTT ) then
!!$ select case ( IDTTMethod )
!!$ case ( IDTTMethodSp )
!!$ case ( IDTTMethodSL )
! セミラグランジュ法による物質移流
! Semi-Lagrangian method for tracer transport
!
call SLTTInit
!!$ case ( IDTTMethodFV )
!!$ ! 物質移流 (有限体積法)
!!$ ! Tracer Transport (finite volume method)
!!$ !
!!$ call TTFVInit
!!$ end select
end if
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, ' TimeIntegScheme = %c', c1 = trim( TimeIntegScheme ) )
call MessageNotify( 'M', module_name, ' HDEFoldTime = %f [%c]', d = (/ HDEFoldTimeValue /), c1 = trim(HDEFoldTimeUnit) )
call MessageNotify( 'M', module_name, ' HDOrder = %d', i = (/ HDOrder /) )
call MessageNotify( 'M', module_name, ' VisCoef = %f', d = (/ VisCoef /) )
call MessageNotify( 'M', module_name, ' FlagSpongeLayer = %b', l = (/ FlagSpongeLayer /) )
call MessageNotify( 'M', module_name, ' FlagSpongeLayerforZonalMean = %b', l = (/ FlagSpongeLayerforZonalMean /) )
call MessageNotify( 'M', module_name, ' FlagSpongeLayerforHeat = %b', l = (/ FlagSpongeLayerforHeat /) )
call MessageNotify( 'M', module_name, ' SLEFoldTime = %f [%c]', d = (/ SLEFoldTimeValue /), c1 = trim(SLEFoldTimeUnit) )
call MessageNotify( 'M', module_name, ' SLOrder = %d', i = (/ SLOrder /) )
call MessageNotify( 'M', module_name, ' SLNumLayer = %d', i = (/ SLNumLayer /) )
call MessageNotify( 'M', module_name, ' RefTemp = %f', d = (/ RefTemp /) )
call MessageNotify( 'M', module_name, ' FlagDivDamp = %b', l = (/ FlagDivDamp /) )
call MessageNotify( 'M', module_name, ' DivDampPeriod = %f', d = (/ DivDampPeriod /) )
call MessageNotify( 'M', module_name, ' DivDampTime0 = %f', d = (/ DivDampTime0 /) )
call MessageNotify( 'M', module_name, ' FlagTotMassFix = %b', l = (/ FlagTotMassFix /) )
call MessageNotify( 'M', module_name, ' FlagMassHorDifCor = %b', l = (/ FlagMassHorDifCor /) )
call MessageNotify( 'M', module_name, ' FlagSLTT = %b', l = (/ FlagSLTT /) )
!!$ call MessageNotify( 'M', module_name, ' TracerTransportMethod = %c', c1 = trim( TracerTransportMethod ) )
call MessageNotify( 'M', module_name, ' FlagAdvTest = %b', l = (/ FlagAdvTest /) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
dynamics_hspl_vas83_inited = .true.
end subroutine DynamicsHSplVAS83Init
| Variable : | |||
| DelTimeSave : | real(DP), save
|
| Variable : | |||
| DivDampPeriod : | real(DP), save
|
| Variable : | |||
| DivDampTime0 : | real(DP), save
|
| Subroutine : | |||
| wz_DivA(lmax, 1:kmax) : | real(DP), intent(inout)
|
発散の減衰項の付加 (場がつり合っていない計算初期の擾乱の減衰を想定)
Addition of divergence damping term
subroutine DivergenceDamping( wz_DivA )
!
! 発散の減衰項の付加 (場がつり合っていない計算初期の擾乱の減衰を想定)
!
! Addition of divergence damping term
!
! モジュール引用 ; USE statements
!
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN ! ステップ $ t $ の時刻. Time of step $ t $.
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
implicit none
real(DP), intent(inout):: wz_DivA (lmax, 1:kmax)
! $ D (t+\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
! 作業変数
! Work variables
!
real(DP) :: DivDampCoef
! 実行文 ; Executable statement
!
if ( FlagDivDamp ) then
DivDampCoef = 1.0_DP / DelTime * ( DivDampPeriod - ( TimeN - DivDampTime0 ) ) / DivDampPeriod
if ( DivDampCoef < 0.0_DP ) DivDampCoef = 0.0_DP
wz_DivA = wz_DivA / ( 1.0_DP + 2.0_DP * DelTime * DivDampCoef )
end if
end subroutine DivergenceDamping
| Variable : | |||
| FlagAdvTest : | logical , save
|
| Variable : | |||
| FlagMassHorDifCor : | logical , save
|
| Subroutine : | |||
| xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_Phi(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
|
格子点データである温度 $ T $ から, 静水圧の式を用いて 格子点データのジオポテンシャル高度 $ Phi $ を求めます.
subroutine HydroGrid( xyz_Temp, xyz_Phi )
!
! 格子点データである温度 $ T $ から, 静水圧の式を用いて
! 格子点データのジオポテンシャル高度 $ \Phi $ を求めます.
!
! モジュール引用 ; USE statements
!
use constants, only: CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 格子点設定
! Grid points settings
!
use gridset, only: imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $ T $ . 温度. Temperature
real(DP), intent(out):: xyz_Phi (0:imax-1, 1:jmax, 1:kmax)
! $ \Phi $ . ジオポテンシャル高度.
! Getpotential height
! 作業変数
! Work variables
!
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
! 実行文 ; Executable statement
!
xyz_Phi(:,:,1) = CpDry * z_HydroAlpha(1) * xyz_Temp(:,:,1)
do k = 2, kmax
xyz_Phi(:,:,k) = xyz_Phi(:,:,k-1) + CpDry * z_HydroAlpha(k ) * xyz_Temp(:,:,k ) + CpDry * z_HydroBeta (k-1) * xyz_Temp(:,:,k-1)
enddo
end subroutine HydroGrid
| Variable : | |||
| IDTimeIntegScheme : | integer , save
|
| Subroutine : | |||
| xyz_UCosLatN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_VCosLatN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_VorN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DivN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_TempN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_QVapN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_GradLambdaPiN(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_GradMuPiN(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_PiAdv(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_UAdvN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_VAdvN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_TempNonLinearN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_KinEngyN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_TempUAdvN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_TempVAdvN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyr_SigDotN(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out)
| ||
| xy_DPiDtNG(0:imax-1, 1:jmax) : | real(DP), intent(out)
|
非線形項 (非重力波項) を格子点上で計算します.
Non-linear terms (non gravitational terms) are calculated on grid points
subroutine NonLinearOnGrid( xyz_UCosLatN, xyz_VCosLatN, xyz_VorN, xyz_DivN, xyz_TempN, xyz_QVapN, xy_GradLambdaPiN, xy_GradMuPiN, xyz_PiAdv, xyz_UAdvN, xyz_VAdvN, xyz_TempNonLinearN, xyz_KinEngyN, xyz_TempUAdvN, xyz_TempVAdvN, xyr_SigDotN, xy_DPiDtNG )
!
! 非線形項 (非重力波項) を格子点上で計算します.
!
! Non-linear terms (non gravitational terms) are calculated on
! grid points
!
! モジュール引用 ; USE statements
!
! 座標データ設定
! Axes data settings
!
use axesset, only: r_Sigma, z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
! 物理定数設定
! Physical constants settings
!
use constants, only: CpDry, EpsV
! $ \epsilon_v $ .
! 水蒸気分子量比.
! Molecular weight of water vapor
! 文字列操作
! Character handling
!
use dc_string, only: LChar
! 格子点設定
! Grid points settings
!
use gridset, only: imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: ncmax, IndexH2OVap
implicit none
real(DP), intent(in):: xyz_UCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ U (t) = u (t) \cos \varphi $ .
real(DP), intent(in):: xyz_VCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ V (t) = v (t) \cos \varphi $ .
real(DP), intent(in):: xyz_VorN (0:imax-1, 1:jmax, 1:kmax)
! $ \zeta (t) $ . 渦度. Vorticity
real(DP), intent(in):: xyz_DivN (0:imax-1, 1:jmax, 1:kmax)
! $ D (t) $ . 発散. Divergence
real(DP), intent(in):: xyz_TempN (0:imax-1, 1:jmax, 1:kmax)
! $ T (t) $ . 温度. Temperature
real(DP), intent(in):: xyz_QVapN(0:imax-1, 1:jmax, 1:kmax)
! $ q (t) $ . 比湿. Specific humidity
real(DP), intent(in):: xy_GradLambdaPiN (0:imax-1, 1:jmax)
! $ \DP{\pi}{\lambda} $
real(DP), intent(in):: xy_GradMuPiN (0:imax-1, 1:jmax)
! $ (1-\mu^2) \DP{\pi}{\mu} $
real(DP), intent(out):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
! $ \Dvect{v} \cdot \nabla \pi $ .
! $ \pi $ の移流. Advection of $ \pi $
real(DP), intent(out):: xyz_UAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ U_A (t) $ . 東西運動量移流項.
! Eastward advection of momentum
real(DP), intent(out):: xyz_VAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ V_A (t) $ . 南北運動量移流項.
! Northward advection of momentum
real(DP), intent(out):: xyz_TempNonLinearN (0:imax-1, 1:jmax, 1:kmax)
! $ \hat{H} (t) $ . 温度時間変化項.
! Temperature tendency
real(DP), intent(out):: xyz_KinEngyN (0:imax-1, 1:jmax, 1:kmax)
! $ KE (t) $ . 運動エネルギー項.
! Kinetic energy
real(DP), intent(out):: xyz_TempUAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ UT' (t) $ . 温度東西移流項.
! Eastward advection of temperature
real(DP), intent(out):: xyz_TempVAdvN (0:imax-1, 1:jmax, 1:kmax)
! $ VT' (t) $ . 温度南北移流項.
! Northward advection of temperature
real(DP), intent(out):: xyr_SigDotN (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} (t) $ .
! 鉛直流. Vertical flow
real(DP), intent(out):: xy_DPiDtNG (0:imax-1, 1:jmax)
! $ Z (t) $ . 地表面気圧時間変化非重力波項.
! Non-gravity wave component of surface pressure tendency
!-----------------------------------
! 作業変数
! Work variables
!
real(DP):: xyz_PiAdvSum (0:imax-1, 1:jmax, 1:kmax)
! $ \sum_k^K(\Dvect{v}\cdot\nabla\pi)\Delta\sigma $ .
! $ \pi $ 移流の積分値. Integral downward of advection of $ \pi $
real(DP):: xyz_DivSum (0:imax-1, 1:jmax, 1:kmax)
! $ \sum_k^K D\Delta\sigma $ .
! 発散の積分値. Integral downward of divergence
real(DP):: xyr_SigDotNonG (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} $ .
! 鉛直流 (非重力波). Vertical flow (non gravitational)
real(DP):: xyz_TempEdd (0:imax-1, 1:jmax, 1:kmax)
! $ T' = T - \bar{T} $ .
! 温度の擾乱 (整数レベル). Temperature eddy (full level)
real(DP):: xyr_TempEdd (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T}' $ .
! 温度の擾乱 (半整数レベル). Temperature eddy (half level)
real(DP):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
! $ T_v $ .
! 仮温度. Virtual temperature
real(DP):: xyz_VirTempEdd (0:imax-1, 1:jmax, 1:kmax)
! $ {T_v}' = T_v - \bar{T} $ .
! 仮温度の擾乱. Virtual temperature eddy
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! $ \pi $ の移流, $ \pi $ 移流と発散の大気上端から下向きの積分
! Calculate advection of $ \pi $, integral advection of $ \pi $ and divergence
! from the top of the model downward
!
do k = 1, kmax
xyz_PiAdv(:,:,k) = ( xyz_UCosLatN(:,:,k) * xy_GradLambdaPiN + xyz_VCosLatN(:,:,k) * xy_GradMuPiN ) / ( 1.0_DP - xy_SinLat**2 )
enddo
xyz_PiAdvSum(:,:,kmax) = xyz_PiAdv(:,:,kmax) * z_DelSigma(kmax)
do k = kmax-1, 1, -1
xyz_PiAdvSum(:,:,k) = xyz_PiAdvSum(:,:,k+1) + xyz_PiAdv(:,:,k) * z_DelSigma(k)
enddo
xyz_DivSum(:,:,kmax) = xyz_DivN(:,:,kmax) * z_DelSigma(kmax)
do k = kmax-1, 1, -1
xyz_DivSum(:,:,k) = xyz_DivSum(:,:,k+1) + xyz_DivN(:,:,k) * z_DelSigma(k)
enddo
! 地表面気圧時間変化の非重力波項 $ Z $ の計算
! Calculate non-gravity wave component of surface pressure tendency $ Z $
!
xy_DPiDtNG = - xyz_PiAdvSum(:,:,1)
! $ \dot{\sigma} $ の計算
! Calculate $ \dot{\sigma} $
!
do k = 1, kmax-1
xyr_SigDotN(:,:,k) = r_Sigma(k) * ( xyz_PiAdvSum(:,:,1) + xyz_DivSum(:,:,1) ) - ( xyz_PiAdvSum(:,:,k+1) + xyz_DivSum(:,:,k+1) )
xyr_SigDotNonG(:,:,k) = r_Sigma(k) * xyz_PiAdvSum(:,:,1) - xyz_PiAdvSum(:,:,k+1)
enddo
! $ \dot{\sigma} $ の上下境界値
! $ \dot{\sigma} $ on upper and lower boundary
!
xyr_SigDotN (:,:,0 ) = 0.
xyr_SigDotN (:,:,kmax) = 0.
xyr_SigDotNonG(:,:,0 ) = 0.
xyr_SigDotNonG(:,:,kmax) = 0.
! 温度の擾乱 (整数レベル), 仮温度, 仮温度の擾乱の計算
! Calculate temperature eddy (full level), virtual temperature,
! virtual temperature eddy
!
do k = 1, kmax
xyz_VirTemp(:,:,k) = xyz_TempN(:,:,k) * ( 1.0_DP + ((( 1.0_DP / EpsV ) - 1.0_DP ) * xyz_QVapN(:,:,k)) )
xyz_TempEdd(:,:,k) = xyz_TempN(:,:,k) - z_RefTemp(k)
xyz_VirTempEdd(:,:,k) = xyz_VirTemp(:,:,k) - z_RefTemp(k)
enddo
! 温度の擾乱 (半整数レベル) の計算
! Calculate temperature eddy (half level)
!
xyr_TempEdd(:,:,0 ) = 0.0_DP
xyr_TempEdd(:,:,kmax) = 0.0_DP
do k = 1, kmax-1
xyr_TempEdd(:,:,k) = z_TInpCoefA(k+1) * xyz_TempN(:,:,k+1) + z_TInpCoefB(k ) * xyz_TempN(:,:,k ) - r_RefTemp(k)
enddo
! 東西運動量移流項の計算
! Calculate advection of eastward momentum
!
xyz_UAdvN(:,:,1) = ( xyz_VorN(:,:,1) + xy_Cori ) * xyz_VCosLatN(:,:,1) - 1.0_DP / ( 2.0_DP * z_DelSigma(1) ) * xyr_SigDotN(:,:,1) * ( xyz_UCosLatN(:,:,1) - xyz_UCosLatN(:,:,2) ) - CpDry * z_TInpCoefK(1) * xyz_VirTempEdd(:,:,1) * xy_GradLambdaPiN
do k = 2, kmax-1
xyz_UAdvN(:,:,k) = ( xyz_VorN(:,:,k) + xy_Cori ) * xyz_VCosLatN(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * ( xyr_SigDotN(:,:,k-1) * ( xyz_UCosLatN(:,:,k-1) - xyz_UCosLatN(:,:,k) ) + xyr_SigDotN(:,:,k) * ( xyz_UCosLatN(:,:,k) - xyz_UCosLatN(:,:,k+1) ) ) - CpDry * z_TInpCoefK(k) * xyz_VirTempEdd(:,:,k) * xy_GradLambdaPiN
end do
xyz_UAdvN(:,:,kmax) = ( xyz_VorN(:,:,kmax) + xy_Cori ) * xyz_VCosLatN(:,:,kmax) - 1.0_DP / ( 2.0_DP * z_DelSigma(kmax) ) * xyr_SigDotN(:,:,kmax-1) * ( xyz_UCosLatN(:,:,kmax-1) - xyz_UCosLatN(:,:,kmax) ) - CpDry * z_TInpCoefK(kmax) * xyz_VirTempEdd(:,:,kmax) * xy_GradLambdaPiN
! 南北運動量移流項の計算
! Calculate advection of northward momentum
!
xyz_VAdvN(:,:,1) = - ( xyz_VorN(:,:,1) + xy_Cori ) * xyz_UCosLatN(:,:,1) - 1.0_DP / ( 2.0_DP * z_DelSigma(1) ) * xyr_SigDotN(:,:,1) * ( xyz_VCosLatN(:,:,1) - xyz_VCosLatN(:,:,2) ) - CpDry * z_TInpCoefK(1) * xyz_VirTempEdd(:,:,1) * xy_GradMuPiN
do k = 2, kmax-1
xyz_VAdvN(:,:,k) = - ( xyz_VorN(:,:,k) + xy_Cori ) * xyz_UCosLatN(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * ( xyr_SigDotN(:,:,k-1) * ( xyz_VCosLatN(:,:,k-1) - xyz_VCosLatN(:,:,k) ) + xyr_SigDotN(:,:,k) * ( xyz_VCosLatN(:,:,k) - xyz_VCosLatN(:,:,k+1) ) ) - CpDry * z_TInpCoefK(k) * xyz_VirTempEdd(:,:,k) * xy_GradMuPiN
end do
xyz_VAdvN(:,:,kmax) = - ( xyz_VorN(:,:,kmax) + xy_Cori ) * xyz_UCosLatN(:,:,kmax) - 1.0_DP / ( 2.0_DP * z_DelSigma(kmax) ) * xyr_SigDotN(:,:,kmax-1) * ( xyz_VCosLatN(:,:,kmax-1) - xyz_VCosLatN(:,:,kmax) ) - CpDry * z_TInpCoefK(kmax) * xyz_VirTempEdd(:,:,kmax) * xy_GradMuPiN
! 運動エネルギー項 (仮温度補正含む) の計算
! Calculate kinematic energy term
! (including virtual temperature correction)
!
call HydroGrid( xyz_VirTemp - xyz_TempN, xyz_KinEngyN ) ! (out)
do k = 1, kmax
xyz_KinEngyN(:,:,k) = xyz_KinEngyN(:,:,k) + ( xyz_UCosLatN(:,:,k)**2 + xyz_VCosLatN(:,:,k)**2 ) / ( 2.0_DP * ( 1.0_DP - xy_SinLat**2 ) )
end do
! 温度東西移流項, 温度南北移流項の計算
! Calculate eastward and northward advection of temperature
!
do k = 1, kmax
xyz_TempUAdvN(:,:,k) = xyz_UCosLatN(:,:,k) * xyz_TempEdd(:,:,k)
xyz_TempVAdvN(:,:,k) = xyz_VCosLatN(:,:,k) * xyz_TempEdd(:,:,k)
end do
! 温度の時間変化項 $ \hat{H} $ の計算
! Calculate temperature tendency term $ \hat{H} $
!
do k = 1, kmax-1
xyz_TempNonLinearN(:,:,k) = xyz_TempEdd(:,:,k) * xyz_DivN(:,:,k) - 1.0_DP / z_DelSigma(k) * ( xyr_SigDotN(:,:,k-1) * ( xyr_TempEdd(:,:,k-1) - xyz_TempEdd(:,:,k) ) + xyr_SigDotN(:,:,k) * ( xyz_TempEdd(:,:,k) - xyr_TempEdd(:,:,k) ) ) - 1.0_DP / z_DelSigma(k) * ( xyr_SigDotNonG(:,:,k-1) * ( r_RefTemp(k-1) - z_RefTemp(k) ) + xyr_SigDotNonG(:,:,k) * ( z_RefTemp(k) - r_RefTemp(k) ) ) + z_TInpCoefK(k) * xyz_VirTemp(:,:,k) * xyz_PiAdv(:,:,k) - z_HydroAlpha(k) / z_DelSigma(k) * ( xyz_VirTemp(:,:,k) * xyz_PiAdvSum(:,:,k) + xyz_VirTempEdd(:,:,k) * xyz_DivSum(:,:,k) ) - z_HydroBeta(k) / z_DelSigma(k) * ( xyz_VirTemp(:,:,k) * xyz_PiAdvSum(:,:,k+1) + xyz_VirTempEdd(:,:,k) * xyz_DivSum(:,:,k+1) )
enddo
xyz_TempNonLinearN(:,:,kmax) = xyz_TempEdd(:,:,kmax) * xyz_DivN(:,:,kmax) - 1.0_DP / z_DelSigma(kmax) * ( xyr_SigDotN(:,:,kmax-1) * ( xyr_TempEdd(:,:,kmax-1) - xyz_TempEdd(:,:,kmax) ) + xyr_SigDotN(:,:,kmax) * ( xyz_TempEdd(:,:,kmax) - xyr_TempEdd(:,:,kmax) ) ) - 1.0_DP / z_DelSigma(kmax) * ( xyr_SigDotNonG(:,:,kmax-1) * ( r_RefTemp(kmax-1) - z_RefTemp(kmax) ) + xyr_SigDotNonG(:,:,kmax) * ( z_RefTemp(kmax) - r_RefTemp(kmax) ) ) + z_TInpCoefK(kmax) * xyz_VirTemp(:,:,kmax) * xyz_PiAdv(:,:,kmax) - z_HydroAlpha(kmax) / z_DelSigma(kmax) * ( xyz_VirTemp(:,:,kmax) * xyz_PiAdvSum(:,:,kmax) + xyz_VirTempEdd(:,:,kmax) * xyz_DivSum(:,:,kmax) )
end subroutine NonLinearOnGrid
| Subroutine : | |||
| xyz_UCosLatN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_VCosLatN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DivN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_SigDotN(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in )
| ||
| xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xyzf_QMixNonLinearN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xyzf_QMixUAdvN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xyzf_QMixVAdvN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(out)
|
非線形項 (非重力波項) を格子点上で計算します.
Non-linear terms (non gravitational terms) are calculated on grid points
subroutine NonLinearOnGridQMix( xyz_UCosLatN, xyz_VCosLatN, xyz_DivN, xyr_SigDotN, xyzf_QMixN, xyzf_QMixNonLinearN, xyzf_QMixUAdvN, xyzf_QMixVAdvN )
!
! 非線形項 (非重力波項) を格子点上で計算します.
!
! Non-linear terms (non gravitational terms) are calculated on
! grid points
!
! モジュール引用 ; USE statements
!
! 座標データ設定
! Axes data settings
!
use axesset, only: r_Sigma, z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
! 格子点設定
! Grid points settings
!
use gridset, only: imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: ncmax
! 成分の数
! Number of composition
implicit none
real(DP), intent(in):: xyz_UCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ U (t) = u (t) \cos \varphi $ .
real(DP), intent(in):: xyz_VCosLatN (0:imax-1, 1:jmax, 1:kmax)
! $ V (t) = v (t) \cos \varphi $ .
real(DP), intent(in):: xyz_DivN (0:imax-1, 1:jmax, 1:kmax)
! $ D (t) $ . 発散. Divergence
real(DP), intent(in ):: xyr_SigDotN (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} (t) $ .
! 鉛直流. Vertical flow
real(DP), intent(in):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q (t) $ . 比湿. Specific humidity
real(DP), intent(out):: xyzf_QMixNonLinearN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ R (t) $ . 比湿時間変化項.
! Specific humidity tendency
real(DP), intent(out):: xyzf_QMixUAdvN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ Uq (t) $ . 比湿東西移流項.
! Eastward advection of specific humidity
real(DP), intent(out):: xyzf_QMixVAdvN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ Vq (t) $ . 比湿南北移流項.
! Northward advection of specific humidity
!-----------------------------------
! 作業変数
! Work variables
!
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 比湿東西移流項, 比湿南北移流項の計算
! Calculate eastward and northward advection of specific humidity
!
do n = 1, ncmax
do k = 1, kmax
xyzf_QMixUAdvN(:,:,k,n) = xyz_UCosLatN(:,:,k) * xyzf_QMixN(:,:,k,n)
xyzf_QMixVAdvN(:,:,k,n) = xyz_VCosLatN(:,:,k) * xyzf_QMixN(:,:,k,n)
end do
end do
! 比湿時間変化項 $ R $ の計算
! Calculate specific humidity tendency $ R $
!
do n = 1, ncmax
xyzf_QMixNonLinearN(:,:,1,n) = xyzf_QMixN(:,:,1,n) * xyz_DivN(:,:,1) - 1.0_DP / ( 2.0_DP * z_DelSigma(1) ) * xyr_SigDotN(:,:,1) * ( xyzf_QMixN(:,:,1,n) - xyzf_QMixN(:,:,2,n) )
do k = 2, kmax - 1
xyzf_QMixNonLinearN(:,:,k,n) = xyzf_QMixN(:,:,k,n) * xyz_DivN(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * ( xyr_SigDotN(:,:,k-1) * ( xyzf_QMixN(:,:,k-1,n) - xyzf_QMixN(:,:,k ,n) ) + xyr_SigDotN(:,:,k) * ( xyzf_QMixN(:,:,k ,n) - xyzf_QMixN(:,:,k+1,n) ) )
end do
xyzf_QMixNonLinearN(:,:,kmax,n) = xyzf_QMixN(:,:,kmax,n) * xyz_DivN(:,:,kmax) - 1.0_DP / ( 2.0_DP * z_DelSigma(kmax) ) * xyr_SigDotN(:,:,kmax-1) * ( xyzf_QMixN(:,:,kmax-1,n) - xyzf_QMixN(:,:,kmax,n) )
end do
end subroutine NonLinearOnGridQMix
| Subroutine : | |||
| xyz_UB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_VB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_TempB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_PsB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_UN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_VN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_TempN(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_PsN(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_UA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_VA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_TempA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_PsA(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_DUDtPhy(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DVDtPhy(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_DTempDtPhy(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyzf_DQMixDtPhy(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xyr_SigDot(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_DPiDt(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_PiAdv(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyz_OMG(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
|
診断量の出力を行います.
Diagnostic variables are output.
subroutine OutputDiagnosedVariables( xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA, xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy, xyr_SigDot, xy_DPiDt, xyz_PiAdv, xyz_OMG )
!
! 診断量の出力を行います.
!
! Diagnostic variables are output.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: RPlanet, Grav, CpDry, LatentHeat
! $ L $ [J kg-1] .
! 凝結の潜熱.
! Latent heat of condensation
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN ! ステップ $ t $ の時刻. Time of step $ t $.
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: ncmax, a_QMixName, IndexH2OVap
! 座標データ設定
! Axes data settings
!
use axesset, only: z_Sigma, r_Sigma, z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
#ifdef LIB_MPI
#ifdef SJPACK
use wa_mpi_module_sjpack, only: w_xy => w_xv, xy_GradLon_w => xv_GradLon_w, xy_GradLat_w => xv_GradLat_w, wa_DivLambda_xya => wa_DivLambda_xva, wa_DivMu_xya => wa_DivMu_xva, xya_wa => xva_wa
#else
use wa_mpi_module, only: w_xy => w_xv, xy_GradLon_w => xv_GradLon_w, xy_GradLat_w => xv_GradLat_w, wa_DivLambda_xya => wa_DivLambda_xva, wa_DivMu_xya => wa_DivMu_xva, xya_wa => xva_wa
#endif
#elif AXISYMMETRY
use wa_zonal_module, only: w_xy, xy_GradLon_w, xy_GradLat_w, wa_DivLambda_xya, wa_DivMu_xya, xya_wa
#elif SJPACK
use wa_module_sjpack, only: w_xy, xy_GradLon_w, xy_GradLat_w, wa_DivLambda_xya, wa_DivMu_xya, xya_wa
#elif AXISYMMETRY_SJPACK
use wa_zonal_module_sjpack, only: w_xy, xy_GradLon_w, xy_GradLat_w, wa_DivLambda_xya, wa_DivMu_xya, xya_wa
#else
use wa_module, only: w_xy, xy_GradLon_w, xy_GradLat_w, wa_DivLambda_xya, wa_DivMu_xya, xya_wa
#endif
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_UB (0:imax-1, 1:jmax, 1:kmax)
! $ u (t-\Delta t) $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_VB (0:imax-1, 1:jmax, 1:kmax)
! $ v (t-\Delta t) $ . 南北風速. Northward wind
real(DP), intent(in):: xyz_TempB (0:imax-1, 1:jmax, 1:kmax)
! $ T (t-\Delta t) $ . 温度. Temperature
real(DP), intent(in):: xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q (t-\Delta t) $ . 比湿. Specific humidity
real(DP), intent(in):: xy_PsB (0:imax-1, 1:jmax)
! $ p_s (t-\Delta t) $ . 地表面気圧. Surface pressure
real(DP), intent(in):: xyz_UN (0:imax-1, 1:jmax, 1:kmax)
! $ u (t) $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_VN (0:imax-1, 1:jmax, 1:kmax)
! $ v (t) $ . 南北風速. Northward wind
real(DP), intent(in):: xyz_TempN (0:imax-1, 1:jmax, 1:kmax)
! $ T (t) $ . 温度. Temperature
real(DP), intent(in):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q (t) $ . 比湿. Specific humidity
real(DP), intent(in):: xy_PsN (0:imax-1, 1:jmax)
! $ p_s (t) $ . 地表面気圧. Surface pressure
real(DP), intent(in):: xyz_UA (0:imax-1, 1:jmax, 1:kmax)
! $ u (t+\Delta t) $ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_VA (0:imax-1, 1:jmax, 1:kmax)
! $ v (t+\Delta t) $ . 南北風速. Northward wind
real(DP), intent(in):: xyz_TempA (0:imax-1, 1:jmax, 1:kmax)
! $ T (t+\Delta t) $ . 温度. Temperature
real(DP), intent(in):: xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q (t+\Delta t) $ . 比湿. Specific humidity
real(DP), intent(in):: xy_PsA (0:imax-1, 1:jmax)
! $ p_s (t+\Delta t) $ . 地表面気圧. Surface pressure
real(DP), intent(in):: xyz_DUDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{u}{t}\right)^{phy} $ .
! 外力項 (物理過程) による東西風速変化.
! Eastward wind tendency by external force terms (physical processes)
real(DP), intent(in):: xyz_DVDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{v}{t}\right)^{phy} $ .
! 外力項 (物理過程) による南北風速変化.
! Northward wind tendency by external force terms (physical processes)
real(DP), intent(in):: xyz_DTempDtPhy (0:imax-1, 1:jmax, 1:kmax)
! $ \left(\DP{T}{t}\right)^{phy} $ .
! 外力項 (物理過程) による温度変化.
! Temperature tendency by external force terms (physical processes)
real(DP), intent(in):: xyzf_DQMixDtPhy(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \left(\DP{q}{t}\right)^{phy} $ .
! 外力項 (物理過程) による比湿変化.
! Temperature tendency by external force terms (physical processes)
real(DP), intent(in):: xyr_SigDot (0:imax-1, 1:jmax, 0:kmax)
! $ \dot{\sigma} $ .
! 鉛直流. Vertical flow
real(DP), intent(in):: xy_DPiDt (0:imax-1, 1:jmax)
! $ Z $ . 地表面気圧時間変化項.
! Surface pressure tendency
real(DP), intent(in):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
! $ \Dvect{v} \cdot \nabla \pi $ .
! $ \pi $ の移流. Advection of $ \pi $
real(DP), intent(out):: xyz_OMG (0:imax-1, 1:jmax, 1:kmax)
!
! OMEGA, DP/Dt
! 作業変数
! Work variables
!
real(DP):: xyz_DUDtDyn (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_DVDtDyn (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_DTempDtDyn (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyzf_DQMixDtDyn(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP):: xy_DPsDtDyn (0:imax-1, 1:jmax)
real(DP):: xyz_UCosLat (0:imax-1, 1:jmax, 1:kmax)
! $ U = u \cos \varphi $ .
real(DP):: xyz_VCosLat (0:imax-1, 1:jmax, 1:kmax)
! $ V = v \cos \varphi $ .
real(DP):: xyz_Vor (0:imax-1, 1:jmax, 1:kmax)
! $ \zeta $ . 渦度. Vorticity
real(DP):: xyz_Div (0:imax-1, 1:jmax, 1:kmax)
! $ D $ . 発散. Divergence
real(DP):: xyr_PiAdvDivSum(0:imax-1, 1:jmax, 0:kmax)
real(DP):: xy_Mass (0:imax-1, 1:jmax)
! 質量.
! Mass
real(DP):: xyz_KinEngy (0:imax-1, 1:jmax, 1:kmax)
! $ KE $ . 運動エネルギー.
! Kinetic energy
real(DP):: xyz_IntEngy (0:imax-1, 1:jmax, 1:kmax)
! $ IE $ . 内部エネルギー.
! Internal energy
real(DP):: xyz_PotEngy (0:imax-1, 1:jmax, 1:kmax)
! $ PE $ . ポテンシャルエネルギー.
! Potential energy
real(DP):: xyz_LatEngy (0:imax-1, 1:jmax, 1:kmax)
! $ LE $ . 潜熱エネルギー.
! Latent heat energy
real(DP):: xyz_TotEngy (0:imax-1, 1:jmax, 1:kmax)
! $ TE $ . 全エネルギー.
! Total energy
real(DP):: xyz_Enstro (0:imax-1, 1:jmax, 1:kmax)
! エンストロフィー.
! Enstrophy
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n
! 実行文 ; Executable statement
!
! Calculate tendencies
!
xyz_DUDtDyn = ( xyz_UA - xyz_UB ) / ( 2.0_DP * DelTime ) - xyz_DUDtPhy
xyz_DVDtDyn = ( xyz_VA - xyz_VB ) / ( 2.0_DP * DelTime ) - xyz_DVDtPhy
xyz_DTempDtDyn = ( xyz_TempA - xyz_TempB ) / ( 2.0_DP * DelTime ) - xyz_DTempDtPhy
xyzf_DQMixDtDyn = ( xyzf_QMixA - xyzf_QMixB ) / ( 2.0_DP * DelTime ) - xyzf_DQMixDtPhy
xy_DPsDtDyn = ( xy_PsA - xy_PsB ) / ( 2.0_DP * DelTime )
call HistoryAutoPut( TimeN, 'DUDtDyn' , xyz_DUDtDyn )
call HistoryAutoPut( TimeN, 'DVDtDyn' , xyz_DVDtDyn )
call HistoryAutoPut( TimeN, 'DTempDtDyn', xyz_DTempDtDyn )
do n = 1, ncmax
call HistoryAutoPut( TimeN, 'D'//trim(a_QMixName(n))//'DtDyn', xyzf_DQMixDtDyn(:,:,:,n) )
end do
call HistoryAutoPut( TimeN, 'DPsDtDyn' , xy_DPsDtDyn )
! 鉛直流と地表面気圧時間変化項の出力
! Output vertical flow and surface pressure tendency
!
call HistoryAutoPut( TimeN, 'SigDot', xyr_SigDot )
call HistoryAutoPut( TimeN, 'DPiDt', xy_DPiDt )
! 風速から渦度発散の計算
! Calculate vorticity and divergence from wind velocity
!
do k = 1, kmax
xyz_UCosLat(:,:,k) = xyz_UN(:,:,k) * xy_CosLat
xyz_VCosLat(:,:,k) = xyz_VN(:,:,k) * xy_CosLat
end do
xyz_Vor = xya_wa( wa_DivLambda_xya( xyz_VCosLat ) - wa_DivMu_xya( xyz_UCosLat ) ) / RPlanet
xyz_Div = xya_wa( wa_DivLambda_xya( xyz_UCosLat ) + wa_DivMu_xya( xyz_VCosLat ) ) / RPlanet
call HistoryAutoPut( TimeN, 'Vor', xyz_Vor )
call HistoryAutoPut( TimeN, 'Div', xyz_Div )
! Calculation of Omega (DPressDt)
! It should be recognized that the value of xy_Ps here may have been modified
! by mass fixer.
!
! Integration from top of the model to k's layer upper interface
xyr_PiAdvDivSum(:,:,kmax) = 0.0_DP
do k = kmax-1, 0, -1
xyr_PiAdvDivSum(:,:,k) = xyr_PiAdvDivSum(:,:,k+1) + ( xyz_PiAdv(:,:,k+1) + xyz_Div(:,:,k+1) ) * z_DelSigma(k+1)
end do
do k = 1, kmax
xyz_OMG(:,:,k) = xy_PsN * ( z_Sigma(k) * xyz_PiAdv(:,:,k) - xyr_PiAdvDivSum(:,:,k) - ( xyz_PiAdv(:,:,k) + xyz_Div(:,:,k) ) * ( z_Sigma(k) - r_Sigma(k) ) )
end do
call HistoryAutoPut( TimeN, 'OMG', xyz_OMG )
! 質量の計算
! Calculate mass
!
xy_Mass = xy_PsN / Grav
! エネルギー, エンストロフィーの計算
! Calculate energy and enstrophy
!
call HydroGrid( xyz_TempN, xyz_PotEngy )
do k = 1, kmax
xyz_KinEngy(:,:,k) = ( xyz_UN(:,:,k) ** 2 + xyz_VN(:,:,k) ** 2 ) / 2.0_DP * xy_Mass
xyz_IntEngy(:,:,k) = CpDry * xyz_TempN(:,:,k) * xy_Mass
xyz_PotEngy(:,:,k) = xyz_PotEngy(:,:,k) * xy_Mass
xyz_LatEngy(:,:,k) = LatentHeat * xyzf_QMixN(:,:,k,IndexH2OVap) * xy_Mass
end do
xyz_TotEngy = xyz_KinEngy + xyz_IntEngy + xyz_PotEngy + xyz_LatEngy
do k = 1, kmax
xyz_Enstro(:,:,k) = xyz_Vor(:,:,k) ** 2 * xy_Mass
end do
call HistoryAutoPut( TimeN, 'Mass', xy_Mass )
call HistoryAutoPut( TimeN, 'KinEngy', xyz_KinEngy )
call HistoryAutoPut( TimeN, 'IntEngy', xyz_IntEngy )
call HistoryAutoPut( TimeN, 'PotEngy', xyz_PotEngy )
call HistoryAutoPut( TimeN, 'LatEngy', xyz_LatEngy )
call HistoryAutoPut( TimeN, 'TotEngy', xyz_TotEngy )
call HistoryAutoPut( TimeN, 'Enstro', xyz_Enstro )
end subroutine OutputDiagnosedVariables
| Subroutine : |
TimeIntegration で使用する係数の設定を行います. 初回および $ Delta t $ が変更された場合以外は, 前回に設定した値をそのまま用います.
Configure coefficients for "TimeIntegration". Setting values that are set last time are used, except when first time and $ Delta t $ are changed.
subroutine SemiImplMatrix
!
! TimeIntegration で使用する係数の設定を行います.
! 初回および $ \Delta t $ が変更された場合以外は,
! 前回に設定した値をそのまま用います.
!
! Configure coefficients for "TimeIntegration".
! Setting values that are set last time are used,
! except when first time and $ \Delta t $ are changed.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: RPlanet, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
! 座標データ設定
! Axes data settings
!
use axesset, only: r_Sigma, z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
#ifdef LIB_MPI
#ifdef SJPACK
use wa_mpi_module_sjpack, only: l_nm
#else
use wa_mpi_module, only: l_nm
#endif
#elif AXISYMMETRY
use wa_zonal_module, only: l_nm
#elif SJPACK
use wa_module_sjpack, only: l_nm
#elif AXISYMMETRY_SJPACK
use wa_zonal_module_sjpack, only: l_nm
#else
use wa_module, only: l_nm
#endif
! LU 分解法により連立 1 次方程式を解くための関数 (spml 同梱モジュール)
! Functions to solve linear simultaneous equation by LU decomposition
! (a module included in spml)
!
use lumatrix, only: LUDecomp
! デバッグ用ユーティリティ
! Utilities for debug
!
use dc_trace, only: DbgMessage, BeginSub, EndSub, Debug
! 宣言文 ; Declaration statements
!
implicit none
! TimeIntegration 等で使用する係数の設定のための作業変数
! Work variable for coefficients for "TimeIntegration", etc.
!
integer:: k, l, kk ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n
! 実行文 ; Executable statement
!
! $ \Delta t $ [s] のチェック・保存
! Check and save $ \Delta t $ [s]
!
if ( DelTimeSave == DelTime ) return
DelTimeSave = DelTime
call DbgMessage( '%c: %c: (DelTime=%f [sec]) coefficients for "TimeIntegration" is generated. ', c1 = module_name, c2 = 'SemiImplMatrix', d = (/ DelTime /) )
! TimeIntegration で使用する係数の設定
! Configure coefficients for "TimeIntegration"
!
if ( .not. allocated( z_siMtxC ) ) allocate( z_siMtxC (1:kmax) )
if ( .not. allocated( z_siMtxG ) ) allocate( z_siMtxG (1:kmax) )
if ( .not. allocated( zz_siMtxH ) ) allocate( zz_siMtxH (1:kmax, 1:kmax) )
if ( .not. allocated( wz_siMtxDi ) ) allocate( wz_siMtxDi (lmax,1:kmax) )
if ( .not. allocated( zz_siMtxDiH ) ) allocate( zz_siMtxDiH (1:kmax, 1:kmax) )
if ( .not. allocated( wzz_siMtxWDiH ) ) allocate( wzz_siMtxWDiH(lmax,1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxGCt ) ) allocate( zz_siMtxGCt (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxW ) ) allocate( zz_siMtxW (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxQ ) ) allocate( zz_siMtxQ (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxS ) ) allocate( zz_siMtxS (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxR ) ) allocate( zz_siMtxR (1:kmax, 1:kmax) )
z_siMtxC = z_DelSigma
z_siMtxG = CpDry * z_TInpCoefK * z_RefTemp
do k = 1, kmax
do l = 1, kmax
zz_siMtxGCt(k,l) = z_siMtxG(k) * z_siMtxC(l)
end do
end do
zz_siMtxW = 0.
do k = 1, kmax
do l = 1, k
zz_siMtxW(k,l) = CpDry * z_HydroAlpha(l)
enddo
do l = 1, k-1
zz_siMtxW(k,l) = zz_siMtxW(k,l) + CpDry * z_HydroBeta(l)
enddo
enddo
zz_siMtxS = 0.
do k = 1, kmax
do l = 1, kmax
zz_siMtxS(k,l) = r_Sigma(k-1) * z_DelSigma(l)
enddo
do l = k, kmax
zz_siMtxS(k,l) = zz_siMtxS(k,l) - z_DelSigma(l)
enddo
enddo
zz_siMtxQ = 0.
do k = 1, kmax
zz_siMtxQ(k,k) = ( r_RefTemp(k-1) - z_RefTemp(k) ) / z_DelSigma(k)
enddo
do k = 1, kmax-1
zz_siMtxQ(k,k+1) = ( z_RefTemp(k) - r_RefTemp(k) ) / z_DelSigma(k)
enddo
zz_siMtxR = 0.
do k = 1, kmax
do l = k, kmax
zz_siMtxR(k,l) = - z_HydroAlpha(k) / z_DelSigma(k) * z_DelSigma(l) * z_RefTemp(k)
enddo
do l = k + 1, kmax
zz_siMtxR(k,l) = zz_siMtxR(k,l) - z_HydroBeta(k) / z_DelSigma(k) * z_DelSigma(l) * z_RefTemp(k)
enddo
enddo
zz_siMtxH = matmul(zz_siMtxQ, zz_siMtxS) - zz_siMtxR
! Check of coefficients for horizontal diffusion and sponge layer
! A threshold value used below is arbitrary.
!
do n = 1, lmax
do k = 1, kmax
if ( abs( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM(n,k) ) < 1.0e-10_DP ) then
call MessageNotify( 'E', module_name, 'Dissipation coefficient for momentum is inappropriate.' )
end if
if ( abs( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefH(n,k) ) < 1.0e-10_DP ) then
call MessageNotify( 'E', module_name, 'Dissipation coefficient for heat is inappropriate.' )
end if
end do
if ( abs( 1.0_DP - 2.0_DP * DelTime * w_DisCoefQ(n) ) < 1.0e-10_DP ) then
call MessageNotify( 'E', module_name, 'Dissipation coefficient for composition is inappropriate.' )
end if
end do
wz_siMtxDi = 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefH )
do n = 1, lmax
! NOTE: wz_siMtiDi is a diagonal matrix.
!
do k = 1, kmax
do kk = 1, kmax
zz_siMtxDiH(k,kk) = wz_siMtxDi(n,k) * zz_siMtxH(k,kk)
end do
end do
zz_siMtxDiH = matmul( zz_siMtxW, zz_siMtxDiH )
do k = 1, kmax
do kk = 1, kmax
wzz_siMtxWDiH(n,k,kk) = zz_siMtxDiH(k,kk)
end do
end do
end do
if ( .not. allocated(wzz_siMtxLU) ) allocate( wzz_siMtxLU(lmax, 1:kmax, 1:kmax) )
if ( .not. allocated(wz_siMtxPiv) ) allocate( wz_siMtxPiv(lmax, 1:kmax) )
! 行列 $ \underline{M} $ の計算
! Calculate matrix $ \underline{M} $
!
do k = 1, kmax
do kk = 1, kmax
wzz_siMtxLU ( :,k,kk ) = - DelTime**2 * w_LaplaEigVal(:) * ( wzz_siMtxWDiH(:,k,kk) + zz_siMtxGCt(k,kk) )
if ( k == kk ) then
wzz_siMtxLU ( :,k,kk ) = wzz_siMtxLU ( :,k,kk ) + ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM (:,k) )
endif
end do
end do
! LU 行列計算
! LU matrix calculation
!
call LUDecomp( wzz_siMtxLU, wz_siMtxPiv ) ! (out)
end subroutine SemiImplMatrix
| Subroutine : | |||
| w_SurfGeoPot(lmax) : | real(DP), intent(in)
| ||
| wz_DVorDtNG(lmax, 1:kmax) : | real(DP), intent(in)
| ||
| wz_DDivDtNG(lmax, 1:kmax) : | real(DP), intent(in)
| ||
| wz_DTempDtNG(lmax, 1:kmax) : | real(DP), intent(in)
| ||
| w_DPiDtNG(lmax) : | real(DP), intent(in)
| ||
| wz_DivN(lmax, 1:kmax) : | real(DP), intent(in)
| ||
| wz_TempN(lmax, 1:kmax) : | real(DP), intent(in)
| ||
| w_PiN(lmax) : | real(DP), intent(in)
| ||
| wz_VorB(lmax, 1:kmax) : | real(DP), intent(in)
| ||
| wz_DivB(lmax, 1:kmax) : | real(DP), intent(in)
| ||
| wz_TempB(lmax, 1:kmax) : | real(DP), intent(in)
| ||
| w_PiB(lmax) : | real(DP), intent(in)
| ||
| wz_VorA(lmax, 1:kmax) : | real(DP), intent(out)
| ||
| wz_DivA(lmax, 1:kmax) : | real(DP), intent(out)
| ||
| wz_TempA(lmax, 1:kmax) : | real(DP), intent(out)
| ||
| w_PiA(lmax) : | real(DP), intent(out)
|
時間積分を行い, 時刻 $ t $ の物理量の時間変化と $ t-\Delta t$ の物理量から 時刻 $ t+Delta t $ の物理量を計算します.
時間積分法にはリープフロッグスキームを用いています. ただし拡散項には後方差分を用いています. デフォルトでは, $ Delta t $ を大きくとるために, 重力波項に セミインプリシット法を適用しています. NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで, 重力波項をエクスプリシット法によって解くことも可能です.
With time integration, physical values at $ t+Delta t $ is calculated from tendency at $ t $ and physical values at $ t-\Delta t $ .
Leap-frog scheme is used as time integration scheme. And backward scheme is used for diffusion terms. As default setting, semi-implicit scheme is applied to gravitational terms in order to enlarge the value of $ Delta t $ . Explicit scheme can be applied to gravitational terms by changing "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
subroutine TimeIntegration( w_SurfGeoPot, wz_DVorDtNG, wz_DDivDtNG, wz_DTempDtNG, w_DPiDtNG, wz_DivN, wz_TempN, w_PiN, wz_VorB, wz_DivB, wz_TempB, w_PiB, wz_VorA, wz_DivA, wz_TempA, w_PiA )
!
! 時間積分を行い,
! 時刻 $ t $ の物理量の時間変化と $ t-\Delta t$ の物理量から
! 時刻 $ t+\Delta t $ の物理量を計算します.
!
! 時間積分法にはリープフロッグスキームを用いています.
! ただし拡散項には後方差分を用いています.
! デフォルトでは, $ \Delta t $ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! With time integration, physical values at $ t+\Delta t $ is calculated
! from tendency at $ t $ and physical values at $ t-\Delta t $ .
!
! Leap-frog scheme is used as time integration scheme.
! And backward scheme is used for diffusion terms.
! As default setting, semi-implicit scheme is applied to gravitational terms
! in order to enlarge the value of $ \Delta t $ .
! Explicit scheme can be applied to gravitational terms by changing
! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: RPlanet, CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! 座標データ設定
! Axes data settings
!
use axesset, only: z_DelSigma ! $ \Delta \sigma $ (整数).
! $ \Delta \sigma $ (Full)
! LU 分解法により連立 1 次方程式を解くための関数 (spml 同梱モジュール)
! Functions to solve linear simultaneous equation by LU decomposition
! (a module included in spml)
!
use lumatrix, only: LUSolve
! 文字列操作
! Character handling
!
use dc_string, only: LChar
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: ncmax, IndexH2OVap
implicit none
real(DP), intent(in):: w_SurfGeoPot (lmax)
! $ \Phi_s $ . 地表ジオポテンシャル.
! Surface geo-potential
real(DP), intent(in):: wz_DVorDtNG (lmax, 1:kmax)
! $ \left( \DD{\zeta}{t} (t) \right)^{NG}$ . 渦度変化の非重力波成分 (スペクトル).
! Non-gravity wave component of vorticity tendency (spectral)
real(DP), intent(in):: wz_DDivDtNG (lmax, 1:kmax)
! $ \DD{D}{t} (t) $ . 発散変化の非重力波成分 (スペクトル).
! Divergence tendency (spectral)
real(DP), intent(in):: wz_DTempDtNG(lmax, 1:kmax)
! $ \left( \DD{T}{t} (t) \right)^{NG}$ . 温度変化の非重力波成分 (スペクトル).
! Temperature tendency (spectral)
real(DP), intent(in):: w_DPiDtNG(lmax)
! $ \left( \DD{p_s}{t} (t) \right) $ . 地表面気圧変化の非重力波項 (スペクトル).
! Non-gravity wave component of surface pressure tendency (spectral)
real(DP), intent(in):: wz_DivN (lmax, 1:kmax)
! $ D (t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP), intent(in):: wz_TempN (lmax, 1:kmax)
! $ T (t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP), intent(in):: w_PiN (lmax)
! $ \pi = \ln p_s (t) $ . 地表面気圧 (スペクトル).
real(DP), intent(in):: wz_VorB (lmax, 1:kmax)
! $ \zeta (t-\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP), intent(in):: wz_DivB (lmax, 1:kmax)
! $ D (t-\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP), intent(in):: wz_TempB (lmax, 1:kmax)
! $ T (t-\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP), intent(in):: w_PiB (lmax)
! $ \pi = \ln p_s (t-\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP), intent(out):: wz_VorA (lmax, 1:kmax)
! $ \zeta (t+\Delta t) $ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP), intent(out):: wz_DivA (lmax, 1:kmax)
! $ D (t+\Delta t) $ . 発散 (スペクトル).
! Divergence (spectral)
real(DP), intent(out):: wz_TempA (lmax, 1:kmax)
! $ T (t+\Delta t) $ . 温度 (スペクトル).
! Temperature (spectral)
real(DP), intent(out):: w_PiA (lmax)
! $ \pi = \ln p_s (t+\Delta t) $ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
! 作業変数
! Work variables
!
real(DP):: wz_HDiv (lmax, 1:kmax)
real(DP):: w_CtDiv (lmax)
real(DP):: wz_WT (lmax, 1:kmax)
real(DP):: wz_siTemp (lmax, 1:kmax)
! 温度 (セミインプリシット法のための作業変数).
! Temperature (work variable for semi-implicit scheme)
real(DP):: w_siPi (lmax)
! $ \pi $ (セミインプリシット法のための作業変数).
! $ \pi $ (work variable for semi-implicit scheme)
real(DP):: wz_siPhi (lmax, 1:kmax)
! $ \Phi = \underline{W} ( 1 - 2 \Delta t \underline{D_H} ) \overline{ \Dvect{T} }^{t}$ .
! (セミインプリシット法のための作業変数).
! (Work variable for semi-implicit scheme)
real(DP):: wz_siVectF (lmax, 1:kmax)
! $ \Dvect{f} $ .
! 発散項に関するセミインプリシット方程式の右辺.
! Right-hand side of a semi-implicit equation of a divergence term.
real(DP):: wz_siDivAvrTime (lmax, 1:kmax)
! $ \overline{\Dvect{D}}^{t} $ .
! 時間平均の $ \Dvect{D} $ (セミインプリシット法のための作業変数).
! Time average $ \Dvect{D} $ (a work variable for semi-implicit scheme)
integer:: k, kk ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 時間積分. 拡散は後方差分
! Time integration. Backward difference is applied to diffusion
!
! 渦度 ; Vorticity
!
wz_VorA = ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM ) ) * ( wz_VorB + 2.0_DP * DelTime * wz_DVorDtNG )
! 発散 ; Divergence
!
select case ( IDTimeIntegScheme )
case ( IDTimeIntegSchemeSemiImplicit )
! セミインプリシット法で用いる行列の計算に使われる項の計算
! Calculate terms used in making a matrix for semi-implicit method
!
wz_siTemp = ( 1.0_DP - DelTime * wz_DisCoefH ) * wz_TempB + DelTime * wz_DTempDtNG
! NOTE:
! The matrix wz_siMtxDi = ( 1 - 2 \Delta t D_H )^{-1} is diagonal matrix.
!
wz_siPhi (:,1) = CpDry * z_HydroAlpha(1) * wz_siMtxDi(:,1) * wz_siTemp(:,1)
do k = 2, kmax
wz_siPhi (:,k) = wz_siPhi(:,k-1) + CpDry * z_HydroAlpha(k ) * wz_siMtxDi(:,k ) * wz_siTemp(:,k ) + CpDry * z_HydroBeta (k-1) * wz_siMtxDi(:,k-1) * wz_siTemp(:,k-1)
end do
w_siPi = w_PiB + DelTime * w_DPiDtNG
! 発散方程式の右辺の計算
! Calculate right side of divergence equation
!
do k = 1, kmax
wz_siVectF(:,k) = ( 1.0_DP - DelTime * wz_DisCoefM (:,k) ) * wz_DivB(:,k) + DelTime * wz_DDivDtNG(:,k) - DelTime * w_LaplaEigVal(:) * ( w_SurfGeoPot + wz_siPhi(:,k) + z_siMtxG(k) * w_siPi )
end do
! 時間平均の $ \Dvect{D} $ を LU 行列で解く
! Solve time average $ \Dvect{D} $ with LU matrix
!
wz_siDivAvrTime = LUSolve( wzz_siMtxLU, wz_siMtxPiv, wz_siVectF )
wz_DivA = 2. * wz_siDivAvrTime - wz_DivB
case ( IDTimeIntegSchemeExplicit )
wz_WT = 0.0_DP
do k = 1, kmax
do kk = 1, kmax
wz_WT(:,k) = wz_WT(:,k) + zz_siMtxW(k,kk) * wz_TempN(:,kk)
end do
end do
do k = 1, kmax
wz_DivA(:,k) = ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM(:,k) ) ) * ( wz_DivB(:,k) + 2.0_DP * DelTime * ( wz_DDivDtNG(:,k) - w_LaplaEigVal(:) * ( w_SurfGeoPot + wz_WT(:,k) + z_siMtxG(k) * w_PiN ) ) )
end do
end select
! 温度 ; Temperature
!
select case ( IDTimeIntegScheme )
case ( IDTimeIntegSchemeSemiImplicit )
wz_HDiv = 0.
do k = 1, kmax
do kk = 1, kmax
wz_HDiv(:,k) = wz_HDiv(:,k) + zz_siMtxH(k,kk) * wz_siDivAvrTime(:,kk)
end do
end do
case ( IDTimeIntegSchemeExplicit )
wz_HDiv = 0.
do k = 1, kmax
do kk = 1, kmax
wz_HDiv(:,k) = wz_HDiv(:,k) + zz_siMtxH(k,kk) * wz_DivN(:,kk)
end do
end do
end select
wz_TempA = ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefH ) ) * ( wz_TempB + 2.0_DP * DelTime * ( wz_DTempDtNG - wz_HDiv ) )
! 地表面気圧 ; Surface pressure
!
select case ( IDTimeIntegScheme )
case ( IDTimeIntegSchemeSemiImplicit )
w_CtDiv = 0.0_DP
do k = 1, kmax
w_CtDiv = w_CtDiv + z_siMtxC(k) * wz_siDivAvrTime(:,k)
end do
case ( IDTimeIntegSchemeExplicit )
w_CtDiv = 0.0_DP
do k = 1, kmax
w_CtDiv = w_CtDiv + z_siMtxC(k) * wz_DivN(:,k)
end do
end select
w_PiA = w_PiB + 2.0_DP * DelTime * ( w_DPiDtNG - w_CtDiv )
end subroutine TimeIntegration
| Subroutine : | |||
| wzf_DQMixDtN(lmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| wzf_QMixB(lmax, 1:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| wzf_QMixA(lmax, 1:kmax, 1:ncmax) : | real(DP), intent(out)
|
時間積分を行い, 時刻 $ t $ の物理量の時間変化と $ t-\Delta t$ の物理量から 時刻 $ t+Delta t $ の物理量を計算します.
時間積分法にはリープフロッグスキームを用いています. ただし拡散項には後方差分を用いています. デフォルトでは, $ Delta t $ を大きくとるために, 重力波項に セミインプリシット法を適用しています. NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで, 重力波項をエクスプリシット法によって解くことも可能です.
With time integration, physical values at $ t+Delta t $ is calculated from tendency at $ t $ and physical values at $ t-\Delta t $ .
Leap-frog scheme is used as time integration scheme. And backward scheme is used for diffusion terms. As default setting, semi-implicit scheme is applied to gravitational terms in order to enlarge the value of $ Delta t $ . Explicit scheme can be applied to gravitational terms by changing "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
subroutine TimeIntegrationQMix( wzf_DQMixDtN, wzf_QMixB, wzf_QMixA )
!
! 時間積分を行い,
! 時刻 $ t $ の物理量の時間変化と $ t-\Delta t$ の物理量から
! 時刻 $ t+\Delta t $ の物理量を計算します.
!
! 時間積分法にはリープフロッグスキームを用いています.
! ただし拡散項には後方差分を用いています.
! デフォルトでは, $ \Delta t $ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! With time integration, physical values at $ t+\Delta t $ is calculated
! from tendency at $ t $ and physical values at $ t-\Delta t $ .
!
! Leap-frog scheme is used as time integration scheme.
! And backward scheme is used for diffusion terms.
! As default setting, semi-implicit scheme is applied to gravitational terms
! in order to enlarge the value of $ \Delta t $ .
! Explicit scheme can be applied to gravitational terms by changing
! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
!
! モジュール引用 ; USE statements
!
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! 格子点設定
! Grid points settings
!
use gridset, only: lmax, imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: ncmax
! 成分の数
! Number of composition
implicit none
real(DP), intent(in):: wzf_DQMixDtN(lmax, 1:kmax, 1:ncmax)
! $ \DD{q}{t} (t) $ . 比湿変化 (スペクトル).
! Specific humidity tendency (spectral)
real(DP), intent(in):: wzf_QMixB(lmax, 1:kmax, 1:ncmax)
! $ q (t-\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP), intent(out):: wzf_QMixA(lmax, 1:kmax, 1:ncmax)
! $ q (t+\Delta t) $ . 比湿 (スペクトル).
! Specific humidity (spectral)
! 作業変数
! Work variables
!
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 時間積分. 拡散は後方差分
! Time integration. Backward difference is applied to diffusion
!
! 比湿 ; Specific humidity
!
do n = 1, ncmax
do k = 1, kmax
wzf_QMixA(:,k,n) = ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * w_DisCoefQ ) ) * ( wzf_QMixB(:,k,n) + 2.0_DP * DelTime * wzf_DQMixDtN(:,k,n) )
end do
end do
end subroutine TimeIntegrationQMix
| Variable : | |||
| dynamics_hspl_vas83_inited = .false. : | logical, save
|
| Constant : | |||
| module_name = ‘dynamics_hspl_vas83‘ : | character(*), parameter
|
| Variable : | |||
| r_RefTemp(:) : | real(DP), allocatable
|
| Constant : | |||
| version = ’$Name: $’ // ’$Id: dynamics_hspl_vas83.F90,v 1.83 2014/05/07 09:39:17 murashin Exp $’ : | character(*), parameter
|
| Variable : | |||
| w_DisCoefQ(:) : | real(DP), save, allocatable
|
| Variable : | |||
| w_LaplaEigVal(:) : | real(DP), save, allocatable
|
| Variable : | |||
| wz_DisCoefH(:,:) : | real(DP), save, allocatable
|
| Variable : | |||
| wz_DisCoefM(:,:) : | real(DP), save, allocatable
|
| Variable : | |||
| wz_siMtxDi(:,:) : | real(DP), save, allocatable
|
| Variable : | |||
| wz_siMtxPiv(:,:) : | integer , save, allocatable
|
| Variable : | |||
| wzz_siMtxLU(:,:,:) : | real(DP), save, allocatable
|
| Variable : | |||
| xy_Cori(:,:) : | real(DP), allocatable
|
| Variable : | |||
| z_HydroAlpha(:) : | real(DP), allocatable
|
| Variable : | |||
| z_HydroBeta(:) : | real(DP), allocatable
|
| Variable : | |||
| z_RefTemp(:) : | real(DP), allocatable
|
| Variable : | |||
| z_TInpCoefA(:) : | real(DP), allocatable
|
| Variable : | |||
| z_TInpCoefB(:) : | real(DP), allocatable
|
| Variable : | |||
| z_TInpCoefK(:) : | real(DP), allocatable
|
| Variable : | |||
| zz_siMtxGCt(:,:) : | real(DP), save, allocatable
|
| Variable : | |||||||||
| zz_siMtxR(:,:) : | real(DP), save, allocatable
|
| Variable : | |||
| zz_siMtxW(:,:) : | real(DP), save, allocatable
|