Class dynamics_hspl_vas83 In: dynamics/dynamics_hspl_vas83.F90

Dynamical Core (Spectral method, Arakawa and Suarez (1983))

Note that Japanese and English are described in parallel.

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 for extension of $Delta t$ . Explicit scheme can be applied to gravitational terms by changing "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".

Procedures List

 Dynamics : 力学計算 ———— : ———— Dynamics : Calculate dynamics

References

• Bourke, W. P., 1988: Spectral methods in global climate and weather prediction models. Physically Based Modelling and Simulation of Climates and Climatic Change. Part I., M.E. Schlesinger (ed.), Kluwer Academic Publishers, Dordrecht, 169--220.
• Arakawa, A., Suarez, M. J., 1983: Vertical differencing of the primitive equations in sigma coordinates. Mon. Wea. Rev., 111, 34--35.

Included Modules

gridset dc_types constants timeset wa_mpi_module wa_module gtool_historyauto dc_string axesset lumatrix w_module namelist_util gtool_history dc_iounit dc_date_types dc_date dc_message dc_present dc_trace

Public Instance methods

Subroutine :
xyz_UB(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $u (t-\Delta t)$ . 東西風速. Eastward wind
xyz_VB(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $v (t-\Delta t)$ . 南北風速. Northward wind
xyz_TempB(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $T (t-\Delta t)$ . 温度. Temperature
xyz_QVapB(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $q (t-\Delta t)$ . 比湿. Specific humidity
xy_PsB(0:imax-1, 1:jmax) :real(DP), intent(in)
 : $p_s (t-\Delta t)$ . 地表面気圧. Surface pressure
xyz_UN(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $u (t)$ . 東西風速. Eastward wind
xyz_VN(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $v (t)$ . 南北風速. Northward wind
xyz_TempN(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $T (t)$ . 温度. Temperature
xyz_QVapN(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $q (t)$ . 比湿. Specific humidity
xy_PsN(0:imax-1, 1:jmax) :real(DP), intent(in)
 : $p_s (t)$ . 地表面気圧. Surface pressure
xyz_DUDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $DP{u}{t}$ . 東西風速変化. Eastward wind tendency
xyz_DVDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $DP{v}{t}$ . 南北風速変化. Northward wind tendency
xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $DP{T}{t}$ . 温度変化. Temperature tendency
xyz_DQVapDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $DP{q}{t}$ . 比湿変化. Temperature tendency
xy_SurfHeight(0:imax-1, 1:jmax) :real(DP), intent(in)
 : $z_s$ . 地表面高度. Surface height.
xyz_UA(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
 : $u (t+Delta t)$ . 東西風速. Eastward wind
xyz_VA(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
 : $v (t+Delta t)$ . 南北風速. Northward wind
xyz_TempA(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
 : $T (t+Delta t)$ . 温度. Temperature
xyz_QVapA(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
 : $q (t+Delta t)$ . 比湿. Specific humidity
xy_PsA(0:imax-1, 1:jmax) :real(DP), intent(out)
 : $p_s (t+Delta t)$ . 地表面気圧. Surface pressure

Calculate dynamical processes. Eastward wind (xyz_UB, xyz_UN), northward wind (xyz_VB, xyz_VN), temperature (xyz_TempB, xyz_TempN), specific humidity (xyz_QVapB, xyz_QVapN), surface pressure (xyz_PsB, xyz_PsN) at $t-\Delta t$ and $t$, and surface height (xy_SurfHeight) are input, and eastward wind (xyz_UA), northward wind (xyz_VA), temperature (xyz_TempA), specific humidity (xyz_QVapA), surface pressure (xyz_PsA) are returned.

In order to add tendencies of vorticity, divergence, temperature and specific humidity by another physical process to tendencies of this dynamical process and to calculate the values at next time step, give these tendencies to "xyz_DUDt", "xyz_DVDt", "xyz_DTempDt", "xyz_DQVapDt"

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".

[Source]

  subroutine Dynamics( xyz_UB,   xyz_VB,   xyz_TempB,   xyz_QVapB,   xy_PsB, xyz_UN,   xyz_VN,   xyz_TempN,   xyz_QVapN,   xy_PsN, xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyz_DQVapDt, xy_SurfHeight, xyz_UA,   xyz_VA,   xyz_TempA,   xyz_QVapA,   xy_PsA )
!
! 力学過程の演算を行い, 与えられた $t-\Delta t$ および $t$ の
! 東西風速 (xyz_UB, xyz_UN), 南北風速 (xyz_VB, xyz_VN),
! 温度 (xyz_TempB, xyz_TempN), 比湿 (xyz_QVapB, xyz_QVapN),
! 地表面気圧 (xyz_PsB, xyz_PsN), および地表面高度 (xy_SurfHeight) から,
! $t+\Delta t$ の
! 東西風速 (xyz_UA), 南北風速 (xyz_VA), 温度 (xyz_TempA),
! 比湿 (xyz_QVapA), 地表面気圧 (xyz_PsA) を返します.
!
! 別の物理プロセスによる渦度, 発散, 温度, 比湿の変化を,
! 力学過程の変化に足して次のステップを計算する場合には,
! それらの変化を xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyz_DQVapDt
! に与えてください.
!
! 時間積分法にはリープフロッグスキームを用いています.
! デフォルトでは, $\Delta t$ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! Calculate dynamical processes.
! Eastward wind (xyz_UB, xyz_UN),
! northward wind (xyz_VB, xyz_VN),
! temperature (xyz_TempB, xyz_TempN),
! specific humidity (xyz_QVapB, xyz_QVapN),
! surface pressure (xyz_PsB, xyz_PsN) at $t-\Delta t$ and $t$,
! and surface height (xy_SurfHeight) are input,
! and
! eastward wind (xyz_UA), northward wind (xyz_VA),
! temperature (xyz_TempA),
! specific humidity (xyz_QVapA), surface pressure (xyz_PsA)
! are returned.
!
! In order to add tendencies of vorticity, divergence,
! temperature and specific humidity by another physical process to
! tendencies of this dynamical process
! and to calculate the values at next time step,
! give these tendencies to
! "xyz_DUDt", "xyz_DVDt", "xyz_DTempDt", "xyz_DQVapDt"
!
! 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: imax, jmax, kmax    ! 鉛直層数.
! Number of vertical level

! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

#ifdef LIB_MPI
! MPI 版 SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library MPI version, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
#else
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
#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.

! 宣言文 ; 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):: xyz_QVapB (0:imax-1, 1:jmax, 1:kmax)
! $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):: xyz_QVapN (0:imax-1, 1:jmax, 1:kmax)
! $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_DUDt    (0:imax-1, 1:jmax, 1:kmax)
! $\DP{u}{t}$ . 東西風速変化.
! Eastward wind tendency
real(DP), intent(in):: xyz_DVDt    (0:imax-1, 1:jmax, 1:kmax)
! $\DP{v}{t}$ . 南北風速変化.
! Northward wind tendency
real(DP), intent(in):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $\DP{T}{t}$ . 温度変化.
! Temperature tendency
real(DP), intent(in):: xyz_DQVapDt (0:imax-1, 1:jmax, 1:kmax)
! $\DP{q}{t}$ . 比湿変化.
! Temperature tendency

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):: xyz_QVapA (0:imax-1, 1:jmax, 1:kmax)
! $q (t+\Delta t)$ .   比湿. Specific humidity
real(DP), intent(out):: xy_PsA    (0:imax-1, 1:jmax)
! $p_s (t+\Delta t)$ . 地表面気圧. Surface pressure

! 作業変数
! Work variables
!
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):: w_PiN ((nmax+1)**2)
! $\pi$ スペクトル
! $\frac{1}{\cos \varphi} \DP{\pi}{\lambda}$
! $\DP{\pi}{\varphi}$

! $u_A (t)$ . 東西運動量移流項.
! $v_A (t)$ . 南北運動量移流項.
real(DP):: xyz_DTempDtN (0:imax-1, 1:jmax, 1:kmax)
! $H (t)$ . 温度時間変化項.
! Temperature tendency
real(DP):: xyz_DQVapDtN (0:imax-1, 1:jmax, 1:kmax)
! $R (t)$ . 比湿時間変化項.
! Specific humidity tendency
real(DP):: xyz_KinEngyN (0:imax-1, 1:jmax, 1:kmax)
! $KE (t)$ . 運動エネルギー項.
! Kinetic energy
! $uT (t)$ . 温度東西移流項.
! $vT (t)$ . 温度南北移流項.
! $\dot{\sigma}$ .
! 鉛直流. Vertical flow
real(DP):: xy_DPiDtN (0:imax-1, 1:jmax)
! $Z$ . 地表面気圧時間変化項.
! Surface pressure tendency
! $uq (t)$ . 比湿東西移流項.
! Eastward advection of specific humidity
! $vq (t)$ . 比湿南北移流項.
! Northward advection of specific humidity
real(DP):: wz_DVorDtN ((nmax+1)**2, 1:kmax)
! $\DD{\zeta}{t} (t)$ . 渦度変化 (スペクトル).
! Vorticity tendency (spectral)
real(DP):: wz_DDivDtN ((nmax+1)**2, 1:kmax)
! $\DD{D}{t} (t)$ . 発散変化 (スペクトル).
! Divergence tendency (spectral)
real(DP):: wz_DTempDtN ((nmax+1)**2, 1:kmax)
! $\DD{T}{t} (t)$ . 温度変化 (スペクトル).
! Temperature tendency (spectral)
real(DP):: wz_DQVapDtN ((nmax+1)**2, 1:kmax)
! $\DD{q}{t} (t)$ . 比湿変化 (スペクトル).
! Specific humidity tendency (spectral)
real(DP):: w_DPiDtN ((nmax+1)**2)
! $\DD{p_s}{t} (t)$ . 地表面気圧変化 (スペクトル).
! Surface pressure tendency (spectral)
real(DP):: wz_VorB ((nmax+1)**2, 1:kmax)
! $\zeta (t-\Delta t)$ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP):: wz_DivB ((nmax+1)**2, 1:kmax)
! $D (t-\Delta t)$ . 発散 (スペクトル).
! Divergence (spectral)
real(DP):: wz_TempB ((nmax+1)**2, 1:kmax)
! $T (t-\Delta t)$ . 温度 (スペクトル).
! Temperature (spectral)
real(DP):: w_PiB ((nmax+1)**2)
! $\pi = \ln p_s (t-\Delta t)$ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP):: wz_QVapB ((nmax+1)**2, 1:kmax)
! $q (t-\Delta t)$ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP):: wz_VorA ((nmax+1)**2, 1:kmax)
! $\zeta (t+\Delta t)$ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP):: wz_DivA ((nmax+1)**2, 1:kmax)
! $D (t+\Delta t)$ . 発散 (スペクトル).
! Divergence (spectral)
real(DP):: wz_TempA ((nmax+1)**2, 1:kmax)
! $T (t+\Delta t)$ . 温度 (スペクトル).
! Temperature (spectral)
real(DP):: w_PiA ((nmax+1)**2)
! $\pi = \ln p_s (t+\Delta t)$ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP):: wz_QVapA ((nmax+1)**2, 1:kmax)
! $q (t+\Delta t)$ . 比湿 (スペクトル).
! Specific humidity (spectral)

real(DP):: wz_Psi ((nmax+1)**2, 1:kmax)
! $\psi$ . 流線関数. Streamline function
real(DP):: wz_Chi ((nmax+1)**2, 1:kmax)
! $\chi$ . ポテンシャル. Potential

real(DP):: wz_VorDiffA ((nmax+1)**2, 1:kmax)
! $\mathscr{D}(\zeta)$ .
! 運動量水平拡散による渦度変化 (スペクトル).
! Vorticity tendency by
! horizontal momentum diffusion (spectral)
real(DP):: wz_DivDiffA ((nmax+1)**2, 1:kmax)
! $\mathscr{D}(D)$ .
! 運動量水平拡散による発散変化 (スペクトル).
! Divergence tendency by
! horizontal momentum diffusion (spectral)
real(DP):: wz_PsiDiff ((nmax+1)**2, 1:kmax)
! 運動量水平拡散による流線関数 $\psi$ 変化
! Streamline function tendency by
! horizontal momentum diffusion
real(DP):: wz_ChiDiff ((nmax+1)**2, 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 ((nmax+1)**2)
! $\Phi_s$ . 地表ジオポテンシャル.
! Surface geo-potential

! 外部のプロセスの時間変化項のための作業変数
! Work variables for tendency of external processes
!
real(DP):: wz_DVorDtWork ((nmax+1)**2, 1:kmax)
! $\DD{\zeta}{t} (t)$ . 渦度変化 (スペクトル).
! Vorticity tendency (spectral)
real(DP):: wz_DDivDtWork ((nmax+1)**2, 1:kmax)
! $\DD{D}{t} (t)$ . 発散変化 (スペクトル).
! Divergence tendency (spectral)
real(DP):: wz_DTempDtWork ((nmax+1)**2, 1:kmax)
! $\DD{T}{t} (t)$ . 温度変化 (スペクトル).
! Temperature tendency (spectral)
real(DP):: wz_DQVapDtWork ((nmax+1)**2, 1:kmax)
! $\DD{q}{t} (t)$ . 比湿変化 (スペクトル).
! Specific humidity tendency (spectral)

! エクスプリシット法のための作業変数
! Work variables for explicit scheme
!
real(DP):: xyz_exWTGPi (0:imax-1, 1:jmax, 1:kmax)
! $\underline{W} \Dvect{T} + \Dvect{G} \pi$ .
real(DP):: xyz_exWT (0:imax-1, 1:jmax, 1:kmax)
! $\underline{W} \Dvect{T}$ .
real(DP):: xyz_exGPi (0:imax-1, 1:jmax, 1:kmax)
! $\Dvect{G} \pi$ .

real(DP):: xyz_exHDiv (0:imax-1, 1:jmax, 1:kmax)
! $\underline{h} \Dvect{D}$ .

integer:: k, kk           ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction

! 実行文 ; Executable statement
!

! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )

! 初期化
! Initialization
!
if ( .not. dynamics_hspl_vas83_inited ) call DynamicsInit

! 地表面気圧の空間変化項の計算
! Calculate spatial surface pressure tendency
!
xy_PiN = log( xy_PsN )
w_PiN = w_xy( xy_PiN )

! 風速から渦度発散の計算
! Calculate vorticity and divergence from wind velocity
!
xyz_VorN = xya_wa( wa_Div_xya_xya( xyz_VN , - xyz_UN ) / RPlanet )
xyz_DivN = xya_wa( wa_Div_xya_xya( xyz_UN ,   xyz_VN ) / RPlanet )

! 格子点上での非線形力学項の計算
! Calculate non-linear dynamical terms on grid points
!

! スペクトル時間変化項の計算
! Calculate spectral tendency terms
!

! 渦度の時間変化 (スペクトル) の計算
! Calculate vorticity tendency (spectral)
!

! 発散の時間変化 (スペクトル) の計算
! Calculate divergence tendency (spectral)
!
wz_DDivDtN = wa_Div_xya_xya( xyz_UAdvN, xyz_VAdvN ) / RPlanet - wa_Lapla_wa( wa_xya(xyz_KinEngyN) ) / RPlanet**2

! 温度の時間変化 (スペクトル) の計算
! Calculate temperature tendency (spectral)
!
wz_DTempDtN = - wa_Div_xya_xya( xyz_TempUAdvN, xyz_TempVAdvN ) / RPlanet + wa_xya( xyz_DTempDtN )

! 地表面気圧の時間変化 (スペクトル) の計算
! Calculate surface pressure tendency (spectral)
!
w_DPiDtN = w_xy( xy_DPiDtN )

! 比湿の時間変化 (スペクトル) の計算
! Calculate specific humidity tendency (spectral)
!
wz_DQVapDtN = - wa_Div_xya_xya( xyz_QVapUAdvN, xyz_QVapVAdvN ) / RPlanet + wa_xya( xyz_DQVapDtN )

! 地表ジオポテンシャルの計算
! Calculate surface geo-potential
!
w_SurfGeoPot = w_xy( Grav * xy_SurfHeight )

! エクスプリシット法を用いる際の計算
! Calculate for explicit scheme
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('explicit')

! $\underline{W} \Dvect{T} + \Dvect{G} \pi$ の格子点値の計算
! Calculate $\underline{W} \Dvect{T} + \Dvect{G} \pi$ on grid
!

! $\underline{W} \Dvect{T}$ の計算
! Calculate $\underline{W} \Dvect{T}$
!
call HydroGrid( xyz_TempN, xyz_exWT )   ! (out)

! $\Dvect{G} \pi$ の計算
! Calculate $\Dvect{G} \pi$
!
do k = 1, kmax
xyz_exGPi(:,:,k) = CpDry * z_TempInpolKappa(k) * z_RefTemp(k) * xy_PiN
enddo

! $\underline{W} \Dvect{T} + \Dvect{G} \pi$ の計算
! Calculate $\underline{W} \Dvect{T} + \Dvect{G} \pi$
!
xyz_exWTGPi = xyz_exWT + xyz_exGPi

! $\underline{h} \Dvect{D}$ の格子点値の計算
! Calculate $\underline{h} \Dvect{D}$ on grid
!
xyz_exHDiv = 0.0_DP

do k = 1, kmax
do kk = 1, kmax
xyz_exHDiv (:,:,k) = xyz_exHDiv (:,:,k) + zz_siMtxH (k,kk) * xyz_DivN (:,:,kk)
enddo
enddo

! 発散, 温度 (スペクトル) の時間変化項の修正
! Modify divergence and temperature tencency (spectral)
!

! 発散の時間変化 (スペクトル) の計算
! Calculate divergence tendency (spectral)
!
wz_DDivDtN = wz_DDivDtN - wa_Lapla_wa( wa_xya( xyz_exWTGPi ) ) / RPlanet**2

do k = 1, kmax
wz_DDivDtN(:,k) = wz_DDivDtN(:,k) - w_Lapla_w( w_SurfGeoPot ) / RPlanet**2
end do

! 温度の時間変化 (スペクトル) の計算
! Calculate temperature tendency (spectral)
!
wz_DTempDtN = wz_DTempDtN - wa_xya( xyz_exHDiv )

end select

! 格子点値をスペクトル値へ ( $t-\Delta t$ )
! Exchange grid values to spectral values ( $t-\Delta t$ )
!
wz_VorB  = wa_Div_xya_xya( xyz_VB, - xyz_UB ) / RPlanet
wz_DivB  = wa_Div_xya_xya( xyz_UB,   xyz_VB ) / RPlanet
wz_TempB = wa_xya( xyz_TempB )
wz_QVapB = wa_xya( xyz_QVapB )
w_PiB    = w_xy( log( xy_PsB ) )

! 外部プロセスによる時間変化格子データをスペクトルデータへ変換
! Convert tendency data on grid by external processes into spectral data
!
wz_DVorDtWork  = wa_Div_xya_xya( xyz_DVDt, - xyz_DUDt ) / RPlanet
wz_DDivDtWork  = wa_Div_xya_xya( xyz_DUDt,   xyz_DVDt ) / RPlanet
wz_DTempDtWork = wa_xya( xyz_DTempDt )
wz_DQVapDtWork = wa_xya( xyz_DQVapDt )

! TimeIntegration で使用する係数の設定
! Configure coefficients for "TimeIntegration"
!
call ImplCoef

! 時間積分
! Time integration
!
call TimeIntegration( wz_DVorDtN + wz_DVorDtWork, wz_DDivDtN + wz_DDivDtWork, wz_DTempDtN + wz_DTempDtWork, wz_DQVapDtN + wz_DQVapDtWork, w_DPiDtN, wz_VorB, wz_DivB, wz_TempB, wz_QVapB, w_PiB, w_SurfGeoPot, wz_VorA, wz_DivA, wz_TempA, wz_QVapA, w_PiA )  ! (out)

! スペクトル値を格子点値へ ( $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_UA = (   xya_GradLon_wa( wz_Chi ) - xya_GradLat_wa( wz_Psi )    ) / RPlanet

xyz_VA = (   xya_GradLon_wa( wz_Psi ) + xya_GradLat_wa( wz_Chi )    ) / RPlanet

xyz_TempA = xya_wa( wz_TempA )
xyz_QVapA = xya_wa( wz_QVapA )
xy_PsA = exp( xy_w( w_PiA ) )

! 拡散による補正
! Correction by diffusion
!

! 運動量水平拡散による渦度発散の時間変化
! Vorticity and divergence tendency by horizontal diffusion of momentum
!
wz_VorDiffA = - wz_VorA * wz_DiffVorDiv
wz_DivDiffA = - wz_DivA * wz_DiffVorDiv

! 運動量水平拡散による摩擦熱補正
! Frictional thermal correction by horizontal momentum diffusion
!
wz_PsiDiff = wa_LaplaInv_wa( wz_VorDiffA ) * RPlanet**2
wz_ChiDiff = wa_LaplaInv_wa( wz_DivDiffA ) * RPlanet**2

xyz_UDiff = (   xya_GradLon_wa( wz_ChiDiff ) - xya_GradLat_wa( wz_PsiDiff )    ) / RPlanet

xyz_VDiff = (   xya_GradLon_wa( wz_PsiDiff ) + xya_GradLat_wa( wz_ChiDiff )    ) / RPlanet

xyz_TempA = xyz_TempA - (   xyz_UA * xyz_UDiff + xyz_VA * xyz_VDiff   ) / CpDry * 2.0_DP * DelTime

! ヒストリデータ出力
! History data output
!
call HistoryAutoPut( TimeN, 'DPiDt',    xy_DPiDtN )
call DiagOutput( xyz_UA, xyz_VA, xyz_TempA, xyz_QVapA, xy_PsA ) ! (in)

! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )

end subroutine Dynamics

dynamics_hspl_vas83_inited
Variable :
dynamics_hspl_vas83_inited = .false. :logical, save, public
 : 初期設定フラグ. Initialization flag

Private Instance methods

DelTimeSave
Variable :
DelTimeSave :real(DP), save
 : 前回の $Delta t$ [s]. 陰解法のための係数の再生成の必要性の チェックに使用する. $Delta t$ [s] at previous step This is used to check necessity of regeneration of coefficients for implicit method.
Subroutine :
xyz_U(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $u$ . 東西風速. Eastward wind
xyz_V(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $v$ . 南北風速. Northward wind
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $T$ . 温度. Temperature
xyz_QVap(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $q$ . 比湿. Specific humidity
xy_Ps(0:imax-1, 1:jmax) :real(DP), intent(in)
 : $p_s$ . 地表面気圧. Surface pressure

Diagnostic variables are output.

[Source]

  subroutine DiagOutput( xyz_U, xyz_V, xyz_Temp, xyz_QVap, xy_Ps )
!
! 診断量の出力を行います.
!
! 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: TimeN  ! ステップ $t$ の時刻. Time of step $t$.

#ifdef LIB_MPI
! MPI 版 SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library MPI version, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_mpi_module, only: wa_Div_xya_xya => wa_Div_xva_xva, xya_wa => xva_wa
#else
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_module, only: wa_Div_xya_xya, xya_wa
#endif

! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut

! 宣言文 ; Declaration statements
!
implicit none
real(DP), intent(in):: xyz_U    (0:imax-1, 1:jmax, 1:kmax)
! $u$ .   東西風速. Eastward wind
real(DP), intent(in):: xyz_V    (0:imax-1, 1:jmax, 1:kmax)
! $v$ .   南北風速. Northward wind
real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $T$ .   温度. Temperature
real(DP), intent(in):: xyz_QVap (0:imax-1, 1:jmax, 1:kmax)
! $q$ .   比湿. Specific humidity
real(DP), intent(in):: xy_Ps    (0:imax-1, 1:jmax)
! $p_s$ . 地表面気圧. Surface pressure

! 作業変数
! Work variables
!
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):: 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

! 実行文 ; Executable statement
!

! 風速から渦度発散の計算
! Calculate vorticity and divergence from wind velocity
!
xyz_Vor = xya_wa( wa_Div_xya_xya( xyz_V , - xyz_U ) / RPlanet )
xyz_Div = xya_wa( wa_Div_xya_xya( xyz_U ,   xyz_V ) / RPlanet )

call HistoryAutoPut( TimeN, 'Vor', xyz_Vor )
call HistoryAutoPut( TimeN, 'Div', xyz_Div )

! 質量の計算
! Calculate mass
!
xy_Mass = xy_Ps / Grav

! エネルギー, エンストロフィーの計算
! Calculate energy and enstrophy
!
call HydroGrid( xyz_Temp, xyz_PotEngy ) ! (out)

do k = 1, kmax
xyz_KinEngy(:,:,k) = ( xyz_U(:,:,k) ** 2 + xyz_V(:,:,k) ** 2 ) / 2.0_DP * xy_Mass

xyz_IntEngy(:,:,k) = CpDry * xyz_Temp(:,:,k) * xy_Mass

xyz_PotEngy(:,:,k) = xyz_PotEngy(:,:,k) * xy_Mass

xyz_LatEngy(:,:,k) = LatentHeat * xyz_QVap(:,:,k) * 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 DiagOutput

 Subroutine :

This procedure input/output NAMELIST#dynamics_hspl_vas83_nml .

[Source]

  subroutine DynamicsInit

! モジュール引用 ; 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

! 座標データ設定
! Axes data settings
!
use axesset, only: z_Sigma, r_Sigma, z_DelSigma            ! $\Delta \sigma$ (整数).
! $\Delta \sigma$ (Full)

#ifdef LIB_MPI
! MPI 版 SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library MPI version, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_mpi_module, only: xy_Lat => xv_Lat, w_xy => w_xv
#else

! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
use wa_module, only: xy_Lat, w_xy
#endif
use w_module, only: rn

! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg

! ヒストリデータ出力
! History data output
!

! 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_date_types, only: DC_DIFFTIME
! 日時の差を表現するデータ型.
! Data type for difference about date and time
use dc_date, only: DCDiffTimeCreate, EvalSec

! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify

! 組み込み関数 PRESENT の拡張版関数
! Extended functions of intrinsic function "PRESENT"
!
use dc_present, only: present_and_not_empty

! 文字列操作
! Character handling
!
use dc_string, only: LChar

! 宣言文 ; Declaration statements
!
implicit none

! 基準温度の設定のための作業変数
! Work variable for reference temperature
!
real(DP):: RefTemp
! $\overline{T}$ . 基準温度.
! Reference temperature

! 水平拡散係数の設定のための作業変数
! Work variable for coefficient of horizontal diffusion
!
integer:: VisOrder        ! 超粘性の次数.  Order of hyper-viscosity
real(DP):: VisCoef        ! 超粘性係数. Hyper-viscosity coefficient
real(DP):: EFoldTimeValue ! 最大波数に対する e-folding time.
! 負の値を与えると, 水平拡散係数をゼロにします.
!
! E-folding time for maximum wavenumber.
! If negative value is given,
! coefficients of horizontal diffusion become zero.
character(TOKEN):: EFoldTimeUnit
! 最大波数に対する e-folding time の単位.
! Unit of e-folding time for maximum wavenumber
real(DP):: EFoldTime      ! 最大波数に対する e-folding time [単位: 秒].
! E-folding time for maximum wavenumber [Unit: sec]
type(DC_DIFFTIME):: dcdiff_efold

! 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:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
! Unit number for NAMELIST file open
integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.

! NAMELIST 変数群
! NAMELIST group name
!
namelist /dynamics_hspl_vas83_nml/ TimeIntegScheme, VisOrder, EFoldTimeValue, EFoldTimeUnit, RefTemp
!
! デフォルト値については初期化手続 "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
call InitCheck

! デフォルト値の設定
! Default values settings
!
TimeIntegScheme = 'Semi-implicit'

VisOrder = 8
EFoldTimeValue = 8640.0_DP
EFoldTimeUnit  = 'sec'

RefTemp = 300.0_DP

! 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')
case ('explicit')
case default
call MessageNotify( 'E', module_name, '"TimeIntegScheme" must be "Semi-Implicit" or "Explicit".' )
end select

! NonLinearOnGrid 等で使用する係数の設定
! Configure coefficients for "NonLinearOnGrid", etc.
!

! コリオリパラメータの計算計算
! Calculate Coriolis parameter
!
allocate( xy_Cori (0:imax-1, 1:jmax) )
xy_Cori = 2.0_DP * Omega * sin( xy_Lat )

! 静水圧の式の係数 $\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_TempInpolA     (1:kmax) )
allocate( z_TempInpolB     (1:kmax) )
allocate( z_TempInpolKappa (1:kmax) )

do k = 1, kmax
z_TempInpolKappa(k) = (   r_Sigma(k-1) * z_HydroAlpha(k) + r_Sigma(k  ) * z_HydroBeta(k)    ) / z_DelSigma(k)
enddo

z_TempInpolA = 0.0_DP
do k = 2, kmax
z_TempInpolA(k) = z_HydroAlpha(k) / ( 1.0_DP - (z_Sigma(k) / z_Sigma(k-1)) ** Kappa )
end do

z_TempInpolB = 0.0_DP
do k = 1, kmax - 1
z_TempInpolB(k) = z_HydroBeta(k) / ( (z_Sigma(k) / z_Sigma(k+1) ) ** Kappa - 1.0_DP )
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.0_DP
r_RefTemp(kmax) = 0.0_DP

do k = 1, kmax - 1
r_RefTemp(k) = z_TempInpolA(k+1) * z_RefTemp(k+1) + z_TempInpolB( k ) * z_RefTemp( k )
enddo

! 水平拡散係数の設定
! Configure coefficient of horizontal diffusion
!
allocate( wz_rn         ((nmax+1)**2, 1:kmax) )
allocate( wz_DiffVorDiv ((nmax+1)**2, 1:kmax) )
allocate( wz_DiffTherm  ((nmax+1)**2, 1:kmax) )

do k = 1, kmax
wz_rn(:,k) = rn(:,1)
enddo

! 粘性係数の計算 (最大波数の e-folding time が EFoldTime となるように)
! Calculate viscosity coefficient
!
call DCDiffTimeCreate( dcdiff_efold, EFoldTimeValue, EFoldTimeUnit )       ! (in)
EFoldTime = EvalSec( dcdiff_efold )

if ( EFoldTimeValue > 0.0_DP ) then
VisCoef = ( (nmax*(nmax+1)) / RPlanet**2 )**(-VisOrder / 2) / EFoldTime
else
VisCoef = 0.0_DP
end if

wz_DiffTherm = VisCoef * ( (-wz_rn / RPlanet**2)**(VisOrder / 2) )

wz_DiffVorDiv = wz_DiffTherm + VisCoef * ( - (2.0_DP / RPlanet**2)**(VisOrder / 2))

! TimeIntegration で使用する係数の設定
! Configure coefficients for "TimeIntegration"
!
call ImplCoef

! ヒストリデータ出力のためのへの変数登録
! Register of variables for history data output
!
call HistoryAutoAddVariable( 'SigmaDot', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'sigma-vertical velocity', '1 s-1' )
call HistoryAutoAddVariable( 'DPiDt', (/ 'lon ', 'lat ', 'time' /), 'Pi (log Ps) tendency)', 'Pa s-1' )

call HistoryAutoAddVariable( 'Vor', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'vorticity', 's-1' )
call HistoryAutoAddVariable( 'Div', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'divergence', 's-1' )

call HistoryAutoAddVariable( 'Mass', (/ 'lon ', 'lat ', 'time' /), 'mass', 'kg' )
call HistoryAutoAddVariable( 'KinEngy', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'kinetic energy', 'J' )
call HistoryAutoAddVariable( 'IntEngy', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'internal energy', 'J' )
call HistoryAutoAddVariable( 'PotEngy', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'potential energy', 'J' )
call HistoryAutoAddVariable( 'LatEngy', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'latent energy', 'J' )
call HistoryAutoAddVariable( 'TotEngy', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'total energy', 'J' )
call HistoryAutoAddVariable( 'Enstro', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'enstrophy', 'kg' )

! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, '  TimeIntegScheme  = %c', c1 = trim( TimeIntegScheme ) )
call MessageNotify( 'M', module_name, '  EFoldTime = %f [%c]', d = (/ EFoldTimeValue /), c1 = trim(EFoldTimeUnit) )
call MessageNotify( 'M', module_name, '  VisOrder  = %d', i = (/ VisOrder /) )
call MessageNotify( 'M', module_name, '  VisCoef   = %f', d = (/ VisCoef /) )
call MessageNotify( 'M', module_name, '  RefTemp   = %f', d = (/ RefTemp /) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

dynamics_hspl_vas83_inited = .true.
end subroutine DynamicsInit

Subroutine :
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $T$ . 温度. Temperature
xyz_Phi(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
 : $Phi$ . ジオポテンシャル高度. Getpotential height

[Source]

  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

! 宣言文 ; 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

 Subroutine :

TimeIntegration で使用する係数の設定.

Configure coefficients for "TimeIntegration"

[Source]

  subroutine ImplCoef
!
! TimeIntegration で使用する係数の設定.
!
! Configure coefficients for "TimeIntegration"
!

! モジュール引用 ; USE statements
!

! 物理定数設定
! Physical constants settings
!
use constants, only: RPlanet, CpDry
! $C_p$ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure

! 座標データ設定
! 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)
!
use wa_module, only: l_nm

! LU 分解法により連立 1 次方程式を解くための関数
! Functions to solve linear simultaneous equation by LU decomposition
!
use lumatrix, only: LUDecomp

! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify

! デバッグ用ユーティリティ
! Utilities for debug
!
use dc_trace, only: DbgMessage, BeginSub, EndSub, Debug

! 宣言文 ; Declaration statements
!
implicit none

! TimeIntegration 等で使用する係数の設定のための作業変数
! Work variable for coefficients for "TimeIntegration", etc.
!
real(DP), allocatable:: zz_siMtxW (:,:)
! $W$ .
! 発散の式での線形重力波項の効果による温度補正の係数.
! Coefficient for correction of temperature
! by effort of linear gravitational terms

real(DP), allocatable:: zz_siMtxQ (:,:)
! $Q = \DD{T}{\sigma}$
real(DP), allocatable:: zz_siMtxS (:,:)
! $S = \DD{\sigma}{t}$
real(DP), allocatable:: zz_siMtxQS (:,:)
! $QS$ .
! この変数は $\sigma$ 移流の効果に相当.
! This variable  corresponds to effort of $\sigma$ advection
real(DP), allocatable:: zz_siMtxR (:,:)
! $R$ .
! 発散と掛け合わせることで気圧変化の効果となる.
! Product R and divergence become effort of
! surface pressure tendency.
! $RD = \kappa T ! (\DD{\pi}{t} + \Dinv{\sigma}\DD{\sigma}{t})$ .

integer, allocatable:: nmo (:,:,:)
! スペクトルの添字順番.
! Spectral subscript expression
real(DP), allocatable:: zz_siMtxM (:,:)
! 行列 $\underline{M}$.
! Matrix $\underline{M}$
integer, allocatable:: z_siMtxPivWork(:)
! 行列のピボット作成の作業変数.
! Work variable for pivot
real(DP):: flapla         ! $\nabla_{\sigma}^{2}$

integer:: k, l, kk        ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction

integer:: mmax            ! 最大東西波数.
! Maximum truncated eastward wavenumber
integer:: lmax            ! 最大南北波数.
! Maximum truncated northward wavenumber
integer:: n, m, nm, mxnm
! 波数方向に回る DO ループ用作業変数
! Work variables for DO loop in wavenumber direction

logical:: lmax_err        ! 最大南北波数に関するエラーフラグ.
! Error flag for maximum truncated northward wavenumber

! 実行文 ; Executable statement
!

! $\Delta t$ [s] のチェック・保存
! Check and save $\Delta t$ [s]
!
if ( .not. dynamics_hspl_vas83_inited ) then
DelTimeSave = DelTime
else
if ( DelTimeSave == DelTime ) then
return
else
DelTimeSave = DelTime
end if
end if

call DbgMessage( '%c:%c: (DelTime=%f) coefficients for "TimeIntegration" is generated. ', c1 = module_name, c2 = 'ImplCoef', d = (/ DelTime /) )

! TimeIntegration で使用する係数の設定
! Configure coefficients for "TimeIntegration"
!
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( zz_siMtxWH ) )  allocate( zz_siMtxWH  (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_siMtxQS ) )  allocate( zz_siMtxQS (1:kmax, 1:kmax) )
if ( .not. allocated( zz_siMtxR  ) )  allocate( zz_siMtxR  (1:kmax, 1:kmax) )

z_siMtxG = CpDry * z_TempInpolKappa * z_RefTemp

do k = 1, kmax
do l = 1, kmax
zz_siMtxGCt(k,l) = z_siMtxG(k) * z_DelSigma(l)
end do
end do

zz_siMtxW = 0.0_DP
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.0_DP
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.0_DP
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_siMtxQS = 0.0_DP
zz_siMtxQS = matmul(zz_siMtxQ, zz_siMtxS)

zz_siMtxR = 0.0_DP
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 = 0.0_DP
zz_siMtxH = zz_siMtxQS - zz_siMtxR

zz_siMtxWH = 0.0_DP
zz_siMtxWH = matmul(zz_siMtxW, zz_siMtxH)

if ( .not. allocated(nmo           ) )  allocate( nmo           (1:2, 0:nmax, 0:nmax) )
if ( .not. allocated(zz_siMtxM     ) )  allocate( zz_siMtxM     (1:kmax, 1:kmax) )
if ( .not. allocated(z_siMtxPivWork) )  allocate( z_siMtxPivWork(1:kmax) )
if ( .not. allocated(wzz_siMtxLU   ) )  allocate( wzz_siMtxLU   ((nmax+1)**2, 1:kmax, 1:kmax) )
if ( .not. allocated(wz_siMtxPiv   ) )  allocate( wz_siMtxPiv   ((nmax+1)**2, 1:kmax) )

mmax = nmax
lmax = nmax
mxnm = 0

! スペクトル添字順序の取り出し
! Fetch spectral subscript expression
!
nmo = 0
do l = 0, lmax
do m = 0, min(mmax, nmax-l)
nmo(1,m,l) = l_nm(m+l,  m)
nmo(2,m,l) = l_nm(m+l, -m)
enddo
enddo

Loop_n: do n = 0, nmax
flapla = - real(n) * real(n+1) / RPlanet**2

! スペクトル添字順序の取り出し
! Fetch spectral subscript expression
!
lmax_err = .true.
do m = 0, min(n,mmax)
if ( n-m <= lmax ) then
nm = nmo(1,m,n-m)
lmax_err = .false.
exit
endif
end do
if ( lmax_err ) then
call MessageNotify( 'E', module_name, 'n-m=<%d> must be less than or equal to lmax=<%d>.', i = (/ n-m, lmax /) )
end if

! 行列 $\underline{M}$ の計算
! Calculate matrix $\underline{M}$
!
do k = 1, kmax
do kk = 1, kmax
zz_siMtxM ( k,kk ) = - DelTime**2 * flapla * (   zz_siMtxWH( k,kk ) + zz_siMtxGCt( k,kk ) * ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm(nm,1) )  )
if ( k == kk ) then
zz_siMtxM ( k,kk ) = zz_siMtxM ( k,kk ) +   ( 1.0_DP + DelTime * 2.0_DP * wz_DiffVorDiv(nm,1) ) * ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm(nm,1) )
endif
end do
end do

! LU 行列計算
! LU matrix calculation
!
call LUDecomp( zz_siMtxM, z_siMtxPivWork ) ! (out)

! ダミー値の代入. (位置 kmax はまだ未定義なため).
! Dummy value is subtituted. (Because position kmax is undefined yet).
!
z_siMtxPivWork(kmax) = 0

!!$write(*,*) 'n= ', n !!$
!!$do k = 1, kmax !!$        do kk = 1, kmax
!!$write(*,*) 'zz_siMtxM(', k+1, ',', kk+1, ')= ', zz_siMtxM(k,kk) !!$        end do
!!$end do !!$      do k = 1, kmax
!!$write(*,*) 'z_siMtxPivWork(', k, ')= ', z_siMtxPivWork(k) !!$      end do

! 配列の詰め替え
! Repack matrices
!
do m = 0, mmax
l = n - m
if ( ( l >= 0 ) .and. ( l <= lmax ) ) then
do k = 1, kmax
do kk = 1, kmax
wzz_siMtxLU ( nmo(1,m,l),k,kk ) = zz_siMtxM ( k,kk )
wzz_siMtxLU ( nmo(2,m,l),k,kk ) = zz_siMtxM ( k,kk )
end do
wz_siMtxPiv ( nmo(1,m,l),k ) = z_siMtxPivWork ( k )
wz_siMtxPiv ( nmo(2,m,l),k ) = z_siMtxPivWork ( k )
mxnm = max( mxnm, nmo(1,m,l) )
mxnm = max( mxnm, nmo(2,m,l) )
end do
endif
end do

end do Loop_n

do nm = mxnm+1, (nmax+1)**2
do k = 1, kmax
do kk = 1, kmax
wzz_siMtxLU ( nm,k,kk ) = zz_siMtxM ( k,kk )
end do
wz_siMtxPiv ( nm,k ) = z_siMtxPivWork ( k )
end do
end do

end subroutine ImplCoef

 Subroutine :

Check initialization of dependency modules

[Source]

  subroutine InitCheck
!
! 依存モジュールの初期化チェック
!
! Check initialization of dependency modules

! モジュール引用 ; USE statements
!

! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_util_inited

! 格子点設定
! Grid points settings
!
use gridset, only: gridset_inited

! 物理定数設定
! Physical constants settings
!
use constants, only: constants_inited

! 座標データ設定
! Axes data settings
!
use axesset, only: axesset_inited

! 時刻管理
! Time control
!
use timeset, only: timeset_inited

! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify

! 実行文 ; Executable statement
!

if ( .not. namelist_util_inited ) call MessageNotify( 'E', module_name, '"namelist_util" module is not initialized.' )

if ( .not. gridset_inited ) call MessageNotify( 'E', module_name, '"gridset" module is not initialized.' )

if ( .not. constants_inited ) call MessageNotify( 'E', module_name, '"constants" module is not initialized.' )

if ( .not. axesset_inited ) call MessageNotify( 'E', module_name, '"axesset" module is not initialized.' )

if ( .not. timeset_inited ) call MessageNotify( 'E', module_name, '"timeset" module is not initialized.' )

end subroutine InitCheck

Subroutine :
xyz_U(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $u$ . 東西風速. Eastward wind
xyz_V(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $v$ . 南北風速. Northward wind
xyz_Vor(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $zeta$ . 渦度. Vorticity
xyz_Div(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $D$ . 発散. Divergence
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $T$ . 温度. Temperature
xyz_QVap(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
 : $q$ . 比湿. Specific humidity
 : $frac{1}{cos varphi} DP{pi}{lambda}$
 : $DP{pi}{varphi}$
 : $u_A$ . 東西運動量移流項. Eastward advection of momentum
 : $v_A$ . 南北運動量移流項. Northward advection of momentum
xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
 : $hat{H}$ . 温度時間変化項. Temperature tendency
xyz_DQVapDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
 : $R$ . 比湿時間変化項. Specific humidity tendency
xyz_KinEngy(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
 : $KE$ . 運動エネルギー項. Kinetic energy
 : $uT$ . 温度東西移流項. Eastward advection of temperature
 : $vT$ . 温度南北移流項. Northward advection of temperature
 : $dot{sigma}$ . 鉛直流. Vertical flow
xy_DPiDt(0:imax-1, 1:jmax) :real(DP), intent(out)
 : $Z$ . 地表面気圧時間変化項. Surface pressure tendency
 : $uq$ . 比湿東西移流項. Eastward advection of specific humidity
 : $vq$ . 比湿南北移流項. Northward advection of specific humidity

Non-linear terms (non gravitational terms) are calculated on grid points

[Source]

  subroutine NonLinearOnGrid( xyz_U,        xyz_V, xyz_Vor,      xyz_Div, xyz_Temp,     xyz_QVap, xy_GradLonPi, xy_GradLatPi, xyz_UAdv,     xyz_VAdv, xyz_DTempDt,  xyz_DQVapDt, xyz_KinEngy, xyz_TempUAdv, xyz_TempVAdv, xyr_SigmaDot, xy_DPiDt, xyz_QVapUAdv, xyz_QVapVAdv )
!
! 非線形項 (非重力波項) を格子点上で計算します.
!
! 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: RPlanet, CpDry, EpsV
! $\epsilon_v$ .
! 水蒸気分子量比.
! Molecular weight of water vapor

! 文字列操作
! Character handling
!
use dc_string, only: LChar

implicit none
real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
! $u$ . 東西風速. Eastward wind
real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
! $v$ . 南北風速. Northward wind
real(DP), intent(in):: xyz_Vor (0:imax-1, 1:jmax, 1:kmax)
! $\zeta$ . 渦度. Vorticity
real(DP), intent(in):: xyz_Div (0:imax-1, 1:jmax, 1:kmax)
! $D$ . 発散. Divergence
real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $T$ . 温度. Temperature
real(DP), intent(in):: xyz_QVap (0:imax-1, 1:jmax, 1:kmax)
! $q$ . 比湿. Specific humidity
! $\frac{1}{\cos \varphi} \DP{\pi}{\lambda}$
! $\DP{\pi}{\varphi}$
real(DP), intent(out):: xyz_UAdv (0:imax-1, 1:jmax, 1:kmax)
! $u_A$ . 東西運動量移流項.
real(DP), intent(out):: xyz_VAdv (0:imax-1, 1:jmax, 1:kmax)
! $v_A$ . 南北運動量移流項.
real(DP), intent(out):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $\hat{H}$ . 温度時間変化項.
! Temperature tendency
real(DP), intent(out):: xyz_DQVapDt (0:imax-1, 1:jmax, 1:kmax)
! $R$ . 比湿時間変化項.
! Specific humidity tendency
real(DP), intent(out):: xyz_KinEngy (0:imax-1, 1:jmax, 1:kmax)
! $KE$ . 運動エネルギー項.
! Kinetic energy
real(DP), intent(out):: xyz_TempUAdv (0:imax-1, 1:jmax, 1:kmax)
! $uT$ . 温度東西移流項.
real(DP), intent(out):: xyz_TempVAdv (0:imax-1, 1:jmax, 1:kmax)
! $vT$ . 温度南北移流項.
real(DP), intent(out):: xyr_SigmaDot (0:imax-1, 1:jmax, 0:kmax)
! $\dot{\sigma}$ .
! 鉛直流. Vertical flow
real(DP), intent(out):: xy_DPiDt (0:imax-1, 1:jmax)
! $Z$ . 地表面気圧時間変化項.
! Surface pressure tendency
real(DP), intent(out):: xyz_QVapUAdv (0:imax-1, 1:jmax, 1:kmax)
! $uq$ . 比湿東西移流項.
! Eastward advection of specific humidity
real(DP), intent(out):: xyz_QVapVAdv (0:imax-1, 1:jmax, 1:kmax)
! $vq$ . 比湿南北移流項.
! Northward advection of specific humidity

!-----------------------------------
!  作業変数
!  Work variables
! $\Dvect{v} \cdot \nabla \pi$ .
! $\pi$ の移流. Advection of $\pi$
! $\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
! $\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_TempVir (0:imax-1, 1:jmax, 1:kmax)
! $T_v$ .
! 仮温度. Virtual temperature
real(DP):: xyz_TempVirEdd (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

! 実行文 ; Executable statement
!

! $\pi$ の移流, $\pi$ 移流の積下げ, 発散の積下げの計算
! Calculate advection of $\pi$, integral downward of advection
!   of $\pi$, integral downward of divergence
!
do k = 1, kmax
enddo

do k = kmax-1, 1, -1
enddo

xyz_DivSum(:,:,kmax) = xyz_Div(:,:,kmax) * z_DelSigma(kmax)
do k = kmax-1, 1, -1
xyz_DivSum(:,:,k) = xyz_DivSum(:,:,k+1) + xyz_Div(:,:,k) * z_DelSigma(k)
enddo

! 地表面気圧時間変化項 $Z$ の計算
! Calculate surface pressure tendency $Z$
!

select case ( LChar( trim( TimeIntegScheme ) ) )
case ('explicit')

xy_DPiDt = xy_DPiDt - xyz_DivSum(:,:,1)

end select

! $\dot{\sigma}$ の計算
! Calculate $\dot{\sigma}$
!
do k = 1, kmax-1

enddo

! $\dot{\sigma}$ の上下境界値
! $\dot{\sigma}$ on upper and lower boundary
!

! 温度の擾乱 (整数レベル), 仮温度, 仮温度の擾乱の計算
! Calculate temperature eddy (full level), virtual temperature,
!   virtual temperature eddy
!
do k = 1, kmax
xyz_TempVir(:,:,k) = xyz_Temp(:,:,k) * (   1.0_DP + ( ( ( 1.0_DP / EpsV ) - 1.0_DP ) * xyz_QVap(:,:,k) ) )

xyz_TempEdd(:,:,k) = xyz_Temp(:,:,k) - z_RefTemp(k)
xyz_TempVirEdd(:,:,k) = xyz_TempVir(:,:,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_TempInpolA(k+1) * xyz_Temp(:,:,k+1) + z_TempInpolB( k ) * xyz_Temp(:,:, k ) - r_RefTemp(k)
enddo

! 東西運動量移流項, 南北運動量移流項の計算
! Calculate advection of eastward momentum and northward momentum
!
do k = 1, kmax
xyz_UAdv(:,:,k) = ( xyz_Vor(:,:,k) + xy_Cori ) * xyz_V(:,:,k) - CpDry * z_TempInpolKappa(k) * xyz_TempVirEdd(:,:,k) * xy_GradLonPi / RPlanet

xyz_VAdv(:,:,k) = - ( xyz_Vor(:,:,k) + xy_Cori ) * xyz_U(:,:,k) - CpDry * z_TempInpolKappa(k) * xyz_TempVirEdd(:,:,k) * xy_GradLatPi / RPlanet
end do

do k = 2, kmax
xyz_UAdv(:,:,k) = xyz_UAdv(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * xyr_SigmaDot(:,:,k-1) * ( xyz_U(:,:,k-1) - xyz_U(:,:,k) )

xyz_VAdv(:,:,k) = xyz_VAdv(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * xyr_SigmaDot(:,:,k-1) * ( xyz_V(:,:,k-1) - xyz_V(:,:,k) )
end do

do k = 1, kmax-1
xyz_UAdv(:,:,k) = xyz_UAdv(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * xyr_SigmaDot(:,:,k) * ( xyz_U(:,:,k) - xyz_U(:,:,k+1) )

xyz_VAdv(:,:,k) = xyz_VAdv(:,:,k) - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) * xyr_SigmaDot(:,:,k) * ( xyz_V(:,:,k) - xyz_V(:,:,k+1) )
end do

! 運動エネルギー項 (仮温度補正含む) の計算
! Calculate kinematic energy term
!   (including virtual temperature correction)
!
call HydroGrid( xyz_TempVir - xyz_Temp, xyz_KinEngy )             ! (out)

xyz_KinEngy = xyz_KinEngy + ( xyz_U**2 + xyz_V**2 ) / 2.0_DP

! 温度東西移流項, 温度南北移流項の計算
! Calculate eastward and northward advection of temperature
!
do k = 1, kmax
end do

! 温度の時間変化項 $\hat{H}$ の計算
! Calculate temperature tendency term $\hat{H}$
!
do k = 1, kmax
xyz_DTempDt(:,:,k) = xyz_TempEdd(:,:,k) * xyz_Div(:,:,k) + z_TempInpolKappa(k) * xyz_TempVir(:,:,k) * xyz_PiAdv(:,:,k) - z_HydroAlpha(k) / z_DelSigma(k) * ( xyz_TempVir(:,:,k) * xyz_PiAdvSum(:,:,k) + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k) )
enddo

do k = 2, kmax
xyz_DTempDt(:,:,k) = xyz_DTempDt(:,:,k) - xyr_SigmaDot(:,:,k-1) / z_DelSigma(k) * ( xyr_TempEdd(:,:,k-1) - xyz_TempEdd(:,:,k) ) - xyr_SigmaDotNonGrav(:,:,k-1) / z_DelSigma(k) * ( r_RefTemp(k-1) - z_RefTemp(k) )
enddo

do k = 1, kmax-1
xyz_DTempDt(:,:,k) = xyz_DTempDt(:,:,k) - xyr_SigmaDot(:,:,k) / z_DelSigma(k) * ( xyz_TempEdd(:,:,k) - xyr_TempEdd(:,:,k) ) - xyr_SigmaDotNonGrav(:,:,k) / z_DelSigma(k) * ( z_RefTemp(k) - r_RefTemp(k) ) - z_HydroBeta(k) / z_DelSigma(k) * (   xyz_TempVir(:,:,k) * xyz_PiAdvSum(:,:,k+1) + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k+1) )
enddo

! 比湿東西移流項, 比湿南北移流項の計算
! Calculate eastward and northward advection of specific humidity
!
do k = 1, kmax
end do

! 比湿時間変化項 $R$ の計算
! Calculate specific humidity tendency $R$
!
xyz_DQVapDt = xyz_QVap * xyz_Div

do k = 2, kmax
xyz_DQVapDt(:,:,k) = xyz_DQVapDt(:,:,k) - xyr_SigmaDot(:,:,k-1) / ( 2.0_DP * z_DelSigma(k) ) * ( xyz_QVap(:,:,k-1)  - xyz_QVap(:,:,k) )
enddo

do k = 1, kmax-1
xyz_DQVapDt(:,:,k) = xyz_DQVapDt(:,:,k) - xyr_SigmaDot(:,:,k) / ( 2.0_DP * z_DelSigma(k) ) * ( xyz_QVap(:,:,k) - xyz_QVap(:,:,k+1) )
enddo

end subroutine NonLinearOnGrid

TimeIntegScheme
Variable :
TimeIntegScheme :character(TOKEN), save
 : 時間積分法. 以下の方法を選択可能. Time integration scheme. Available schemes are as follows. "Semi-implicit" "Explicit"
Subroutine :
wz_DVorDtN((nmax+1)**2, 1:kmax) :real(DP), intent(in)
 : $DD{zeta}{t} (t)$ . 渦度変化 (スペクトル). Vorticity tendency (spectral)
wz_DDivDtN((nmax+1)**2, 1:kmax) :real(DP), intent(in)
 : $DD{D}{t} (t)$ . 発散変化 (スペクトル). Divergence tendency (spectral)
wz_DTempDtN((nmax+1)**2, 1:kmax) :real(DP), intent(in)
 : $DD{T}{t} (t)$ . 温度変化 (スペクトル). Temperature tendency (spectral)
wz_DQVapDtN((nmax+1)**2, 1:kmax) :real(DP), intent(in)
 : $DD{q}{t} (t)$ . 比湿変化 (スペクトル). Specific humidity tendency (spectral)
w_DPiDtN((nmax+1)**2) :real(DP), intent(in)
 : $DD{p_s}{t} (t)$ . 地表面気圧変化 (スペクトル). Surface pressure tendency (spectral)
wz_VorB((nmax+1)**2, 1:kmax) :real(DP), intent(in)
 : $zeta (t-\Delta t)$ . 渦度 (スペクトル). Vorticity (spectral)
wz_DivB((nmax+1)**2, 1:kmax) :real(DP), intent(in)
 : $D (t-\Delta t)$ . 発散 (スペクトル). Divergence (spectral)
wz_TempB((nmax+1)**2, 1:kmax) :real(DP), intent(in)
 : $T (t-\Delta t)$ . 温度 (スペクトル). Temperature (spectral)
wz_QVapB((nmax+1)**2, 1:kmax) :real(DP), intent(in)
 : $q (t-\Delta t)$ . 比湿 (スペクトル). Specific humidity (spectral)
w_PiB((nmax+1)**2) :real(DP), intent(in)
 : $pi = ln p_s (t-\Delta t)$ . 地表面気圧 (スペクトル). Surface pressure (spectral)
w_SurfGeoPot((nmax+1)**2) :real(DP), intent(in)
 : $Phi_s$ . 地表ジオポテンシャル. Surface geo-potential
wz_VorA((nmax+1)**2, 1:kmax) :real(DP), intent(out)
 : $zeta (t+Delta t)$ . 渦度 (スペクトル). Vorticity (spectral)
wz_DivA((nmax+1)**2, 1:kmax) :real(DP), intent(out)
 : $D (t+Delta t)$ . 発散 (スペクトル). Divergence (spectral)
wz_TempA((nmax+1)**2, 1:kmax) :real(DP), intent(out)
 : $T (t+Delta t)$ . 温度 (スペクトル). Temperature (spectral)
wz_QVapA((nmax+1)**2, 1:kmax) :real(DP), intent(out)
 : $q (t+Delta t)$ . 比湿 (スペクトル). Specific humidity (spectral)
w_PiA((nmax+1)**2) :real(DP), intent(out)
 : $pi = ln p_s (t+Delta t)$ . 地表面気圧 (スペクトル). Surface pressure (spectral)

With time integration, calculate physical values at $t+Delta t$ from tendency at $t$ and physical values at $t-\Delta t$ .

Leap-frog scheme is used as time integration scheme. And forward difference is used to diffusion terms. 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".

[Source]

  subroutine TimeIntegration( wz_DVorDtN, wz_DDivDtN, wz_DTempDtN, wz_DQVapDtN, w_DPiDtN, wz_VorB,    wz_DivB,    wz_TempB,    wz_QVapB,    w_PiB, w_SurfGeoPot, wz_VorA,    wz_DivA,    wz_TempA,    wz_QVapA,    w_PiA )
!
! 時間積分を行い,
! 時刻 $t$ の物理量の時間変化と $t-\Delta t$ の物理量から
! 時刻 $t+\Delta t$ の物理量を計算します.
!
! 時間積分法にはリープフロッグスキームを用いています.
! ただし拡散項には前方差分を用いています.
! デフォルトでは, $\Delta t$ を大きくとるために, 重力波項に
! セミインプリシット法を適用しています.
! NAMELIST#dynamics_hspl_vas83_nml の TimeIntegScheme を変更することで,
! 重力波項をエクスプリシット法によって解くことも可能です.
!
! With time integration, calculate physical values at $t+\Delta t$
! from tendency at $t$ and physical values at $t-\Delta t$ .
!
! Leap-frog scheme is used as time integration scheme.
! And forward difference is used to diffusion terms.
! 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
!

! 物理定数設定
! 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 次方程式を解くための関数
! Functions to solve linear simultaneous equation by LU decomposition
!
use lumatrix, only: LUSolve

! 文字列操作
! Character handling
!
use dc_string, only: LChar

implicit none
real(DP), intent(in):: wz_DVorDtN ((nmax+1)**2, 1:kmax)
! $\DD{\zeta}{t} (t)$ . 渦度変化 (スペクトル).
! Vorticity tendency (spectral)
real(DP), intent(in):: wz_DDivDtN ((nmax+1)**2, 1:kmax)
! $\DD{D}{t} (t)$ . 発散変化 (スペクトル).
! Divergence tendency (spectral)
real(DP), intent(in):: wz_DTempDtN ((nmax+1)**2, 1:kmax)
! $\DD{T}{t} (t)$ . 温度変化 (スペクトル).
! Temperature tendency (spectral)
real(DP), intent(in):: wz_DQVapDtN ((nmax+1)**2, 1:kmax)
! $\DD{q}{t} (t)$ . 比湿変化 (スペクトル).
! Specific humidity tendency (spectral)
real(DP), intent(in):: w_DPiDtN ((nmax+1)**2)
! $\DD{p_s}{t} (t)$ . 地表面気圧変化 (スペクトル).
! Surface pressure tendency (spectral)
real(DP), intent(in):: wz_VorB ((nmax+1)**2, 1:kmax)
! $\zeta (t-\Delta t)$ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP), intent(in):: wz_DivB ((nmax+1)**2, 1:kmax)
! $D (t-\Delta t)$ . 発散 (スペクトル).
! Divergence (spectral)
real(DP), intent(in):: wz_TempB ((nmax+1)**2, 1:kmax)
! $T (t-\Delta t)$ . 温度 (スペクトル).
! Temperature (spectral)
real(DP), intent(in):: w_PiB ((nmax+1)**2)
! $\pi = \ln p_s (t-\Delta t)$ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP), intent(in):: wz_QVapB ((nmax+1)**2, 1:kmax)
! $q (t-\Delta t)$ . 比湿 (スペクトル).
! Specific humidity (spectral)
real(DP), intent(in):: w_SurfGeoPot ((nmax+1)**2)
! $\Phi_s$ . 地表ジオポテンシャル.
! Surface geo-potential

real(DP), intent(out):: wz_VorA ((nmax+1)**2, 1:kmax)
! $\zeta (t+\Delta t)$ . 渦度 (スペクトル).
! Vorticity (spectral)
real(DP), intent(out):: wz_DivA ((nmax+1)**2, 1:kmax)
! $D (t+\Delta t)$ . 発散 (スペクトル).
! Divergence (spectral)
real(DP), intent(out):: wz_TempA ((nmax+1)**2, 1:kmax)
! $T (t+\Delta t)$ . 温度 (スペクトル).
! Temperature (spectral)
real(DP), intent(out):: w_PiA ((nmax+1)**2)
! $\pi = \ln p_s (t+\Delta t)$ . 地表面気圧 (スペクトル).
! Surface pressure (spectral)
real(DP), intent(out):: wz_QVapA ((nmax+1)**2, 1:kmax)
! $q (t+\Delta t)$ . 比湿 (スペクトル).
! Specific humidity (spectral)

! 作業変数
! Work variables
!
real(DP):: wz_siTempW ((nmax+1)**2, 1:kmax)
! 温度 (セミインプリシット法のための作業変数).
! Temperature (work variable for semi-implicit scheme)
real(DP):: wz_siDTempDtW ((nmax+1)**2, 1:kmax)
! $\DD{T}{t} (t)$ . 温度変化 (スペクトル) の作業変数.
! Temperature tendency (spectral) work variable

real(DP):: w_siPiW ((nmax+1)**2)
! $\pi$ (セミインプリシット法のための作業変数).
! $\pi$ (work variable for semi-implicit scheme)
real(DP):: w_siDPiDtW ((nmax+1)**2)
! $\DD{p_s}{t} (t)$ . 地表面気圧変化 (スペクトル).
! Surface pressure tendency (spectral)
real(DP):: wz_siPhiW ((nmax+1)**2, 1:kmax)
! $\Phi = \underline{W} \overline{ \Dvect{T} }^{t}$ .
! (セミインプリシット法のための作業変数).
! (Work variable for semi-implicit scheme)
real(DP):: wz_siDivAvrTime ((nmax+1)**2, 1:kmax)
! 時間平均の $\Dvect{D}$ (セミインプリシット法のための作業変数).
! Time average $\Dvect{D}$ (work variable for semi-implicit scheme)

integer:: k, kk           ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction

! 実行文 ; Executable statement
!

! 非重力波項の仮積分
! Integration non gravitational terms temporarily
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')

wz_siTempW = wz_TempB * ( 1.0_DP - DelTime * wz_DiffTherm ) + wz_DTempDtN * DelTime

w_siPiW = w_PiB + w_DPiDtN * DelTime

end select

! ジオポテンシャルの計算
! Calculate geopotential
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')

wz_siPhiW (:,1) = CpDry * z_HydroAlpha(1) * wz_siTempW (:,1)

do k = 2, kmax
wz_siPhiW (:,k) = wz_siPhiW(:,k-1) + CpDry * z_HydroAlpha(k) * wz_siTempW (:,k) + CpDry * z_HydroBeta (k-1) * wz_siTempW (:,k-1)
end do

case ('explicit')

end select

! 発散方程式の右辺の計算
! Calculate right side of divergence equation
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')

do k = 1, kmax
wz_siDivAvrTime(:,k) = wz_DivB(:,k) * ( 1.0_DP + DelTime * wz_DiffVorDiv(:,k) ) * ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm(:,k) ) + DelTime * (   wz_DDivDtN(:,k) - wz_rn(:,k) / RPlanet**2 * (   wz_siPhiW(:,k) + ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm(:,k) ) * ( w_SurfGeoPot + w_siPiW * z_siMtxG(k) ) ) )
end do

end select

! 時間平均の $\Dvect{D}$ を LU 行列で解く
! Solve time average $\Dvect{D}$ with LU matrix
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')

wz_siDivAvrTime = LUSolve( wzz_siMtxLU, wz_siMtxPiv, wz_siDivAvrTime )

end select

! 温度, 地表気圧の時間変化項の計算
! Calculate tendency of temperature and surface pressure
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')

wz_siDTempDtW = wz_DTempDtN
do k = 1, kmax
do kk = 1, kmax
wz_siDTempDtW(:,k) = wz_siDTempDtW(:,k) - zz_siMtxH(k,kk) * wz_siDivAvrTime(:,kk)
end do
end do

w_siDPiDtW = w_DPiDtN
do kk = 1, kmax
w_siDPiDtW = w_siDPiDtW - z_DelSigma(kk) * wz_siDivAvrTime(:,kk)
end do

end select

! 時間積分. 拡散は前方差分
! Time integration. Forward difference is applied to diffusion
!

! 渦度 ; Vorticity
!
wz_VorA = ( wz_VorB + wz_DVorDtN * DelTime * 2.0_DP ) / ( 1.0_DP + DelTime * 2.0_DP * wz_DiffVorDiv )

! 発散 ; Divergence
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')

wz_DivA  = wz_siDivAvrTime * 2.0_DP - wz_DivB

case ('explicit')

wz_DivA  = ( wz_DivB + wz_DDivDtN * DelTime * 2.0_DP ) / ( 1.0_DP + DelTime * 2.0_DP * wz_DiffVorDiv )

end select

! 温度 ; Temperature
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')

wz_TempA = ( wz_TempB + wz_siDTempDtW * DelTime * 2.0_DP ) / ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm )

case ('explicit')

wz_TempA = ( wz_TempB + wz_DTempDtN * DelTime * 2.0_DP ) / ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm )

end select

! 比湿 ; Specific humidity
!
wz_QVapA = ( wz_QVapB + wz_DQVapDtN * DelTime * 2.0_DP ) / ( 1.0_DP + DelTime * 2.0_DP * wz_DiffTherm )

! 地表面気圧 ; Surface pressure
!
select case ( LChar( trim( TimeIntegScheme ) ) )
case ('semi-implicit')

w_PiA  = w_PiB + w_siDPiDtW * DelTime * 2.0_DP

case ('explicit')

w_PiA  = w_PiB + w_DPiDtN * DelTime * 2.0_DP

end select

end subroutine TimeIntegration

module_name
Constant :
module_name = ‘dynamics_hspl_vas83 :character(*), parameter
 : モジュールの名称. Module name
r_RefTemp
Variable :
r_RefTemp(:) :real(DP), allocatable
 : $hat{overline{T}}$ . 基準温度 (半整数レベル). Reference temperature on half sigma level
version
Constant :
version = ’$Name: dcpam5-20090126$’ // ’$Id: dynamics_hspl_vas83.F90,v 1.5 2009-01-23 14:41:12 morikawa Exp$’ :character(*), parameter
 : モジュールのバージョン Module version
wz_DiffTherm
Variable :
wz_DiffTherm(:,:) :real(DP), allocatable
 : $-(-1)^{N_D/2}K_{HD}nabla^{N_D}$ . 熱, 水水平拡散係数. Coefficient of horizontal thermal and water diffusion
wz_DiffVorDiv
Variable :
wz_DiffVorDiv(:,:) :real(DP), allocatable
 : $-K_{HD} [(-1)^{N_D/2}nabla^{N_D}- (2/a^2)^{N_D/2}]$ . 運動量水平拡散係数. Coefficient of horizontal momentum diffusion
wz_rn
Variable :
wz_rn(:,:) :real(DP), allocatable
 : $-n times (n+1)$ . ラプラシアンの係数. Laplacian coefficient
wz_siMtxPiv
Variable :
wz_siMtxPiv(:,:) :integer, allocatable
 : セミインプリシット行列のピボット. Pivot for semi-implicit matrix
wzz_siMtxLU
Variable :
wzz_siMtxLU(:,:,:) :real(DP), allocatable
 : セミインプリシット行列の LU 分解. LU decomposition for semi-implicit matrix
xy_Cori
Variable :
xy_Cori(:,:) :real(DP), allocatable
 : $f\equiv 2\Omega\sin\varphi$ . コリオリパラメータ. Coriolis parameter
z_HydroAlpha
Variable :
z_HydroAlpha(:) :real(DP), allocatable
 : $alpha$ . 静水圧の式の係数. Coefficient in hydrostatic equation.
z_HydroBeta
Variable :
z_HydroBeta(:) :real(DP), allocatable
 : $beta$ . 静水圧の式の係数. Coefficient in hydrostatic equation.
z_RefTemp
Variable :
z_RefTemp(:) :real(DP), allocatable
 : $overline{T}$ . 基準温度 (整数レベル). Reference temperature on full sigma level
z_TempInpolA
Variable :
z_TempInpolA(:) :real(DP), allocatable
 : $a$ . 温度鉛直補間の係数. Coefficient for vertical interpolation of temperature
z_TempInpolB
Variable :
z_TempInpolB(:) :real(DP), allocatable
 : $b$ . 温度鉛直補間の係数. Coefficient for vertical interpolation of temperature
z_TempInpolKappa
Variable :
z_TempInpolKappa(:) :real(DP), allocatable
 : $kappa$ . 温度鉛直補間の係数. Coefficient for vertical interpolation of temperature
z_siMtxG
Variable :
z_siMtxG(:) :real(DP), allocatable
 : $G = C_p kappa T$
zz_siMtxGCt
Variable :
zz_siMtxGCt(:,:) :real(DP), allocatable
 : $G C^{T}$ ( $C = Delta sigma$ )
zz_siMtxH
Variable :
zz_siMtxH(:,:) :real(DP), allocatable
 : $h = QS - R$
zz_siMtxWH
Variable :
zz_siMtxWH(:,:) :real(DP), allocatable
 : $W h$