!= 力学過程 (スペクトル法, Arakawa and Suarez (1983))
!
!= Dynamical process (Spectral method, Arakawa and Suarez (1983))
!
! Authors::   Yasuhiro MORIKAWA
! Version::   $Id: dynamics_hspl_vas83.F90,v 1.5 2009-01-23 14:41:12 morikawa Exp $ 
! Tag Name::  $Name: dcpam5-20090126 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module dynamics_hspl_vas83
  !
  != 力学過程 (スペクトル法, Arakawa and Suarez (1983))
  !
  != Dynamical Core (Spectral method, Arakawa and Suarez (1983))
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 力学過程を演算するモジュールです. 
  ! 水平離散化にスペクトル法 (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 
  ! 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
  !
  !== NAMELIST
  !
  ! NAMELIST#dynamics_hspl_vas83_nml
  !
  !== References
  !
  ! * Bourke, W. P., 1988:
  !   Spectral methods in global climate and weather 
  !   prediction models.
  !   <i> Physically Based Modelling and Simulation of
  !   Climates and Climatic Change. Part I.</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.
  !   <i>Mon. Wea. Rev.</i>, <b>111</b>, 34--35.

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

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: nmax, & ! 最大全波数. 
                             ! Maximum truncated wavenumber
    &                imax, & ! 経度格子点数. 
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数. 
                             ! Number of grid points in latitude
    &                kmax    ! 鉛直層数. 
                             ! Number of vertical level

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &  ! 倍精度実数型. Double precision. 
    &                 TOKEN  ! キーワード.   Keywords. 

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public:: Dynamics

  ! 公開変数
  ! Public variables
  !
  logical, save, public:: dynamics_hspl_vas83_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag

  ! 非公開変数
  ! Private variables
  !

  character(TOKEN), save:: TimeIntegScheme
                              ! 時間積分法. 
                              ! 以下の方法を選択可能. 
                              !
                              ! Time integration scheme. 
                              ! Available schemes are as follows.
                              ! 
                              ! * "Semi-implicit"
                              ! * "Explicit"
                              ! 

  real(DP), save:: DelTimeSave
                              ! 前回の $ \Delta t $ [s]. 
                              ! 陰解法のための係数の再生成の必要性の
                              ! チェックに使用する. 
                              ! 
                              ! $ \Delta t $ [s] at previous step
                              ! This is used to check necessity of 
                              ! regeneration of coefficients for 
                              ! implicit method. 

  ! 水平拡散係数
  ! Coefficient of horizontal diffusion
  !
  real(DP), allocatable:: wz_rn (:,:)
                              ! $ -n \times (n+1) $ . ラプラシアンの係数. 
                              ! Laplacian coefficient
  real(DP), allocatable:: wz_DiffVorDiv (:,:)
                              ! $ -K_{HD} [(-1)^{N_D/2}\nabla^{N_D}- (2/a^2)^{N_D/2}] $ . 
                              ! 運動量水平拡散係数. 
                              ! Coefficient of horizontal momentum diffusion
  real(DP), allocatable:: wz_DiffTherm (:,:)
                              ! $ -(-1)^{N_D/2}K_{HD}\nabla^{N_D} $ . 
                              ! 熱, 水水平拡散係数. 
                              ! Coefficient of horizontal thermal and water diffusion

  ! NonLinearOnGrid 等で使用する係数
  ! Coefficients for "NonLinearOnGrid", etc.
  !
  real(DP), allocatable:: xy_Cori (:,:)
                              ! $ f\equiv 2\Omega\sin\varphi $ . 
                              ! コリオリパラメータ. Coriolis parameter
  real(DP), allocatable:: z_HydroAlpha (:)
                              ! $ \alpha $ . 静水圧の式の係数. 
                              ! Coefficient in hydrostatic equation.
  real(DP), allocatable:: z_HydroBeta (:)
                              ! $ \beta $ . 静水圧の式の係数. 
                              ! Coefficient in hydrostatic equation.
  real(DP), allocatable:: z_TempInpolA (:)
                              ! $ a $ . 温度鉛直補間の係数. 
                              ! Coefficient for vertical interpolation of temperature
  real(DP), allocatable:: z_TempInpolB (:)
                              ! $ b $ . 温度鉛直補間の係数. 
                              ! Coefficient for vertical interpolation of temperature
  real(DP), allocatable:: z_TempInpolKappa (:)
                              ! $ \kappa $ . 温度鉛直補間の係数. 
                              ! Coefficient for vertical interpolation of temperature
  real(DP), allocatable:: z_RefTemp (:)
                              ! $ \overline{T} $ . 基準温度 (整数レベル). 
                              ! Reference temperature on full sigma level
  real(DP), allocatable:: r_RefTemp (:)
                              ! $ \hat{\overline{T}} $ . 基準温度 (半整数レベル). 
                              ! Reference temperature on half sigma level

  ! TimeIntegration で使用する係数
  ! Coefficients for "TimeIntegration"
  !
  real(DP), allocatable:: z_siMtxG (:)
                              ! $ G = C_p \kappa T $
  real(DP), allocatable:: zz_siMtxH (:,:)
                              ! $ h = QS - R $
  real(DP), allocatable:: zz_siMtxWH (:,:)
                              ! $ W h $
  real(DP), allocatable:: zz_siMtxGCt (:,:)
                              ! $ G C^{T} $ ( $ C = \Delta \sigma $ )

  real(DP), allocatable:: wzz_siMtxLU(:,:,:)
                              ! セミインプリシット行列の LU 分解. 
                              ! LU decomposition for semi-implicit matrix
  integer, allocatable:: wz_siMtxPiv(:,:)
                              ! セミインプリシット行列のピボット. 
                              ! Pivot for semi-implicit matrix


  character(*), parameter:: module_name = 'dynamics_hspl_vas83'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name: dcpam5-20090126 $' // &
    & '$Id: dynamics_hspl_vas83.F90,v 1.5 2009-01-23 14:41:12 morikawa Exp $'
                              ! モジュールのバージョン
                              ! Module version

  ! INTERFACE 文 ; INTERFACE statements
  !
  interface Dynamics
    module procedure Dynamics
  end interface

contains

  subroutine Dynamics( &
    & xyz_UB,   xyz_VB,   xyz_TempB,   xyz_QVapB,   xy_PsB, &   ! (in)
    & xyz_UN,   xyz_VN,   xyz_TempN,   xyz_QVapN,   xy_PsN, &   ! (in)
    & xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyz_DQVapDt, &           ! (in)
    & xy_SurfHeight, &                                          ! (in)
    & xyz_UA,   xyz_VA,   xyz_TempA,   xyz_QVapA,   xy_PsA  &   ! (out)
    & )
    !
    ! 力学過程の演算を行い, 与えられた $ 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, &
                              ! $ a $ [m]. 
                              ! 惑星半径. 
                              ! Radius of planet
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure
      & Grav
                              ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数. 
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数. 
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数. 
                               ! Number of vertical level

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $ [s]
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

#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_Lapla_wa, w_Lapla_w, wa_LaplaInv_wa, &
      & xy_GradLon_w   => xv_GradLon_w, &
      & xy_GradLat_w   => xv_GradLat_w, &
      & xya_GradLon_wa => xva_GradLon_wa, &
      & xya_GradLat_wa => xva_GradLat_wa, &
      & wa_Div_xya_xya => wa_Div_xva_xva, &
      & w_xy   => w_xv,   xy_w   => xv_w, &
      & wa_xya => wa_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_Lapla_wa, w_Lapla_w, wa_LaplaInv_wa, &
      & xy_GradLon_w, &
      & xy_GradLat_w, &
      & xya_GradLon_wa, &
      & xya_GradLat_wa, &
      & wa_Div_xya_xya, &
      & 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. 

    ! 宣言文 ; 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 $ スペクトル
    real(DP):: xy_GradLonPiN (0:imax-1, 1:jmax)
                              ! $ \frac{1}{\cos \varphi} \DP{\pi}{\lambda} $
    real(DP):: xy_GradLatPiN (0:imax-1, 1:jmax)
                              ! $ \DP{\pi}{\varphi} $

    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_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
    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_SigmaDotN (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \dot{\sigma} $ .
                              ! 鉛直流. Vertical flow
    real(DP):: xy_DPiDtN (0:imax-1, 1:jmax)
                              ! $ Z $ . 地表面気圧時間変化項. 
                              ! Surface pressure tendency
    real(DP):: xyz_QVapUAdvN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ uq (t) $ . 比湿東西移流項. 
                              ! Eastward advection of specific humidity
    real(DP):: xyz_QVapVAdvN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ 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 )
    xy_GradLonPiN = xy_GradLon_w( w_PiN )
    xy_GradLatPiN = xy_GradLat_w( w_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
    !
    call NonLinearOnGrid( &
      & xyz_UN,        xyz_VN, &        ! (in)
      & xyz_VorN,      xyz_DivN, &      ! (in)
      & xyz_TempN,     xyz_QVapN, &     ! (in)
      & xy_GradLonPiN, xy_GradLatPiN, & ! (in)
!
      & xyz_UAdvN,     xyz_VAdvN, &     ! (out)
      & xyz_DTempDtN,  xyz_DQVapDtN, &  ! (out)
      & xyz_KinEngyN, &                 ! (out)
      & xyz_TempUAdvN, xyz_TempVAdvN, & ! (out)
      & xyr_SigmaDotN, xy_DPiDtN, &     ! (out)
      & xyz_QVapUAdvN, xyz_QVapVAdvN )  ! (out)


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

    ! 渦度の時間変化 (スペクトル) の計算
    ! Calculate vorticity tendency (spectral)
    !
    wz_DVorDtN = wa_Div_xya_xya( xyz_VAdvN, - xyz_UAdvN ) / RPlanet

    ! 発散の時間変化 (スペクトル) の計算
    ! 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, & ! (in)
        &             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, &                  ! (in)
      & wz_DDivDtN + wz_DDivDtWork, &                  ! (in)
      & wz_DTempDtN + wz_DTempDtWork, &                ! (in)
      & wz_DQVapDtN + wz_DQVapDtWork, &                ! (in)
      & w_DPiDtN, &                                    ! (in)
      & wz_VorB, wz_DivB, wz_TempB, wz_QVapB, w_PiB, & ! (in)
      & w_SurfGeoPot, &                                ! (in)
      & 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, 'SigmaDot', xyr_SigmaDotN )
    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

  !-------------------------------------------------------------------

  subroutine NonLinearOnGrid( &
    & xyz_U,        xyz_V, &        ! (in)
    & xyz_Vor,      xyz_Div, &      ! (in)
    & xyz_Temp,     xyz_QVap, &     ! (in)
    & xy_GradLonPi, xy_GradLatPi, & ! (in)
!
    & xyz_UAdv,     xyz_VAdv, &     ! (out)
    & xyz_DTempDt,  xyz_DQVapDt, &  ! (out)
    & xyz_KinEngy, &                ! (out)
    & xyz_TempUAdv, xyz_TempVAdv, & ! (out)
    & xyr_SigmaDot, xy_DPiDt, &     ! (out)
    & xyz_QVapUAdv, xyz_QVapVAdv &  ! (out)
    & )
    !
    ! 非線形項 (非重力波項) を格子点上で計算します.
    !
    ! Non-linear terms (non gravitational terms) are calculated on
    ! grid points
    !

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & r_Sigma, &            ! $ \sigma $ レベル (半整数). 
                              ! Half $ \sigma $ level
      & z_DelSigma            ! $ \Delta \sigma $ (整数). 
                              ! $ \Delta \sigma $ (Full)

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & RPlanet, &
                              ! $ a $ [m]. 
                              ! 惑星半径. 
                              ! Radius of planet
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure
      & 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
    real(DP), intent(in):: xy_GradLonPi (0:imax-1, 1:jmax)
                              ! $ \frac{1}{\cos \varphi} \DP{\pi}{\lambda} $
    real(DP), intent(in):: xy_GradLatPi (0:imax-1, 1:jmax)
                              ! $ \DP{\pi}{\varphi} $
    real(DP), intent(out):: xyz_UAdv (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u_A $ . 東西運動量移流項.
                              ! Eastward advection of momentum
    real(DP), intent(out):: xyz_VAdv (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v_A $ . 南北運動量移流項. 
                              ! Northward advection of momentum
    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 $ . 温度東西移流項. 
                              ! Eastward advection of temperature
    real(DP), intent(out):: xyz_TempVAdv (0:imax-1, 1:jmax, 1:kmax)
                              ! $ vT $ . 温度南北移流項. 
                              ! Northward advection of temperature
    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
    real(DP):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \Dvect{v} \cdot \nabla \pi $ . 
                              ! $ \pi $ の移流. Advection of $ \pi $
    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_SigmaDotNonGrav (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_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
      xyz_PiAdv(:,:,k) = &
        & ( xyz_U(:,:,k) * xy_GradLonPi + xyz_V(:,:,k) * xy_GradLatPi ) &
        & / RPlanet
    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_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 $
    !
    xy_DPiDt = - xyz_PiAdvSum(:,:,1)

    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
      xyr_SigmaDot(:,:,k) = &
        & r_Sigma(k) * ( xyz_PiAdvSum(:,:,1) + xyz_DivSum(:,:,1) ) &
        & - ( xyz_PiAdvSum(:,:,k+1) + xyz_DivSum(:,:,k+1) )

      xyr_SigmaDotNonGrav(:,:,k) = &
        & r_Sigma(k) * xyz_PiAdvSum(:,:,1) - xyz_PiAdvSum(:,:,k+1)
    enddo

    ! $ \dot{\sigma} $ の上下境界値
    ! $ \dot{\sigma} $ on upper and lower boundary
    !
    xyr_SigmaDot(:,:,0)    = 0.0_DP
    xyr_SigmaDot(:,:,kmax) = 0.0_DP
    xyr_SigmaDotNonGrav(:,:,0)    = 0.0_DP
    xyr_SigmaDotNonGrav(:,:,kmax) = 0.0_DP

    ! 温度の擾乱 (整数レベル), 仮温度, 仮温度の擾乱の計算
    ! 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, & ! (in)
      &             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
      xyz_TempUAdv(:,:,k) =  xyz_U(:,:,k) * xyz_TempEdd(:,:,k)
      xyz_TempVAdv(:,:,k) =  xyz_V(:,:,k) * xyz_TempEdd(:,:,k)
    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
      xyz_QVapUAdv(:,:,k) =  xyz_U(:,:,k) * xyz_QVap(:,:,k)
      xyz_QVapVAdv(:,:,k) =  xyz_V(:,:,k) * xyz_QVap(:,:,k)
    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

  !-------------------------------------------------------------------

  subroutine HydroGrid( xyz_Temp, & ! (in)
    &                   xyz_Phi   & ! (out)
    & )
    !
    ! 格子点データである温度 $ 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( &
    & wz_DVorDtN, wz_DDivDtN, wz_DTempDtN, wz_DQVapDtN, w_DPiDtN, & ! (in)
    & wz_VorB,    wz_DivB,    wz_TempB,    wz_QVapB,    w_PiB, &    ! (in)
    & w_SurfGeoPot, &                                               ! (in)
    & wz_VorA,    wz_DivA,    wz_TempA,    wz_QVapA,    w_PiA  &    ! (out)
    & )
    !
    ! 時間積分を行い, 
    ! 時刻 $ 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, &
                              ! $ a $ [m]. 
                              ! 惑星半径. 
                              ! Radius of planet
      & 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

  !-------------------------------------------------------------------

  subroutine DiagOutput( &
    & xyz_U, xyz_V, xyz_Temp, xyz_QVap, xy_Ps & ! (in)
    & )
    !
    ! 診断量の出力を行います. 
    !
    ! Diagnostic variables are output. 
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & RPlanet, &
                              ! $ a $ [m]. 
                              ! 惑星半径. 
                              ! Radius of planet
      & Grav, &
                              ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure
      & 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, &   ! (in)
      &             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 DynamicsInit

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & RPlanet, &
                              ! $ a $ [m]. 
                              ! 惑星半径. 
                              ! Radius of planet
      & Omega, &
                              ! $ \Omega $ [s-1]. 
                              ! 回転角速度. 
                              ! Angular velocity
      & GasRDry, &
                              ! $ R $ [J kg-1 K-1]. 
                              ! 乾燥大気の気体定数. 
                              ! Gas constant of air
      & CpDry
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & z_Sigma, &            ! $ \sigma $ レベル (整数). 
                              ! Full $ \sigma $ level
      & r_Sigma, &            ! $ \sigma $ レベル (半整数). 
                              ! Half $ \sigma $ level
      & 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
    !
    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, &             ! 標準出力の装置番号. Unit number of standard output
      & 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. 
                              ! IOSTAT of NAMELIST read

    ! 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, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, &                ! (in)
        & nml = dynamics_hspl_vas83_nml, &  ! (out)
        & 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, &    ! (out)
      & 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 ImplCoef
    !
    ! TimeIntegration で使用する係数の設定. 
    !
    ! Configure coefficients for "TimeIntegration"
    !

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

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & r_Sigma, &            ! $ \sigma $ レベル (半整数). 
                              ! Half $ \sigma $ level
      & 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, &     ! (inout)
        & 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 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

end module dynamics_hspl_vas83
