!= 蜉帛ｭｦ驕守ｨ� (繧ｹ繝壹け繝医Ν豕�, Arakawa and Suarez (1983))
!
!= Dynamical process (Spectral method, Arakawa and Suarez (1983))
!
! Authors::   Yasuhiro MORIKAWA, Satoshi NODA, Yoshiyuki O. TAKAHASHI, Shin-ichi Takehiro
! Version::   $Id: dynamics_hspl_vas83.F90,v 1.83 2014/05/07 09:39:17 murashin Exp $
! Tag Name::  $Name:  $
! 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
  ! in order to enlarge the value of $ \Delta t $ .
  ! Explicit scheme can be applied to gravitational terms by changing
  ! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
  !
  !== Procedures List
  !
  ! DynamicsHSplVAS83         :: 蜉帛ｭｦ險育ｮ�
  ! DynamicsHSplVAS83Init     :: 蛻晄悄蛹�
  ! DynamicsHSplVAS83Finalize :: 邨ゆｺ��逅� (繝｢繧ｸ繝･繝ｼ繝ｫ蜀�Κ縺ｮ螟画焚縺ｮ蜑ｲ繧贋ｻ倥￠隗｣髯､)
  ! ---------------------     :: ------------
  ! DynamicsHSplVAS83         :: Calculate dynamics
  ! DynamicsHSplVAS83Init     :: Initialization
  ! DynamicsHSplVAS83Finalize :: Termination (deallocate variables in this module)
  !
  !== 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
  !

  ! 遞ｮ蛻･蝙九ヱ繝ｩ繝｡繧ｿ
  ! Kind type parameter
  !
  use dc_types, only: DP, &  ! 蛟咲ｲｾ蠎ｦ螳滓焚蝙�. Double precision.
    &                 TOKEN  ! 繧ｭ繝ｼ繝ｯ繝ｼ繝�.   Keywords.

  ! 繝｡繝�そ繝ｼ繧ｸ蜃ｺ蜉�
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 螳｣險譁� ; Declaration statements
  !
  implicit none
  private

  ! 蜈ｬ髢区焔邯壹″
  ! Public procedure
  !
  public :: DynamicsHSplVAS83
  public :: DynamicsHSplVAS83Init
  public :: DynamicsHSplVAS83Finalize

  ! 蜈ｬ髢句､画焚
  ! Public variables
  !
  logical, save, public:: wa_module_inited           = .false.

  ! 髱槫�髢句､画焚
  ! Private variables
  !
  logical, save :: dynamics_hspl_vas83_inited = .false.
                              ! 蛻晄悄險ｭ螳壹ヵ繝ｩ繧ｰ.
                              ! Initialization flag

  integer , save      :: IDTimeIntegScheme
                              ! ID for used time integration scheme

  integer , parameter :: IDTimeIntegSchemeExplicit     = 0
  integer , parameter :: IDTimeIntegSchemeSemiIMplicit = 1

  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), save, allocatable:: w_LaplaEigVal (:)
                              ! $ \nabla^2_{\sigma} = -n (n+1) $ .
                              ! 繝ｩ繝励Λ繧ｷ繧｢繝ｳ縺ｮ蝗ｺ譛牙､.
                              ! Eigen value of laplacian operator.
  real(DP), save, allocatable:: wz_DisCoefM (:,:)
                              ! 驕句虚驥上�豌ｴ蟷ｳ諡｡謨｣菫よ焚縺ｨ繧ｹ繝昴Φ繧ｸ螻､貂幄｡ｰ菫よ焚縺ｮ蜥�
                              ! Damping coefficient for momentum by horizonatl diffusion
                              ! and a sponge layer
  real(DP), save, allocatable:: wz_DisCoefH (:,:)
                              ! 辭ｱ縺ｮ豌ｴ蟷ｳ諡｡謨｣菫よ焚縺ｨ繧ｹ繝昴Φ繧ｸ螻､貂幄｡ｰ菫よ焚縺ｮ蜥�
                              ! Damping coefficient for heat by horizonatl diffusion
                              ! and a sponge layer
  real(DP), save, allocatable:: w_DisCoefQ  (:)
                              ! 謌仙�縺ｮ豌ｴ蟷ｳ諡｡謨｣菫よ焚
                              ! Damping coefficient for composition by horizonatl diffusion

  ! NonLinearOnGrid 遲峨〒菴ｿ逕ｨ縺吶ｋ菫よ焚
  ! Coefficients for "NonLinearOnGrid", etc.
  !
  real(DP), allocatable:: xy_SinLat (:,:)
                              ! $ \sin \varphi $ .
  real(DP), allocatable:: xy_CosLat (:,:)
                              ! $ \cos \varphi $ .
  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_TInpCoefA (:)
                              ! $ a $ . 貂ｩ蠎ｦ驩帷峩陬憺俣縺ｮ菫よ焚.
                              ! Coefficient for vertical interpolation of temperature
  real(DP), allocatable:: z_TInpCoefB (:)
                              ! $ b $ . 貂ｩ蠎ｦ驩帷峩陬憺俣縺ｮ菫よ焚.
                              ! Coefficient for vertical interpolation of temperature
  real(DP), allocatable:: z_TInpCoefK (:)
                              ! $ \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), save, allocatable:: z_siMtxC (:)
                              ! $ C = \sigma $
  real(DP), save, allocatable:: z_siMtxG (:)
                              ! $ G = C_p \kappa T $
  real(DP), save, allocatable:: zz_siMtxH    (:,:)
                              ! $ h = QS - R $
  real(DP), save, allocatable:: wz_siMtxDi   (:,:)
                              ! $ D = ( 1 - 2 \Delta t D_h )^{-1}$
                              ! NOTE: This is a diagonal matrix.
  real(DP), save, allocatable:: zz_siMtxDiH(:,:)
                              ! $ Di h $
  real(DP), save, allocatable:: wzz_siMtxWDiH(:,:,:)
                              ! $ W Di h $
  real(DP), save, allocatable:: zz_siMtxGCt  (:,:)
                              ! $ G C^{T} $ ( $ C = \Delta \sigma $ )

  real(DP), save, allocatable:: wzz_siMtxLU(:,:,:)
                              ! 繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ陦悟�縺ｮ LU 蛻�ｧ｣.
                              ! LU decomposition for semi-implicit matrix
  integer , save, allocatable:: wz_siMtxPiv(:,:)
                              ! 繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ陦悟�縺ｮ繝斐�繝�ヨ.
                              ! Pivot for semi-implicit matrix

  ! 繧ｵ繝悶Ν繝ｼ繝√Φ SemiImplMatrix 荳ｭ縺ｧ菴ｿ逕ｨ縺吶ｋ菫よ焚
  ! Coefficients used in a subroutine "SemiImplMatrix"
  !
  real(DP), save, allocatable:: zz_siMtxW (:,:)
                              ! $ W $ .
                              ! 逋ｺ謨｣縺ｮ蠑上〒縺ｮ邱壼ｽ｢驥榊鴨豕｢鬆��蜉ｹ譫懊↓繧医ｋ貂ｩ蠎ｦ陬懈ｭ｣縺ｮ菫よ焚.
                              ! Coefficient for correction of temperature
                              ! by effort of linear gravitational terms

  real(DP), save, allocatable:: zz_siMtxQ (:,:)
                              ! $ Q = \DD{T}{\sigma} $
  real(DP), save, allocatable:: zz_siMtxS (:,:)
                              ! $ S = \DD{\sigma}{t} $
  real(DP), save, 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, save, allocatable:: nmo (:,:,:)
!!$                              ! 繧ｹ繝壹け繝医Ν縺ｮ豺ｻ蟄鈴��分.
!!$                              ! Spectral subscript expression
!!$  real(DP), save, allocatable:: wzz_siMtxM (:,:,:)
!!$                              ! 陦悟� $ \underline{M} $.
!!$                              ! Matrix $ \underline{M} $
!!$  integer, save, allocatable:: z_siMtxPivWork(:)
!!$                              ! 陦悟�縺ｮ繝斐�繝�ヨ菴懈�縺ｮ菴懈･ｭ螟画焚.
!!$                              ! Work variable for pivot

  logical , save :: FlagDivDamp
                             ! Flag for divergence damping
  real(DP), save :: DivDampPeriod
                             ! Period for divergence damping application in unit of
                             ! second
  real(DP), save :: DivDampTime0
                             ! Time for starting divergence damping in unit of second

  logical , save :: FlagMassHorDifCor
                             ! Flag for correction of horizontal diffusion of mass
  logical , save :: FlagTotMassFix
                             ! Flag for fixing total mass

  logical , save      :: FlagSLTT

!!$  integer , save      :: IDTTMethod
!!$  integer , parameter :: IDTTMethodSp = 1
!!$  integer , parameter :: IDTTMethodSL = 2
!!$                             ! 繧ｻ繝溘Λ繧ｰ繝ｩ繝ｳ繧ｸ繝･豕輔↓繧医ｋ迚ｩ雉ｪ遘ｻ豬√ヵ繝ｩ繧ｰ縲�
!!$                             ! Flag for using Semi-Lagrangian method for tracer transport.
!!$  integer , parameter :: IDTTMethodFV = 3

  logical , save :: FlagAdvTest
                             !
                             ! Flag for advection test.
                             ! When this is true, U, V, T, and Ps are not
                             ! calculated prognostically.

  character(*), parameter:: module_name = 'dynamics_hspl_vas83'
                              ! 繝｢繧ｸ繝･繝ｼ繝ｫ縺ｮ蜷咲ｧｰ.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: dynamics_hspl_vas83.F90,v 1.83 2014/05/07 09:39:17 murashin Exp $'
                              ! 繝｢繧ｸ繝･繝ｼ繝ｫ縺ｮ繝舌�繧ｸ繝ｧ繝ｳ
                              ! Module version

contains

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

  subroutine DynamicsHSplVAS83(                                       &
    & xyz_UB,      xyz_VB,      xyz_TempB,      xyzf_QMixB,   xy_PsB, & ! (in)
    & xyz_UN,      xyz_VN,      xyz_TempN,      xyzf_QMixN,   xy_PsN, & ! (in)
    & xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy,      & ! (in)
    & xy_SurfHeight,                                                  & ! (in)
    & xyz_UA,      xyz_VA,      xyz_TempA,      xyzf_QMixA,   xy_PsA, & ! (out)
    & xyz_ArgOMG                                                      & ! (out)
    & )
    !
    ! 蜉帛ｭｦ驕守ｨ九�貍皮ｮ励ｒ陦後＞, 荳弱∴繧峨ｌ縺� $ t-\Delta t $ 縺翫ｈ縺ｳ $ t $ 縺ｮ
    ! 譚ｱ隘ｿ鬚ｨ騾� (xyz_UB, xyz_UN), 蜊怜圏鬚ｨ騾� (xyz_VB, xyz_VN),
    ! 貂ｩ蠎ｦ (xyz_TempB, xyz_TempN), 雉ｪ驥乗ｷｷ蜷域ｯ� (xyzf_QMixB, xyzf_QMixN),
    ! 蝨ｰ陦ｨ髱｢豌怜悸 (xyz_PsB, xyz_PsN), 縺翫ｈ縺ｳ蝨ｰ陦ｨ髱｢鬮伜ｺｦ (xy_SurfHeight) 縺九ｉ,
    ! $ t+\Delta t $ 縺ｮ
    ! 譚ｱ隘ｿ鬚ｨ騾� (xyz_UA), 蜊怜圏鬚ｨ騾� (xyz_VA), 貂ｩ蠎ｦ (xyz_TempA),
    ! 雉ｪ驥乗ｷｷ蜷域ｯ� (xyzf_QMixA), 蝨ｰ陦ｨ髱｢豌怜悸 (xyz_PsA) 繧定ｿ斐＠縺ｾ縺�.
    !
    ! 蛻･縺ｮ迚ｩ逅��繝ｭ繧ｻ繧ｹ縺ｫ繧医ｋ貂ｦ蠎ｦ, 逋ｺ謨｣, 貂ｩ蠎ｦ, 豈疲ｹｿ縺ｮ螟牙喧繧�,
    ! 蜉帛ｭｦ驕守ｨ九�螟牙喧縺ｫ雜ｳ縺励※谺｡縺ｮ繧ｹ繝�ャ繝励ｒ險育ｮ励☆繧句�ｴ蜷医↓縺ｯ,
    ! 縺昴ｌ繧峨�螟牙喧繧� xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy
    ! 縺ｫ荳弱∴縺ｦ縺上□縺輔＞.
    !
    ! 譎る俣遨榊�豕輔↓縺ｯ繝ｪ繝ｼ繝励ヵ繝ｭ繝�げ繧ｹ繧ｭ繝ｼ繝�繧堤畑縺�※縺�∪縺�.
    ! 繝�ヵ繧ｩ繝ｫ繝医〒縺ｯ, $ \Delta t $ 繧貞､ｧ縺阪￥縺ｨ繧九◆繧√↓, 驥榊鴨豕｢鬆�↓
    ! 繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ豕輔ｒ驕ｩ逕ｨ縺励※縺�∪縺�.
    ! NAMELIST#dynamics_hspl_vas83_nml 縺ｮ TimeIntegScheme 繧貞､画峩縺吶ｋ縺薙→縺ｧ,
    ! 驥榊鴨豕｢鬆�ｒ繧ｨ繧ｯ繧ｹ繝励Μ繧ｷ繝�ヨ豕輔↓繧医▲縺ｦ隗｣縺上％縺ｨ繧ょ庄閭ｽ縺ｧ縺�.
    !
    ! Calculate dynamical processes.
    ! Input data are
    ! Eastward wind (xyz_UB, xyz_UN),
    ! northward wind (xyz_VB, xyz_VN),
    ! temperature (xyz_TempB, xyz_TempN),
    ! mass mixing ratio (xyzf_QMixB, xyzf_QMixN),
    ! surface pressure (xyz_PsB, xyz_PsN) at $ t-\Delta t $ and $ t $,
    ! and surface height (xy_SurfHeight).
    ! Eastward wind (xyz_UA), northward wind (xyz_VA),
    ! temperature (xyz_TempA),
    ! mass mixing ratio (xyzf_QMixA), surface pressure (xyz_PsA)
    ! are returned.
    !
    ! In order to add tendencies of vorticity, divergence,
    ! temperature and specific humidity calculated by
    ! physical processes other than those by
    ! this dynamical process,
    ! these tendencies should be given to
    ! "xyz_DUDtPhy", "xyz_DVDtPhy", "xyz_DTempDtPhy", "xyzf_DQMixDtPhy"
    !
    ! Leap-frog scheme is used as time integration.
    ! By default, semi-implicit scheme is applied to gravitational terms
    ! for extension of $ \Delta t $ .
    ! Explicit scheme can be applied to gravitational terms by changing
    ! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; USE statements
    !

    use constants, only: &
      & RPlanet, &
                              ! $ 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: lmax, & ! 繧ｹ繝壹け繝医Ν繝��繧ｿ縺ｮ驟榊�繧ｵ繧､繧ｺ
                               ! Size of array for spectral data
      &                imax, & ! 邨悟ｺｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in longitude
      &                jmax, & ! 邱ｯ蠎ｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in latitude
      &                kmax    ! 驩帷峩螻､謨ｰ.
                               ! Number of vertical level

    ! 邨��縺ｫ髢｢繧上ｋ驟榊�縺ｮ險ｭ螳�
    ! Settings of array for atmospheric composition
    !
    use composition, only:                              &
      &                    ncmax,                       &
                               ! 謌仙�縺ｮ謨ｰ
                               ! Number of composition
      &                    IndexH2OVap,                 &
      &                    CompositionInqFlagAdv

    ! 譎ょ綾邂｡逅�
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $ [s]
      & TimeN, &              ! 繧ｹ繝�ャ繝� $ t $ 縺ｮ譎ょ綾. Time of step $ t $.
      & TimesetClockStart, TimesetClockStop

    ! 蠎ｧ讓吶ョ繝ｼ繧ｿ險ｭ螳�
    ! Axes data settings
    !
    use axesset, only: &
      & r_Sigma, &
                              ! $ \sigma $ 繝ｬ繝吶Ν (蜊頑紛謨ｰ).
                              ! Half $ \sigma $ level
      & z_Sigma               ! $ \sigma $ 繝ｬ繝吶Ν (謨ｴ謨ｰ).
                              ! Full $ \sigma $ level

    ! SPMODEL 繝ｩ繧､繝悶Λ繝ｪ, 逅�擇荳翫�蝠城｡後ｒ逅�擇隱ｿ蜥悟�謨ｰ螟画鋤縺ｫ繧医ｊ隗｣縺�(螟壼ｱ､蟇ｾ蠢�)
    ! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
    !
#ifdef LIB_MPI
#ifdef SJPACK
    use wa_mpi_module_sjpack, only:              &
      & wa_Lapla_wa,                             &
      & w_Lapla_w,                               &
      & wa_LaplaInv_wa,                          &
      & wa_DivLambda_xya  => wa_DivLambda_xva,   &
      & wa_DivMu_xya      => wa_DivMu_xva,       &
      & xy_GradLambda_w   => xv_GradLambda_w,    &
      & xya_GradLambda_wa => xva_GradLambda_wa,  &
      & xy_GradMu_w       => xv_GradMu_w,        &
      & xya_GradMu_wa     => xva_GradMu_wa,      &
      & w_xy              => w_xv,               &
      & xy_w              => xv_w,               &
      & wa_xya            => wa_xva,             &
      & xya_wa            => xva_wa
#else
    use wa_mpi_module, only:                     &
      & wa_Lapla_wa,                             &
      & w_Lapla_w,                               &
      & wa_LaplaInv_wa,                          &
      & wa_DivLambda_xya  => wa_DivLambda_xva,   &
      & wa_DivMu_xya      => wa_DivMu_xva,       &
      & xy_GradLambda_w   => xv_GradLambda_w,    &
      & xya_GradLambda_wa => xva_GradLambda_wa,  &
      & xy_GradMu_w       => xv_GradMu_w,        &
      & xya_GradMu_wa     => xva_GradMu_wa,      &
      & w_xy              => w_xv,               &
      & xy_w              => xv_w,               &
      & wa_xya            => wa_xva,             &
      & xya_wa            => xva_wa
#endif
#elif AXISYMMETRY
    use wa_zonal_module, only:   &
      & wa_Lapla_wa,       &
      & w_Lapla_w,         &
      & wa_LaplaInv_wa,    &
      & wa_DivLambda_xya,  &
      & wa_DivMu_xya,      &
      & xy_GradLambda_w,   &
      & xya_GradLambda_wa, &
      & xy_GradMu_w,       &
      & xya_GradMu_wa,     &
      & w_xy, xy_w,        &
      & wa_xya, xya_wa
#elif SJPACK
    use wa_module_sjpack, only:   &
      & wa_Lapla_wa,       &
      & w_Lapla_w,         &
      & wa_LaplaInv_wa,    &
      & wa_DivLambda_xya,  &
      & wa_DivMu_xya,      &
      & xy_GradLambda_w,   &
      & xya_GradLambda_wa, &
      & xy_GradMu_w,       &
      & xya_GradMu_wa,     &
      & w_xy, xy_w,        &
      & wa_xya, xya_wa
#elif AXISYMMETRY_SJPACK
    use wa_zonal_module_sjpack, only:   &
      & wa_Lapla_wa,       &
      & w_Lapla_w,         &
      & wa_LaplaInv_wa,    &
      & wa_DivLambda_xya,  &
      & wa_DivMu_xya,      &
      & xy_GradLambda_w,   &
      & xya_GradLambda_wa, &
      & xy_GradMu_w,       &
      & xya_GradMu_wa,     &
      & w_xy, xy_w,        &
      & wa_xya, xya_wa
#else
    use wa_module, only:   &
      & wa_Lapla_wa,       &
      & w_Lapla_w,         &
      & wa_LaplaInv_wa,    &
      & wa_DivLambda_xya,  &
      & wa_DivMu_xya,      &
      & xy_GradLambda_w,   &
      & xya_GradLambda_wa, &
      & xy_GradMu_w,       &
      & xya_GradMu_wa,     &
      & w_xy, xy_w,        &
      & wa_xya, xya_wa
#endif

    ! 繝偵せ繝医Μ繝��繧ｿ蜃ｺ蜉�
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 譁�ｭ怜�謫堺ｽ�
    ! Character handling
    !
    use dc_string, only: LChar

    ! 遞ｮ蛻･蝙九ヱ繝ｩ繝｡繧ｿ
    ! Kind type parameter
    !
    use dc_types, only: DP      ! 蛟咲ｲｾ蠎ｦ螳滓焚蝙�. Double precision.

    ! 遨榊�縺ｨ蟷ｳ蝮��謫堺ｽ�
    ! Operation for integral and average
    !
    use intavr_operate, only: IntLonLat_xy

    ! 雉ｪ驥上�陬懈ｭ｣
    ! Mass fixer
    !
    use mass_fixer, only: MassFixer, MassFixerColumn

    ! 繧ｻ繝溘Λ繧ｰ繝ｩ繝ｳ繧ｸ繝･豕輔↓繧医ｋ迚ｩ雉ｪ遘ｻ豬∬ｨ育ｮ�
    ! Semi-Lagrangian method for tracer transport
    !
    use sltt, only: SLTTMain

!!$    ! 迚ｩ雉ｪ遘ｻ豬� (譛蛾剞菴鍋ｩ肴ｳ�)
!!$    ! Tracer Transport (finite volume method)
!!$    !
!!$    use tt_fv, only: TTFVMain

    !
    ! Utility module for advection test
    !
    use adv_test, only: AdvTestSetHorVels, AdvTestSetVerVel

    ! 螳｣險譁� ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyz_UB    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u (t-\Delta t) $ .   譚ｱ隘ｿ鬚ｨ騾�. Eastward wind
    real(DP), intent(in):: xyz_VB    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v (t-\Delta t) $ .   蜊怜圏鬚ｨ騾�. Northward wind
    real(DP), intent(in):: xyz_TempB (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t-\Delta t) $ .   貂ｩ蠎ｦ. Temperature
    real(DP), intent(in):: xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t-\Delta t) $ .   豈疲ｹｿ. Specific humidity
    real(DP), intent(in):: xy_PsB    (0:imax-1, 1:jmax)
                              ! $ p_s (t-\Delta t) $ . 蝨ｰ陦ｨ髱｢豌怜悸. Surface pressure

    real(DP), intent(in):: xyz_UN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u (t) $ .     譚ｱ隘ｿ鬚ｨ騾�. Eastward wind
    real(DP), intent(in):: xyz_VN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v (t) $ .     蜊怜圏鬚ｨ騾�. Northward wind
    real(DP), intent(in):: xyz_TempN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t) $ .     貂ｩ蠎ｦ. Temperature
    real(DP), intent(in):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t) $ .     豈疲ｹｿ. Specific humidity
    real(DP), intent(in):: xy_PsN    (0:imax-1, 1:jmax)
                              ! $ p_s (t) $ .   蝨ｰ陦ｨ髱｢豌怜悸. Surface pressure

    real(DP), intent(in):: xyz_DUDtPhy    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \left(\DP{u}{t}\right)^{phy} $ .
                              ! 螟門鴨鬆� (迚ｩ逅�℃遞�) 縺ｫ繧医ｋ譚ｱ隘ｿ鬚ｨ騾溷､牙喧.
                              ! Eastward wind tendency by external force terms (physical processes)
    real(DP), intent(in):: xyz_DVDtPhy    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \left(\DP{v}{t}\right)^{phy} $ .
                              ! 螟門鴨鬆� (迚ｩ逅�℃遞�) 縺ｫ繧医ｋ蜊怜圏鬚ｨ騾溷､牙喧.
                              ! Northward wind tendency by external force terms (physical processes)
    real(DP), intent(in):: xyz_DTempDtPhy (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \left(\DP{T}{t}\right)^{phy} $ .
                              ! 螟門鴨鬆� (迚ｩ逅�℃遞�) 縺ｫ繧医ｋ貂ｩ蠎ｦ螟牙喧.
                              ! Temperature tendency by external force terms (physical processes)
    real(DP), intent(in):: xyzf_DQMixDtPhy(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \left(\DP{q}{t}\right)^{phy} $ .
                              ! 螟門鴨鬆� (迚ｩ逅�℃遞�) 縺ｫ繧医ｋ豈疲ｹｿ螟牙喧.
                              ! Temperature tendency by external force terms (physical processes)

    real(DP), intent(in):: xy_SurfHeight (0:imax-1, 1:jmax)
                              ! $ z_s $ . 蝨ｰ陦ｨ髱｢鬮伜ｺｦ.
                              ! Surface height.

    real(DP), intent(out):: xyz_UA    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u (t+\Delta t) $ .   譚ｱ隘ｿ鬚ｨ騾�. Eastward wind
    real(DP), intent(out):: xyz_VA    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v (t+\Delta t) $ .   蜊怜圏鬚ｨ騾�. Northward wind
    real(DP), intent(out):: xyz_TempA (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t+\Delta t) $ .   貂ｩ蠎ｦ. Temperature
    real(DP), intent(out):: xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t+\Delta t) $ .   豈疲ｹｿ. Specific humidity
    real(DP), intent(out):: xy_PsA    (0:imax-1, 1:jmax)
                              ! $ p_s (t+\Delta t) $ . 蝨ｰ陦ｨ髱｢豌怜悸. Surface pressure
    real(DP), intent(out), optional :: xyz_ArgOMG(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! OMEGA, DP/Dt

    ! 菴懈･ｭ螟画焚
    ! Work variables
    !
    real(DP):: xyz_UCosLatN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ U (t) = u (t) \cos \varphi $ .
    real(DP):: xyz_VCosLatN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ V (t) = v (t) \cos \varphi $ .
    real(DP):: xyz_VorN       (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \zeta (t) $ . 貂ｦ蠎ｦ. Vorticity
    real(DP):: xyz_DivN       (0:imax-1, 1:jmax, 1:kmax)
                              ! $ D (t) $ .     逋ｺ謨｣. Divergence

    real(DP):: xy_PiN         (0:imax-1, 1:jmax)
                              ! $ \pi = \ln p_s $
    real(DP):: wz_DivN        (lmax, 1:kmax)
                              ! $ D (t) $ . 逋ｺ謨｣ (繧ｹ繝壹け繝医Ν).
                              ! Divergence (spectral)
    real(DP):: wz_TempN       (lmax, 1:kmax)
                              ! $ T (t) $ . 貂ｩ蠎ｦ (繧ｹ繝壹け繝医Ν).
                              ! Temperature (spectral)
    real(DP):: w_PiN          (lmax)
                              ! $ \pi $ 繧ｹ繝壹け繝医Ν
    real(DP):: xy_GradLambdaPiN (0:imax-1, 1:jmax)
                              ! $ \DP{\pi}{\lambda} $
    real(DP):: xy_GradMuPiN     (0:imax-1, 1:jmax)
                              ! $ (1-\mu^2) \DP{\pi}{\mu} $

    real(DP):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \Dvect{v} \cdot \nabla \pi $ .
                              ! $ \pi $ 縺ｮ遘ｻ豬�. Advection of $ \pi $

    real(DP):: xyz_UAdvN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ U_A (t) $ . 譚ｱ隘ｿ驕句虚驥冗ｧｻ豬���.
                              ! Eastward advection of momentum
    real(DP):: xyz_VAdvN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ V_A (t) $ . 蜊怜圏驕句虚驥冗ｧｻ豬���.
                              ! Northward advection of momentum
    real(DP):: xyz_TempNonLinearN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ H (t) $ . 貂ｩ蠎ｦ譎る俣螟牙喧鬆�.
                              ! Temperature tendency
    real(DP):: xyzf_QMixNonLinearN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ R (t) $ . 豈疲ｹｿ譎る俣螟牙喧鬆�.
                              ! Specific humidity tendency
    real(DP):: xyz_KinEngyN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ KE (t) $ . 驕句虚繧ｨ繝阪Ν繧ｮ繝ｼ鬆�.
                              ! Kinetic energy
    real(DP):: xyz_TempUAdvN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ UT' (t) $ . 貂ｩ蠎ｦ譚ｱ隘ｿ遘ｻ豬���.
                              ! Eastward advection of temperature
    real(DP):: xyz_TempVAdvN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ VT' (t) $ . 貂ｩ蠎ｦ蜊怜圏遘ｻ豬���.
                              ! Northward advection of temperature
    real(DP):: xyr_SigDotN (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \dot{\sigma} $ .
                              ! 驩帷峩豬�. Vertical flow

    real(DP):: xy_DPiDtNG (0:imax-1, 1:jmax)
                              ! $ Z $ . 蝨ｰ陦ｨ髱｢豌怜悸譎る俣螟牙喧髱樣㍾蜉帶ｳ｢鬆�.
                              ! Non-gravity wave component of surface pressure tendency

    real(DP):: xyzf_QMixUAdvN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ Uq (t) $ . 豈疲ｹｿ譚ｱ隘ｿ遘ｻ豬���.
                              ! Eastward advection of specific humidity
    real(DP):: xyzf_QMixVAdvN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ Vq (t) $ . 豈疲ｹｿ蜊怜圏遘ｻ豬���.
                              ! Northward advection of specific humidity
    real(DP):: wz_DVorDtNG (lmax, 1:kmax)
                              ! $ \left( \DD{\zeta}{t} (t) \right)^{NG}$ . 貂ｦ蠎ｦ螟牙喧縺ｮ髱樣㍾蜉帶ｳ｢謌仙� (繧ｹ繝壹け繝医Ν).
                              ! Non-gravity wave component of vorticity tendency (spectral)
    real(DP):: wz_DDivDtNG (lmax, 1:kmax)
                              ! $ \left( \DD{D}{t} (t) \right)^{NG}$ . 逋ｺ謨｣螟牙喧縺ｮ髱樣㍾蜉帶ｳ｢謌仙� (繧ｹ繝壹け繝医Ν).
                              ! Non-gravity wave component of divergence tendency (spectral)
    real(DP):: wz_DTempDtNG (lmax, 1:kmax)
                              ! $ \left( \DD{T}{t} (t) \right)^{NG}$ . 貂ｩ蠎ｦ螟牙喧縺ｮ髱樣㍾蜉帶ｳ｢謌仙� (繧ｹ繝壹け繝医Ν).
                              ! Non-gravity wave component of temperature tendency (spectral)
    real(DP):: wzf_DQMixDtN(lmax, 1:kmax, 1:ncmax)
                              ! $ \DD{q}{t} (t) $ . 豈疲ｹｿ螟牙喧 (繧ｹ繝壹け繝医Ν).
                              ! Specific humidity tendency (spectral)
    real(DP):: w_DPiDtNG(lmax)
                              ! $ \left( \DD{p_s}{t} (t) \right)^{NG}$ . 蝨ｰ陦ｨ髱｢豌怜悸螟牙喧雋ｻ驥榊鴨豕｢鬆� (繧ｹ繝壹け繝医Ν).
                              ! Non-gravity wave component of surface pressure tendency (spectral)
    real(DP):: xyz_UCosLatB   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ U (t-\Delta t) = u (t-\Delta t) \cos \varphi $ .
    real(DP):: xyz_VCosLatB   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ V (t-\Delta t) = v (t-\Delta t) \cos \varphi $ .
    real(DP):: wz_VorB (lmax, 1:kmax)
                              ! $ \zeta (t-\Delta t) $ . 貂ｦ蠎ｦ (繧ｹ繝壹け繝医Ν).
                              ! Vorticity (spectral)
    real(DP):: wz_DivB (lmax, 1:kmax)
                              ! $ D (t-\Delta t) $ . 逋ｺ謨｣ (繧ｹ繝壹け繝医Ν).
                              ! Divergence (spectral)
    real(DP):: wz_TempB (lmax, 1:kmax)
                              ! $ T (t-\Delta t) $ . 貂ｩ蠎ｦ (繧ｹ繝壹け繝医Ν).
                              ! Temperature (spectral)
    real(DP):: w_PiB (lmax)
                              ! $ \pi = \ln p_s (t-\Delta t) $ . 蝨ｰ陦ｨ髱｢豌怜悸 (繧ｹ繝壹け繝医Ν).
                              ! Surface pressure (spectral)
    real(DP):: wzf_QMixB(lmax, 1:kmax, 1:ncmax)
                              ! $ q (t-\Delta t) $ . 豈疲ｹｿ (繧ｹ繝壹け繝医Ν).
                              ! Specific humidity (spectral)
    real(DP):: xyz_DUDtPhyCosLat   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \left(\DP{u}{t}\right)^{phy} \cos \varphi $ .
    real(DP):: xyz_DVDtPhyCosLat   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \left(\DP{v}{t}\right)^{phy} \cos \varphi $ .
    real(DP):: wz_VorA (lmax, 1:kmax)
                              ! $ \zeta (t+\Delta t) $ . 貂ｦ蠎ｦ (繧ｹ繝壹け繝医Ν).
                              ! Vorticity (spectral)
    real(DP):: wz_DivA (lmax, 1:kmax)
                              ! $ D (t+\Delta t) $ . 逋ｺ謨｣ (繧ｹ繝壹け繝医Ν).
                              ! Divergence (spectral)
    real(DP):: wz_TempA (lmax, 1:kmax)
                              ! $ T (t+\Delta t) $ . 貂ｩ蠎ｦ (繧ｹ繝壹け繝医Ν).
                              ! Temperature (spectral)
    real(DP):: w_PiA (lmax)
                              ! $ \pi = \ln p_s (t+\Delta t) $ . 蝨ｰ陦ｨ髱｢豌怜悸 (繧ｹ繝壹け繝医Ν).
                              ! Surface pressure (spectral)
    real(DP):: wzf_QMixA(lmax, 1:kmax, 1:ncmax)
                              ! $ q (t+\Delta t) $ . 豈疲ｹｿ (繧ｹ繝壹け繝医Ν).
                              ! Specific humidity (spectral)

    real(DP):: wz_Psi (lmax, 1:kmax)
                              ! $ \psi $ . 豬∫ｷ夐未謨ｰ. Streamline function
    real(DP):: wz_Chi (lmax, 1:kmax)
                              ! $ \chi $ . 繝昴ユ繝ｳ繧ｷ繝｣繝ｫ. Potential
    real(DP):: xyz_UCosLatA   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ U (t+\Delta t) = u (t+\Delta t) \cos \varphi $ .
    real(DP):: xyz_VCosLatA   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ V (t+\Delta t) = v (t+\Delta t) \cos \varphi $ .

    real(DP):: wz_VorDiffA (lmax, 1:kmax)
                              ! $ \mathscr{D}(\zeta) $ .
                              ! 驕句虚驥乗ｰｴ蟷ｳ諡｡謨｣縺ｫ繧医ｋ貂ｦ蠎ｦ螟牙喧 (繧ｹ繝壹け繝医Ν).
                              ! Vorticity tendency by
                              ! horizontal momentum diffusion (spectral)
    real(DP):: wz_DivDiffA (lmax, 1:kmax)
                              ! $ \mathscr{D}(D) $ .
                              ! 驕句虚驥乗ｰｴ蟷ｳ諡｡謨｣縺ｫ繧医ｋ逋ｺ謨｣螟牙喧 (繧ｹ繝壹け繝医Ν).
                              ! Divergence tendency by
                              ! horizontal momentum diffusion (spectral)
    real(DP):: wz_PsiDiff (lmax, 1:kmax)
                              ! 驕句虚驥乗ｰｴ蟷ｳ諡｡謨｣縺ｫ繧医ｋ豬∫ｷ夐未謨ｰ $ \psi $ 螟牙喧
                              ! Streamline function tendency by
                              ! horizontal momentum diffusion
    real(DP):: wz_ChiDiff (lmax, 1:kmax)
                              ! 驕句虚驥乗ｰｴ蟷ｳ諡｡謨｣縺ｫ繧医ｋ繝昴ユ繝ｳ繧ｷ繝｣繝ｫ $ \chi $ 螟牙喧
                              ! Potential tendency by
                              ! horizontal momentum diffusion
    real(DP):: xyz_UDiff (0:imax-1, 1:jmax, 1:kmax)
                              ! 驕句虚驥乗ｰｴ蟷ｳ諡｡謨｣縺ｫ繧医ｋ譚ｱ隘ｿ鬚ｨ螟牙喧.
                              ! Eastward wind tendency by
                              ! horizontal momentum diffusion
    real(DP):: xyz_VDiff (0:imax-1, 1:jmax, 1:kmax)
                              ! 驕句虚驥乗ｰｴ蟷ｳ諡｡謨｣縺ｫ繧医ｋ蜊怜圏鬚ｨ螟牙喧.
                              ! Northward wind tendency by
                              ! horizontal momentum diffusion

    ! 蝨ｰ陦ｨ髱｢鬮伜ｺｦ髢｢騾｣
    ! Surface height etc.
    !
    real(DP):: w_SurfGeoPot (lmax)
                              ! $ \Phi_s $ . 蝨ｰ陦ｨ繧ｸ繧ｪ繝昴ユ繝ｳ繧ｷ繝｣繝ｫ.
                              ! Surface geo-potential

    ! Variables for mass fixing
    !
    real(DP):: MassB
    real(DP):: MassA
    real(DP):: xyr_PressA(0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_PressB(0:imax-1, 1:jmax, 0:kmax)

    real(DP) :: xyz_PressA       (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_HDifPsA       (0:imax-1, 1:jmax)
    real(DP) :: xyz_HDifPressA   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyzf_DQMixDtHDCor(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    integer  :: kp, kn

    real(DP) :: xyz_OMG          (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! OMEGA, DP/Dt
    real(DP) :: xyzf_QMixTentative(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    logical  :: FlagExistNonAdvComp

    real(DP) :: xyz_UTracerAdv(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_VTracerAdv(0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xyr_AdvTestStreamFunc(0:imax-1, 1:jmax, 0:kmax)

    integer:: k               ! 驩帷峩譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 邨��譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in dimension of constituents

    ! 螳溯｡梧枚 ; Executable statement
    !

    ! 蛻晄悄蛹也｢ｺ隱�
    ! Initialization check
    !
    if ( .not. dynamics_hspl_vas83_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if

    ! 險育ｮ玲凾髢楢ｨ域ｸｬ髢句ｧ�
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )

    ! TimeIntegration 縺ｧ菴ｿ逕ｨ縺吶ｋ菫よ焚縺ｮ險ｭ螳�
    ! Configure coefficients for "TimeIntegration"
    !
    call SemiImplMatrix

    ! 譬ｼ蟄千せ蛟､繧偵せ繝壹け繝医Ν蛟､縺ｸ ( $ t-\Delta t$ )
    ! Exchange grid values to spectral values ( $ t-\Delta t$ )
    !
    !   貂ｦ蠎ｦ逋ｺ謨｣縺ｮ險育ｮ�
    !   Calculate vorticity and divergence
    do k = 1, kmax
      xyz_UCosLatB(:,:,k) = xyz_UB(:,:,k) * xy_CosLat
      xyz_VCosLatB(:,:,k) = xyz_VB(:,:,k) * xy_CosLat
    end do
    wz_VorB  =   (   wa_DivLambda_xya( xyz_VCosLatB )  &
      &            - wa_DivMu_xya    ( xyz_UCosLatB ) ) / RPlanet
    wz_DivB  =   (   wa_DivLambda_xya( xyz_UCosLatB )  &
      &            + wa_DivMu_xya    ( xyz_VCosLatB ) ) / RPlanet
    !
    wz_TempB = wa_xya( xyz_TempB )
    !
    w_PiB    = w_xy( log( xy_PsB ) )

    ! 譬ｼ蟄千せ蛟､繧偵せ繝壹け繝医Ν蛟､縺ｸ ( $ t $ )
    ! Exchange grid values to spectral values ( $ t $ )
    !
    !   貂ｦ蠎ｦ逋ｺ謨｣縺ｮ險育ｮ�
    !   Calculate vorticity and divergence
    !
    do k = 1, kmax
      xyz_UCosLatN(:,:,k) = xyz_UN(:,:,k) * xy_CosLat
      xyz_VCosLatN(:,:,k) = xyz_VN(:,:,k) * xy_CosLat
    end do
    xyz_VorN = xya_wa(   wa_DivLambda_xya( xyz_VCosLatN )  &
      &                - wa_DivMu_xya    ( xyz_UCosLatN ) ) / RPlanet
    wz_DivN  =       (   wa_DivLambda_xya( xyz_UCosLatN )  &
      &                + wa_DivMu_xya    ( xyz_VCosLatN ) ) / RPlanet
    xyz_DivN = xya_wa( wz_DivN )
    !
    select case ( IDTimeIntegScheme )
    case ( IDTimeIntegSchemeExplicit )
      wz_TempN = wa_xya( xyz_TempN )
    end select
    !
    xy_PiN = log( xy_PsN )
    w_PiN = w_xy( xy_PiN )
    xy_GradLambdaPiN = xy_GradLambda_w( w_PiN ) / RPlanet
    xy_GradMuPiN     = xy_GradMu_w    ( w_PiN ) / RPlanet

    ! 蝨ｰ陦ｨ繧ｸ繧ｪ繝昴ユ繝ｳ繧ｷ繝｣繝ｫ縺ｮ險育ｮ�
    ! Calculate surface geo-potential
    !
    w_SurfGeoPot = w_xy( Grav * xy_SurfHeight )

    ! 譬ｼ蟄千せ荳翫〒縺ｮ髱樒ｷ壼ｽ｢蜉帛ｭｦ鬆��險育ｮ�
    ! Calculate non-linear dynamical terms on grid points
    !
    call NonLinearOnGrid(               &
      & xyz_UCosLatN,     xyz_VCosLatN, &  ! (in)
      & xyz_VorN,         xyz_DivN, &      ! (in)
      & xyz_TempN,        xyzf_QMixN(:,:,:,IndexH2OVap),   &  ! (in)
      & xy_GradLambdaPiN, xy_GradMuPiN, &  ! (in)
!
      & xyz_PiAdv,                       & ! (out)
      & xyz_UAdvN,        xyz_VAdvN, &     ! (out)
      & xyz_TempNonLinearN, &              ! (out)
      & xyz_KinEngyN, &                    ! (out)
      & xyz_TempUAdvN,    xyz_TempVAdvN, & ! (out)
      & xyr_SigDotN,      xy_DPiDtNG  &    ! (out)
      & )

    ! 髱樒ｷ壼ｽ｢鬆�→迚ｩ逅�℃遞九↓繧医ｋ譎る俣螟牙喧鬆��蜥後ｒ繧ｹ繝壹け繝医Ν螟画鋤
    ! Spectral transformation of sum of non-linear term and tendencies by physical
    ! processes
    !
    do k = 1, kmax
      xyz_DUDtPhyCosLat(:,:,k) = xyz_DUDtPhy(:,:,k) * xy_CosLat
      xyz_DVDtPhyCosLat(:,:,k) = xyz_DVDtPhy(:,:,k) * xy_CosLat
    end do
    wz_DVorDtNG =   (   wa_DivLambda_xya( xyz_VAdvN + xyz_DVDtPhyCosLat )   &
      &               - wa_DivMu_xya    ( xyz_UAdvN + xyz_DUDtPhyCosLat ) ) &
      &             / RPlanet
    wz_DDivDtNG =   (   wa_DivLambda_xya( xyz_UAdvN + xyz_DUDtPhyCosLat )   &
      &               + wa_DivMu_xya    ( xyz_VAdvN + xyz_DVDtPhyCosLat ) ) &
      &             / RPlanet                           &
      &           - wa_Lapla_wa( wa_xya(xyz_KinEngyN) ) / RPlanet**2
    wz_DTempDtNG = - (   wa_DivLambda_xya( xyz_TempUAdvN )   &
      &                + wa_DivMu_xya    ( xyz_TempVAdvN ) ) &
      &              / RPlanet                               &
      &            + wa_xya( xyz_TempNonLinearN + xyz_DTempDtPhy )
    w_DPiDtNG = w_xy( xy_DPiDtNG )

    ! 譎る俣遨榊�
    ! Time integration
    !
    call TimeIntegration( &
      & w_SurfGeoPot, &                                      ! (in)
      & wz_DVorDtNG, wz_DDivDtNG, wz_DTempDtNG, w_DPiDtNG, & ! (in)
      & wz_DivN,                  wz_TempN,     w_PiN,     & ! (in)
      & wz_VorB,     wz_DivB,     wz_TempB,     w_PiB,     & ! (in)
      & wz_VorA,     wz_DivA,     wz_TempA,     w_PiA      & ! (out)
      & )

    ! Divergence damping
    !
    call DivergenceDamping( &
      & wz_DivA             & ! (inout)
      & )

    ! 繧ｹ繝壹け繝医Ν蛟､繧呈�ｼ蟄千せ蛟､縺ｸ ( $ t+\Delta t$ )
    ! Exchange spectral values to grid values ( $ t+\Delta t$ )
    !
    wz_Psi = wa_LaplaInv_wa( wz_VorA ) * RPlanet**2
    wz_Chi = wa_LaplaInv_wa( wz_DivA ) * RPlanet**2

    xyz_UCosLatA = &
      & (   xya_GradLambda_wa( wz_Chi ) &
      &   - xya_GradMu_wa    ( wz_Psi ) ) / RPlanet

    xyz_VCosLatA = &
      & (   xya_GradLambda_wa( wz_Psi ) &
      &   + xya_GradMu_wa    ( wz_Chi ) ) / RPlanet

    do k = 1, kmax
      xyz_UA(:,:,k) = xyz_UCosLatA(:,:,k) / xy_CosLat
      xyz_VA(:,:,k) = xyz_VCosLatA(:,:,k) / xy_CosLat
    end do
    xyz_TempA = xya_wa( wz_TempA )

    xy_PsA = exp( xy_w( w_PiA ) )

    ! 豌ｴ蟷ｳ諡｡謨｣縺ｨ繧ｹ繝昴Φ繧ｸ螻､縺ｫ繧医ｋ謨｣騾ｸ縺ｮ陬懈ｭ｣
    ! Correction for dissipation of kinetic enery by horizontal diffusion and sponge
    ! layer
    !

    ! 豌ｴ蟷ｳ諡｡謨｣縺ｨ繧ｹ繝昴Φ繧ｸ螻､縺ｫ繧医ｋ貂ｦ蠎ｦ逋ｺ謨｣縺ｮ譎る俣螟牙喧
    ! Vorticity and divergence tendency by horizontal diffusion and damping in a sponge
    ! layer
    !
    wz_VorDiffA = wz_VorA * wz_DisCoefM
    wz_DivDiffA = wz_DivA * wz_DisCoefM

    ! 豌ｴ蟷ｳ諡｡謨｣縺ｨ繧ｹ繝昴Φ繧ｸ螻､縺ｫ繧医ｋ謨｣騾ｸ縺ｮ鞫ｩ謫ｦ辭ｱ陬懈ｭ｣
    ! Correction for internal energy by horizontal diffusion and damping in a sponge
    ! layer
    !
    wz_PsiDiff = wa_LaplaInv_wa( wz_VorDiffA ) * RPlanet**2
    wz_ChiDiff = wa_LaplaInv_wa( wz_DivDiffA ) * RPlanet**2

    xyz_UDiff = &
      & (   xya_GradLambda_wa( wz_ChiDiff ) &
      &   - xya_GradMu_wa    ( wz_PsiDiff ) ) / RPlanet

    xyz_VDiff = &
      & (   xya_GradLambda_wa( wz_PsiDiff ) &
      &   + xya_GradMu_wa    ( wz_ChiDiff ) ) / RPlanet

    do k = 1, kmax
      xyz_TempA(:,:,k) = xyz_TempA(:,:,k) &
        & - (   xyz_UA(:,:,k) * xyz_UDiff(:,:,k) &
        &     + xyz_VA(:,:,k) * xyz_VDiff(:,:,k)   ) &
        &   / xy_CosLat / CpDry * 2.0_DP * DelTime
    end do

    ! Mass fixer
    !   Total Mass
    !   NOTICE: It is assumed that the total atmospheric mass does not change.
    !
    if ( FlagTotMassFix ) then
      MassB = IntLonLat_xy( xy_PsB )
      MassA = IntLonLat_xy( xy_PsA )
      if ( MassA /= 0.0_DP ) then
        xy_PsA = MassB / MassA * xy_PsA
      end if
    end if

    ! Set velocities for tracer advection
    !
    xyz_UTracerAdv = xyz_UN
    xyz_VTracerAdv = xyz_VN

    ! 荳ｻ縺ｫ遘ｻ豬√ユ繧ｹ繝医�縺溘ａ縺ｫ菴ｿ縺� if 譁�
    !
    if ( FlagAdvTest ) then
      xyz_UA      = xyz_UB
      xyz_VA      = xyz_VB
      xyz_TempA   = xyz_TempB
      xy_PsA      = xy_PsB
      xyr_SigDotN = 0.0_DP

      !
      ! Utility module for advection test
      !
      call AdvTestSetHorVels(   &
        & xyz_UTracerAdv, xyz_VTracerAdv  & ! (out)
        & )
      ! calculation of variables required in spectral advection calculation
      do k = 1, kmax
        xyz_UCosLatN(:,:,k) = xyz_UTracerAdv(:,:,k) * xy_CosLat
        xyz_VCosLatN(:,:,k) = xyz_VTracerAdv(:,:,k) * xy_CosLat
      end do
      wz_DivN  =       (   wa_DivLambda_xya( xyz_UCosLatN )  &
        &                + wa_DivMu_xya    ( xyz_VCosLatN ) ) / RPlanet
      xyz_DivN = xya_wa( wz_DivN )

      call AdvTestSetVerVel( &
        & xyr_SigDotN,              & ! (out)
        & xyr_AdvTestStreamFunc     & ! (out) optional
        & )

      call HistoryAutoPut( TimeN, 'AdvTestStreamFunc', xyr_AdvTestStreamFunc )

    end if

    if ( FlagSLTT ) then

!!$    select case ( IDTTMethod )
!!$    case ( IDTTMethodSL )
      ! Semi-Lagrangian advection

      do k = 0, kmax
        xyr_PressA(:,:,k) = xy_PsA * r_Sigma(k)
        xyr_PressB(:,:,k) = xy_PsB * r_Sigma(k)
      end do
      ! Mass fixer is included in SLTTMain
      call SLTTMain(&
        & xyr_PressB, xyr_PressA, &
        & xyz_UTracerAdv, xyz_VTracerAdv, xyr_SigDotN, & ! (in)
        & xyzf_DQMixDtPhy,             & ! (in)
        & xyzf_QMixB,                  & ! (in)
        & xyzf_QMixA                   & ! (out)
        &)

!!$    case ( IDTTMethodFV )
!!$      ! Finite volume method
!!$
!!$      do k = 0, kmax
!!$        xyr_PressA(:,:,k) = xy_PsA * r_Sigma(k)
!!$        xyr_PressB(:,:,k) = xy_PsB * r_Sigma(k)
!!$      end do
!!$      call TTFVMain(             &
!!$        & xyr_PressB, xyr_PressA, &
!!$        & xyz_UTracerAdv, xyz_VTracerAdv, xyr_SigDotN, & ! (in)
!!$        & xyzf_DQMixDtPhy,             & ! (in)
!!$        & xyzf_QMixB,                  & ! (in)
!!$        & xyzf_QMixA                   & ! (out)
!!$        & )

    else

!!$    case ( IDTTMethodSp )
      ! Spectral advection

      do n = 1, ncmax
        wzf_QMixB(:,:,n) = wa_xya( xyzf_QMixB(:,:,:,n) )
      end do
      !
      call NonLinearOnGridQMix( &
        & xyz_UCosLatN, xyz_VCosLatN, xyz_DivN, xyr_SigDotN, & ! (in)
        & xyzf_QMixN,                                        & ! (in)
        !
        & xyzf_QMixNonLinearN, &             ! (out)
        & xyzf_QMixUAdvN, xyzf_QMixVAdvN &   ! (out)
        & )

      do n = 1, ncmax
        wzf_DQMixDtN(:,:,n) = - (   wa_DivLambda_xya( xyzf_QMixUAdvN(:,:,:,n) )   &
          &                       + wa_DivMu_xya    ( xyzf_QMixVAdvN(:,:,:,n) ) ) &
          &                     / RPlanet                                         &
          &                   + wa_xya(   xyzf_QMixNonLinearN(:,:,:,n)            &
          &                             + xyzf_DQMixDtPhy    (:,:,:,n) )
      end do

      call TimeIntegrationQMix( &
        & wzf_DQMixDtN,   & ! (in)
        & wzf_QMixB,      & ! (in)
        & wzf_QMixA       & ! (out)
        & )

      do n = 1, ncmax
        xyzf_QMixA(:,:,:,n) = xya_wa( wzf_QMixA(:,:,n) )
      end do

      if ( FlagMassHorDifCor ) then

        ! Correction for horizontal diffusion of constituents.
        ! This is performed for the aim of correcting the diffusion in the viscinity of
        ! steep terrain slope.
        !

        do k = 1, kmax
          xyz_PressA(:,:,k) = xy_PsA * z_Sigma(k)
        end do
        xy_HDifPsA = xy_w( w_DisCoefQ * w_xy( xy_PsA ) )
        do k = 1, kmax
          xyz_HDifPressA(:,:,k)  = z_Sigma(k) * xy_HDifPsA(:,:)
        end do

        do n = 1, ncmax
          do k = 1, kmax
            kp = k-1
            kn = k+1
            if( k == 1    ) kp = 1
            if( k == kmax ) kn = kmax
            xyzf_DQMixDtHDCor(:,:,k,n) =                          &
              & - ( xyzf_QMixA(:,:,kn,n) - xyzf_QMixA(:,:,kp,n) ) &
              & / ( xyz_PressA(:,:,kn)   - xyz_PressA(:,:,kp)   ) &
              & * xyz_HDifPressA(:,:,k)
          end do
        end do
        xyzf_QMixA = xyzf_QMixA + xyzf_DQMixDtHDCor * 2.0_DP * DelTime
      endif

      ! Mass fixer
      !   Constituents
      !
      do k = 0, kmax
        xyr_PressA(:,:,k) = xy_PsA * r_Sigma(k)
        xyr_PressB(:,:,k) = xy_PsB * r_Sigma(k)
      end do
      call MassFixer(                  &
        & xyr_PressA,                  & ! (in)
        & xyzf_QMixA,                  & ! (inout)
        & xyr_PressRef = xyr_PressB,               & ! (in) optional
        & xyzf_QMixRef = ( xyzf_QMixB+xyzf_DQMixDtPhy*2.0_DP*DelTime ) & ! (in) optional
        & )
!!$    end select

    end if

    ! 荳ｻ縺ｫ遘ｻ豬√ユ繧ｹ繝医�縺溘ａ縺ｫ菴ｿ縺� if 譁�
    !
!!$    if ( FlagAdvTest ) then
!!$      xyr_SigDotN = 0.0_DP
!!$    end if

    ! Calculation for non-advection case
    ! 繧ｪ繧､繝ｩ繝ｼ遘ｻ豬∬ｨ育ｮ励�縺ｨ縺阪√ヴ繝ｼ繧ｭ繝ｼ縺ｪ蛻�ｸ��迚ｩ雉ｪ繧堤ｧｻ豬√＆縺帙↑縺�りｨ育ｮ怜ｮ牙ｮ壹�
    ! 縺溘ａ縲�
    !
    FlagExistNonAdvComp = .false.
    do n = 1, ncmax
      if ( .not. CompositionInqFlagAdv( n ) ) then
        xyzf_QMixA(:,:,:,n) = xyzf_QMixB(:,:,:,n) &
          & + xyzf_DQMixDtPhy(:,:,:,n) * 2.0_DP * DelTime
        xyzf_QMixTentative(:,:,:,n) = xyzf_QMixA(:,:,:,n)
        FlagExistNonAdvComp = .true.
      else
        xyzf_QMixTentative(:,:,:,n) = 0.0_DP
      end if
    end do
    if ( FlagExistNonAdvComp ) then
      ! Mass fixer
      !   Constituents
      call MassFixerColumn(            &
        & xyr_PressB,                  & ! (in)
        & xyzf_QMixTentative,          & ! (inout)
        & xyr_PressRef = xyr_PressB,   & ! (in) optional
        & xyzf_QMixRef = ( xyzf_QMixB+xyzf_DQMixDtPhy*2.0_DP*DelTime ) & ! (in) optional
        & )
      do n = 1, ncmax
        if ( .not. CompositionInqFlagAdv( n ) ) then
          xyzf_QMixA(:,:,:,n) = xyzf_QMixTentative(:,:,:,n)
        end if
      end do
    end if

    ! 繝偵せ繝医Μ繝��繧ｿ蜃ｺ蜉�
    ! History data output
    !
    call OutputDiagnosedVariables(                                 &
      & xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB,             & ! (in)
      & xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN,             & ! (in)
      & xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA,             & ! (in)
      & xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy, & ! (in)
      & xyr_SigDotN, xy_DPiDtNG,                                   & ! (in)
      & xyz_PiAdv,                                                 & ! (in)
      & xyz_OMG                                                    & ! (out)
      & )

    if ( present( xyz_ArgOMG ) ) then
      xyz_ArgOMG = xyz_OMG
    end if

    ! 險育ｮ玲凾髢楢ｨ域ｸｬ荳譎ょ●豁｢
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine DynamicsHSplVAS83

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

  subroutine NonLinearOnGrid( &
    & xyz_UCosLatN,     xyz_VCosLatN,  & ! (in)
    & xyz_VorN,         xyz_DivN,      & ! (in)
    & xyz_TempN,        xyz_QVapN,     & ! (in)
    & xy_GradLambdaPiN, xy_GradMuPiN, &  ! (in)
!
    & xyz_PiAdv,                       & ! (out)
    & xyz_UAdvN,        xyz_VAdvN, &     ! (out)
    & xyz_TempNonLinearN, &              ! (out)
    & xyz_KinEngyN, &                    ! (out)
    & xyz_TempUAdvN,    xyz_TempVAdvN, & ! (out)
    & xyr_SigDotN,      xy_DPiDtNG  &    ! (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: &
      & 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

    ! 譬ｼ蟄千せ險ｭ螳�
    ! 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

    ! 邨��縺ｫ髢｢繧上ｋ驟榊�縺ｮ險ｭ螳�
    ! Settings of array for atmospheric composition
    !
    use composition, only: &
      &                    ncmax, &
                               ! 謌仙�縺ｮ謨ｰ
                               ! Number of composition
      &                    IndexH2OVap

    implicit none
    real(DP), intent(in):: xyz_UCosLatN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ U (t) = u (t) \cos \varphi $ .
    real(DP), intent(in):: xyz_VCosLatN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ V (t) = v (t) \cos \varphi $ .
    real(DP), intent(in):: xyz_VorN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \zeta (t) $ . 貂ｦ蠎ｦ. Vorticity
    real(DP), intent(in):: xyz_DivN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ D (t) $ . 逋ｺ謨｣. Divergence
    real(DP), intent(in):: xyz_TempN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t) $ . 貂ｩ蠎ｦ. Temperature
    real(DP), intent(in):: xyz_QVapN(0:imax-1, 1:jmax, 1:kmax)
                              ! $ q (t) $ . 豈疲ｹｿ. Specific humidity
    real(DP), intent(in):: xy_GradLambdaPiN (0:imax-1, 1:jmax)
                              ! $ \DP{\pi}{\lambda} $
    real(DP), intent(in):: xy_GradMuPiN (0:imax-1, 1:jmax)
                              ! $ (1-\mu^2) \DP{\pi}{\mu} $
    real(DP), intent(out):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \Dvect{v} \cdot \nabla \pi $ .
                              ! $ \pi $ 縺ｮ遘ｻ豬�. Advection of $ \pi $
    real(DP), intent(out):: xyz_UAdvN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ U_A (t) $ . 譚ｱ隘ｿ驕句虚驥冗ｧｻ豬���.
                              ! Eastward advection of momentum
    real(DP), intent(out):: xyz_VAdvN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ V_A (t) $ . 蜊怜圏驕句虚驥冗ｧｻ豬���.
                              ! Northward advection of momentum
    real(DP), intent(out):: xyz_TempNonLinearN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \hat{H} (t) $ . 貂ｩ蠎ｦ譎る俣螟牙喧鬆�.
                              ! Temperature tendency
    real(DP), intent(out):: xyz_KinEngyN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ KE (t) $ . 驕句虚繧ｨ繝阪Ν繧ｮ繝ｼ鬆�.
                              ! Kinetic energy
    real(DP), intent(out):: xyz_TempUAdvN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ UT' (t) $ . 貂ｩ蠎ｦ譚ｱ隘ｿ遘ｻ豬���.
                              ! Eastward advection of temperature
    real(DP), intent(out):: xyz_TempVAdvN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ VT' (t) $ . 貂ｩ蠎ｦ蜊怜圏遘ｻ豬���.
                              ! Northward advection of temperature
    real(DP), intent(out):: xyr_SigDotN (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \dot{\sigma} (t) $ .
                              ! 驩帷峩豬�. Vertical flow
    real(DP), intent(out):: xy_DPiDtNG (0:imax-1, 1:jmax)
                              ! $ Z (t) $ . 蝨ｰ陦ｨ髱｢豌怜悸譎る俣螟牙喧髱樣㍾蜉帶ｳ｢鬆�.
                              ! Non-gravity wave component of surface pressure tendency

    !-----------------------------------
    !  菴懈･ｭ螟画焚
    !  Work variables
    !
    real(DP):: xyz_PiAdvSum (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \sum_k^K(\Dvect{v}\cdot\nabla\pi)\Delta\sigma $ .
                              ! $ \pi $ 遘ｻ豬√�遨榊�蛟､. Integral downward of advection of $ \pi $
    real(DP):: xyz_DivSum (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \sum_k^K D\Delta\sigma $ .
                              ! 逋ｺ謨｣縺ｮ遨榊�蛟､. Integral downward of divergence
    real(DP):: xyr_SigDotNonG (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \dot{\sigma} $ .
                              ! 驩帷峩豬� (髱樣㍾蜉帶ｳ｢). Vertical flow (non gravitational)
    real(DP):: xyz_TempEdd (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T' = T - \bar{T} $ .
                              ! 貂ｩ蠎ｦ縺ｮ謫ｾ荵ｱ (謨ｴ謨ｰ繝ｬ繝吶Ν). Temperature eddy (full level)
    real(DP):: xyr_TempEdd (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T}' $ .
                              ! 貂ｩ蠎ｦ縺ｮ謫ｾ荵ｱ (蜊頑紛謨ｰ繝ｬ繝吶Ν). Temperature eddy (half level)
    real(DP):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T_v $ .
                              ! 莉ｮ貂ｩ蠎ｦ. Virtual temperature
    real(DP):: xyz_VirTempEdd (0:imax-1, 1:jmax, 1:kmax)
                              ! $ {T_v}' = T_v - \bar{T} $ .
                              ! 莉ｮ貂ｩ蠎ｦ縺ｮ謫ｾ荵ｱ. Virtual temperature eddy

    integer:: k               ! 驩帷峩譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 邨��譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in dimension of constituents

    ! 螳溯｡梧枚 ; Executable statement
    !

    ! $ \pi $ 縺ｮ遘ｻ豬�, $ \pi $ 遘ｻ豬√→逋ｺ謨｣縺ｮ螟ｧ豌嶺ｸ顔ｫｯ縺九ｉ荳句髄縺阪�遨榊�
    ! Calculate advection of $ \pi $, integral advection of $ \pi $ and divergence
    !   from the top of the model downward
    !
    do k = 1, kmax
      xyz_PiAdv(:,:,k) = &
        & (   xyz_UCosLatN(:,:,k) * xy_GradLambdaPiN   &
        &   + xyz_VCosLatN(:,:,k) * xy_GradMuPiN     ) &
        & / ( 1.0_DP - xy_SinLat**2 )
    enddo

    xyz_PiAdvSum(:,:,kmax) = xyz_PiAdv(:,:,kmax) * z_DelSigma(kmax)
    do k = kmax-1, 1, -1
      xyz_PiAdvSum(:,:,k) = &
        & xyz_PiAdvSum(:,:,k+1) + xyz_PiAdv(:,:,k) * z_DelSigma(k)
    enddo

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

    ! 蝨ｰ陦ｨ髱｢豌怜悸譎る俣螟牙喧縺ｮ髱樣㍾蜉帶ｳ｢鬆� $ Z $ 縺ｮ險育ｮ�
    ! Calculate non-gravity wave component of surface pressure tendency $ Z $
    !
    xy_DPiDtNG = - xyz_PiAdvSum(:,:,1)

    ! $ \dot{\sigma} $ 縺ｮ險育ｮ�
    ! Calculate $ \dot{\sigma} $
    !
    do k = 1, kmax-1
      xyr_SigDotN(:,:,k) = &
        & r_Sigma(k) * ( xyz_PiAdvSum(:,:,1) + xyz_DivSum(:,:,1) ) &
        & - ( xyz_PiAdvSum(:,:,k+1) + xyz_DivSum(:,:,k+1) )

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

    ! $ \dot{\sigma} $ 縺ｮ荳贋ｸ句｢�阜蛟､
    ! $ \dot{\sigma} $ on upper and lower boundary
    !
    xyr_SigDotN   (:,:,0   ) = 0.
    xyr_SigDotN   (:,:,kmax) = 0.
    xyr_SigDotNonG(:,:,0   ) = 0.
    xyr_SigDotNonG(:,:,kmax) = 0.

    ! 貂ｩ蠎ｦ縺ｮ謫ｾ荵ｱ (謨ｴ謨ｰ繝ｬ繝吶Ν), 莉ｮ貂ｩ蠎ｦ, 莉ｮ貂ｩ蠎ｦ縺ｮ謫ｾ荵ｱ縺ｮ險育ｮ�
    ! Calculate temperature eddy (full level), virtual temperature,
    !   virtual temperature eddy
    !
    do k = 1, kmax
      xyz_VirTemp(:,:,k) = &
        & xyz_TempN(:,:,k) &
        &   * ( 1.0_DP + ((( 1.0_DP / EpsV ) - 1.0_DP ) * xyz_QVapN(:,:,k)) )

      xyz_TempEdd(:,:,k) = xyz_TempN(:,:,k) - z_RefTemp(k)
      xyz_VirTempEdd(:,:,k) = xyz_VirTemp(:,:,k) - z_RefTemp(k)
    enddo

    ! 貂ｩ蠎ｦ縺ｮ謫ｾ荵ｱ (蜊頑紛謨ｰ繝ｬ繝吶Ν) 縺ｮ險育ｮ�
    ! Calculate temperature eddy (half level)
    !
    xyr_TempEdd(:,:,0   ) = 0.0_DP
    xyr_TempEdd(:,:,kmax) = 0.0_DP
    do k = 1, kmax-1
      xyr_TempEdd(:,:,k) = &
        &   z_TInpCoefA(k+1) * xyz_TempN(:,:,k+1) &
        & + z_TInpCoefB(k  ) * xyz_TempN(:,:,k  ) &
        & - r_RefTemp(k)
    enddo

    ! 譚ｱ隘ｿ驕句虚驥冗ｧｻ豬����險育ｮ�
    ! Calculate advection of eastward momentum
    !
    xyz_UAdvN(:,:,1) = &
      &   ( xyz_VorN(:,:,1) + xy_Cori ) * xyz_VCosLatN(:,:,1) &
!
      & - 1.0_DP / ( 2.0_DP * z_DelSigma(1) ) * xyr_SigDotN(:,:,1) &
      &                               * ( xyz_UCosLatN(:,:,1) - xyz_UCosLatN(:,:,2) ) &
!
      & - CpDry * z_TInpCoefK(1) * xyz_VirTempEdd(:,:,1) * xy_GradLambdaPiN

    do k = 2, kmax-1
      xyz_UAdvN(:,:,k) = &
        &   ( xyz_VorN(:,:,k) + xy_Cori ) * xyz_VCosLatN(:,:,k) &
!
        & - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) &
        &   * ( &
        &        xyr_SigDotN(:,:,k-1) * ( xyz_UCosLatN(:,:,k-1) - xyz_UCosLatN(:,:,k) ) &
        &      + xyr_SigDotN(:,:,k)   * ( xyz_UCosLatN(:,:,k)   - xyz_UCosLatN(:,:,k+1) ) &
        &     ) &
!
        & - CpDry * z_TInpCoefK(k) * xyz_VirTempEdd(:,:,k) * xy_GradLambdaPiN
    end do

    xyz_UAdvN(:,:,kmax) = &
      &   ( xyz_VorN(:,:,kmax) + xy_Cori ) * xyz_VCosLatN(:,:,kmax) &
!
      & - 1.0_DP / ( 2.0_DP * z_DelSigma(kmax) ) * xyr_SigDotN(:,:,kmax-1) &
      &                                  * ( xyz_UCosLatN(:,:,kmax-1) - xyz_UCosLatN(:,:,kmax) ) &
!
      & - CpDry * z_TInpCoefK(kmax) * xyz_VirTempEdd(:,:,kmax) * xy_GradLambdaPiN

    ! 蜊怜圏驕句虚驥冗ｧｻ豬����險育ｮ�
    ! Calculate advection of northward momentum
    !
    xyz_VAdvN(:,:,1) = &
      & - ( xyz_VorN(:,:,1) + xy_Cori ) * xyz_UCosLatN(:,:,1) &
!
      & - 1.0_DP / ( 2.0_DP * z_DelSigma(1) ) * xyr_SigDotN(:,:,1)&
      &                               * ( xyz_VCosLatN(:,:,1) - xyz_VCosLatN(:,:,2) ) &
!
      & - CpDry * z_TInpCoefK(1) &
      &          * xyz_VirTempEdd(:,:,1) * xy_GradMuPiN

    do k = 2, kmax-1
      xyz_VAdvN(:,:,k) = &
        & - ( xyz_VorN(:,:,k) + xy_Cori ) * xyz_UCosLatN(:,:,k) &
!
        & - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) &
        &   * ( &
        &        xyr_SigDotN(:,:,k-1) * ( xyz_VCosLatN(:,:,k-1) - xyz_VCosLatN(:,:,k) ) &
        &      + xyr_SigDotN(:,:,k)   * ( xyz_VCosLatN(:,:,k)   - xyz_VCosLatN(:,:,k+1) ) &
        &     ) &
!
        & - CpDry * z_TInpCoefK(k) * xyz_VirTempEdd(:,:,k) * xy_GradMuPiN
    end do

    xyz_VAdvN(:,:,kmax) = &
      & - ( xyz_VorN(:,:,kmax) + xy_Cori ) * xyz_UCosLatN(:,:,kmax) &
!
      & - 1.0_DP / ( 2.0_DP * z_DelSigma(kmax) ) * xyr_SigDotN(:,:,kmax-1) &
      &                                  * ( xyz_VCosLatN(:,:,kmax-1) - xyz_VCosLatN(:,:,kmax) ) &
!
      & - CpDry * z_TInpCoefK(kmax) * xyz_VirTempEdd(:,:,kmax) * xy_GradMuPiN

    ! 驕句虚繧ｨ繝阪Ν繧ｮ繝ｼ鬆� (莉ｮ貂ｩ蠎ｦ陬懈ｭ｣蜷ｫ繧) 縺ｮ險育ｮ�
    ! Calculate kinematic energy term
    !   (including virtual temperature correction)
    !
    call HydroGrid( xyz_VirTemp - xyz_TempN, & ! (in)
      &             xyz_KinEngyN )             ! (out)

    do k = 1, kmax
      xyz_KinEngyN(:,:,k) =                                     &
        &   xyz_KinEngyN(:,:,k)                                 &
        & + ( xyz_UCosLatN(:,:,k)**2 + xyz_VCosLatN(:,:,k)**2 ) &
        &       / ( 2.0_DP * ( 1.0_DP - xy_SinLat**2 ) )
    end do

    ! 貂ｩ蠎ｦ譚ｱ隘ｿ遘ｻ豬���, 貂ｩ蠎ｦ蜊怜圏遘ｻ豬����險育ｮ�
    ! Calculate eastward and northward advection of temperature
    !
    do k = 1, kmax
      xyz_TempUAdvN(:,:,k) =  xyz_UCosLatN(:,:,k) * xyz_TempEdd(:,:,k)
      xyz_TempVAdvN(:,:,k) =  xyz_VCosLatN(:,:,k) * xyz_TempEdd(:,:,k)
    end do

    ! 貂ｩ蠎ｦ縺ｮ譎る俣螟牙喧鬆� $ \hat{H} $ 縺ｮ險育ｮ�
    ! Calculate temperature tendency term $ \hat{H} $
    !
    do k = 1, kmax-1
      xyz_TempNonLinearN(:,:,k) = &
        &   xyz_TempEdd(:,:,k) * xyz_DivN(:,:,k) &
!
        & - 1.0_DP / z_DelSigma(k) &
        &  * (   xyr_SigDotN(:,:,k-1) &
        &          * ( xyr_TempEdd(:,:,k-1) - xyz_TempEdd(:,:,k) ) &
        &      + xyr_SigDotN(:,:,k) &
        &          * ( xyz_TempEdd(:,:,k)   - xyr_TempEdd(:,:,k) ) ) &
!
        & - 1.0_DP / z_DelSigma(k) &
        &  * (   xyr_SigDotNonG(:,:,k-1) &
        &          * ( r_RefTemp(k-1) - z_RefTemp(k) ) &
        &      + xyr_SigDotNonG(:,:,k) &
        &          * ( z_RefTemp(k)   - r_RefTemp(k) ) ) &
!
        & + z_TInpCoefK(k) * xyz_VirTemp(:,:,k) * xyz_PiAdv(:,:,k) &
!
        & - z_HydroAlpha(k) / z_DelSigma(k) &
        &    * (   xyz_VirTemp(:,:,k)    * xyz_PiAdvSum(:,:,k) &
        &        + xyz_VirTempEdd(:,:,k) * xyz_DivSum(:,:,k) ) &
!
        & - z_HydroBeta(k) / z_DelSigma(k) &
        &    * (   xyz_VirTemp(:,:,k)    * xyz_PiAdvSum(:,:,k+1) &
        &        + xyz_VirTempEdd(:,:,k) * xyz_DivSum(:,:,k+1) )
    enddo

    xyz_TempNonLinearN(:,:,kmax) = &
      &   xyz_TempEdd(:,:,kmax) * xyz_DivN(:,:,kmax) &
!
      & - 1.0_DP / z_DelSigma(kmax) &
      &  * (   xyr_SigDotN(:,:,kmax-1) &
      &          * ( xyr_TempEdd(:,:,kmax-1) - xyz_TempEdd(:,:,kmax) ) &
      &      + xyr_SigDotN(:,:,kmax) &
      &          * ( xyz_TempEdd(:,:,kmax)   - xyr_TempEdd(:,:,kmax) ) ) &
!
      & - 1.0_DP / z_DelSigma(kmax) &
      &  * (   xyr_SigDotNonG(:,:,kmax-1) &
      &          * ( r_RefTemp(kmax-1) - z_RefTemp(kmax) ) &
      &      + xyr_SigDotNonG(:,:,kmax) &
      &          * ( z_RefTemp(kmax)   - r_RefTemp(kmax) ) ) &
!
      & + z_TInpCoefK(kmax) * xyz_VirTemp(:,:,kmax) * xyz_PiAdv(:,:,kmax) &
!
      & - z_HydroAlpha(kmax) / z_DelSigma(kmax) &
      &    * (   xyz_VirTemp(:,:,kmax)    * xyz_PiAdvSum(:,:,kmax) &
      &        + xyz_VirTempEdd(:,:,kmax) * xyz_DivSum(:,:,kmax) )

  end subroutine NonLinearOnGrid

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

  subroutine NonLinearOnGridQMix( &
    & xyz_UCosLatN, xyz_VCosLatN, xyz_DivN, xyr_SigDotN, & ! (in)
    & xyzf_QMixN,                                        & ! (in)
!
    & xyzf_QMixNonLinearN, &             ! (out)
    & xyzf_QMixUAdvN, xyzf_QMixVAdvN &   ! (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)

    ! 譬ｼ蟄千せ險ｭ螳�
    ! 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

    ! 邨��縺ｫ髢｢繧上ｋ驟榊�縺ｮ險ｭ螳�
    ! Settings of array for atmospheric composition
    !
    use composition, only: &
      &                    ncmax
                               ! 謌仙�縺ｮ謨ｰ
                               ! Number of composition

    implicit none
    real(DP), intent(in):: xyz_UCosLatN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ U (t) = u (t) \cos \varphi $ .
    real(DP), intent(in):: xyz_VCosLatN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ V (t) = v (t) \cos \varphi $ .
    real(DP), intent(in):: xyz_DivN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ D (t) $ . 逋ｺ謨｣. Divergence
    real(DP), intent(in ):: xyr_SigDotN (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \dot{\sigma} (t) $ .
                              ! 驩帷峩豬�. Vertical flow
    real(DP), intent(in):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t) $ . 豈疲ｹｿ. Specific humidity
    real(DP), intent(out):: xyzf_QMixNonLinearN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ R (t) $ . 豈疲ｹｿ譎る俣螟牙喧鬆�.
                              ! Specific humidity tendency
    real(DP), intent(out):: xyzf_QMixUAdvN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ Uq (t) $ . 豈疲ｹｿ譚ｱ隘ｿ遘ｻ豬���.
                              ! Eastward advection of specific humidity
    real(DP), intent(out):: xyzf_QMixVAdvN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ Vq (t) $ . 豈疲ｹｿ蜊怜圏遘ｻ豬���.
                              ! Northward advection of specific humidity

    !-----------------------------------
    !  菴懈･ｭ螟画焚
    !  Work variables
    !
    integer:: k               ! 驩帷峩譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 邨��譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in dimension of constituents

    ! 螳溯｡梧枚 ; Executable statement
    !

    ! 豈疲ｹｿ譚ｱ隘ｿ遘ｻ豬���, 豈疲ｹｿ蜊怜圏遘ｻ豬����險育ｮ�
    ! Calculate eastward and northward advection of specific humidity
    !
    do n = 1, ncmax
      do k = 1, kmax
        xyzf_QMixUAdvN(:,:,k,n) =  xyz_UCosLatN(:,:,k) * xyzf_QMixN(:,:,k,n)
        xyzf_QMixVAdvN(:,:,k,n) =  xyz_VCosLatN(:,:,k) * xyzf_QMixN(:,:,k,n)
      end do
    end do

    ! 豈疲ｹｿ譎る俣螟牙喧鬆� $ R $ 縺ｮ險育ｮ�
    ! Calculate specific humidity tendency $ R $
    !
    do n = 1, ncmax
      xyzf_QMixNonLinearN(:,:,1,n) =   &
        &    xyzf_QMixN(:,:,1,n) * xyz_DivN(:,:,1) &
        &  - 1.0_DP / ( 2.0_DP * z_DelSigma(1) ) &
        &       * xyr_SigDotN(:,:,1) &
        &       * ( xyzf_QMixN(:,:,1,n) - xyzf_QMixN(:,:,2,n) )

      do k = 2, kmax - 1
        xyzf_QMixNonLinearN(:,:,k,n) =   &
          &    xyzf_QMixN(:,:,k,n) * xyz_DivN(:,:,k) &
          &  - 1.0_DP / ( 2.0_DP * z_DelSigma(k) ) &
          &       * ( &
          &             xyr_SigDotN(:,:,k-1) &
          &               * ( xyzf_QMixN(:,:,k-1,n) - xyzf_QMixN(:,:,k  ,n)   ) &
          &           + xyr_SigDotN(:,:,k) &
          &               * ( xyzf_QMixN(:,:,k  ,n) - xyzf_QMixN(:,:,k+1,n) ) )
      end do

      xyzf_QMixNonLinearN(:,:,kmax,n) =   &
        &    xyzf_QMixN(:,:,kmax,n) * xyz_DivN(:,:,kmax) &
        &  - 1.0_DP / ( 2.0_DP * z_DelSigma(kmax) ) &
        &       * xyr_SigDotN(:,:,kmax-1) &
        &       * ( xyzf_QMixN(:,:,kmax-1,n) - xyzf_QMixN(:,:,kmax,n) )
    end do

  end subroutine NonLinearOnGridQMix

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

  subroutine 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

    ! 譬ｼ蟄千せ險ｭ螳�
    ! 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

    ! 螳｣險譁� ; 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( &
    & w_SurfGeoPot, &                                      ! (in)
    & wz_DVorDtNG, wz_DDivDtNG, wz_DTempDtNG, w_DPiDtNG, & ! (in)
    & wz_DivN,                  wz_TempN,     w_PiN,     & ! (in)
    & wz_VorB,     wz_DivB,     wz_TempB,     w_PiB,     & ! (in)
    & wz_VorA,     wz_DivA,     wz_TempA,     w_PiA      & ! (out)
    & )
    !
    ! 譎る俣遨榊�繧定｡後＞,
    ! 譎ょ綾 $ t $ 縺ｮ迚ｩ逅�㍼縺ｮ譎る俣螟牙喧縺ｨ $ t-\Delta t$ 縺ｮ迚ｩ逅�㍼縺九ｉ
    ! 譎ょ綾 $ t+\Delta t $ 縺ｮ迚ｩ逅�㍼繧定ｨ育ｮ励＠縺ｾ縺�.
    !
    ! 譎る俣遨榊�豕輔↓縺ｯ繝ｪ繝ｼ繝励ヵ繝ｭ繝�げ繧ｹ繧ｭ繝ｼ繝�繧堤畑縺�※縺�∪縺�.
    ! 縺溘□縺玲僑謨｣鬆�↓縺ｯ蠕梧婿蟾ｮ蛻�ｒ逕ｨ縺�※縺�∪縺�.
    ! 繝�ヵ繧ｩ繝ｫ繝医〒縺ｯ, $ \Delta t $ 繧貞､ｧ縺阪￥縺ｨ繧九◆繧√↓, 驥榊鴨豕｢鬆�↓
    ! 繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ豕輔ｒ驕ｩ逕ｨ縺励※縺�∪縺�.
    ! NAMELIST#dynamics_hspl_vas83_nml 縺ｮ TimeIntegScheme 繧貞､画峩縺吶ｋ縺薙→縺ｧ,
    ! 驥榊鴨豕｢鬆�ｒ繧ｨ繧ｯ繧ｹ繝励Μ繧ｷ繝�ヨ豕輔↓繧医▲縺ｦ隗｣縺上％縺ｨ繧ょ庄閭ｽ縺ｧ縺�.
    !
    ! With time integration, physical values at $ t+\Delta t $ is calculated
    ! from tendency at $ t $ and physical values at $ t-\Delta t $ .
    !
    ! Leap-frog scheme is used as time integration scheme.
    ! And backward scheme is used for diffusion terms.
    ! As default setting, semi-implicit scheme is applied to gravitational terms
    ! in order to enlarge the value of $ \Delta t $ .
    ! Explicit scheme can be applied to gravitational terms by changing
    ! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; 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 谺｡譁ｹ遞句ｼ上ｒ隗｣縺上◆繧√�髢｢謨ｰ (spml 蜷梧｢ｱ繝｢繧ｸ繝･繝ｼ繝ｫ)
    ! Functions to solve linear simultaneous equation by LU decomposition
    ! (a module included in spml)
    !
    use lumatrix, only: LUSolve

    ! 譁�ｭ怜�謫堺ｽ�
    ! Character handling
    !
    use dc_string, only: LChar

    ! 譬ｼ蟄千せ險ｭ螳�
    ! Grid points settings
    !
    use gridset, only: lmax, & ! 繧ｹ繝壹け繝医Ν繝��繧ｿ縺ｮ驟榊�繧ｵ繧､繧ｺ
                               ! Size of array for spectral data
      &                imax, & ! 邨悟ｺｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in longitude
      &                jmax, & ! 邱ｯ蠎ｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in latitude
      &                kmax    ! 驩帷峩螻､謨ｰ.
                               ! Number of vertical level

    ! 邨��縺ｫ髢｢繧上ｋ驟榊�縺ｮ險ｭ螳�
    ! Settings of array for atmospheric composition
    !
    use composition, only: &
      &                    ncmax, &
                               ! 謌仙�縺ｮ謨ｰ
                               ! Number of composition
      &                    IndexH2OVap

    implicit none
    real(DP), intent(in):: w_SurfGeoPot (lmax)
                              ! $ \Phi_s $ . 蝨ｰ陦ｨ繧ｸ繧ｪ繝昴ユ繝ｳ繧ｷ繝｣繝ｫ.
                              ! Surface geo-potential

    real(DP), intent(in):: wz_DVorDtNG (lmax, 1:kmax)
                              ! $ \left( \DD{\zeta}{t} (t) \right)^{NG}$ . 貂ｦ蠎ｦ螟牙喧縺ｮ髱樣㍾蜉帶ｳ｢謌仙� (繧ｹ繝壹け繝医Ν).
                              ! Non-gravity wave component of vorticity tendency (spectral)
    real(DP), intent(in):: wz_DDivDtNG (lmax, 1:kmax)
                              ! $ \DD{D}{t} (t) $ . 逋ｺ謨｣螟牙喧縺ｮ髱樣㍾蜉帶ｳ｢謌仙� (繧ｹ繝壹け繝医Ν).
                              ! Divergence tendency (spectral)
    real(DP), intent(in):: wz_DTempDtNG(lmax, 1:kmax)
                              ! $ \left( \DD{T}{t} (t) \right)^{NG}$ . 貂ｩ蠎ｦ螟牙喧縺ｮ髱樣㍾蜉帶ｳ｢謌仙� (繧ｹ繝壹け繝医Ν).
                              ! Temperature tendency (spectral)
    real(DP), intent(in):: w_DPiDtNG(lmax)
                              ! $ \left( \DD{p_s}{t} (t) \right) $ . 蝨ｰ陦ｨ髱｢豌怜悸螟牙喧縺ｮ髱樣㍾蜉帶ｳ｢鬆� (繧ｹ繝壹け繝医Ν).
                              ! Non-gravity wave component of surface pressure tendency (spectral)
    real(DP), intent(in):: wz_DivN (lmax, 1:kmax)
                              ! $ D (t) $ . 逋ｺ謨｣ (繧ｹ繝壹け繝医Ν).
                              ! Divergence (spectral)
    real(DP), intent(in):: wz_TempN (lmax, 1:kmax)
                              ! $ T (t) $ . 貂ｩ蠎ｦ (繧ｹ繝壹け繝医Ν).
                              ! Temperature (spectral)
    real(DP), intent(in):: w_PiN   (lmax)
                              ! $ \pi = \ln p_s (t) $ . 蝨ｰ陦ｨ髱｢豌怜悸 (繧ｹ繝壹け繝医Ν).

    real(DP), intent(in):: wz_VorB (lmax, 1:kmax)
                              ! $ \zeta (t-\Delta t) $ . 貂ｦ蠎ｦ (繧ｹ繝壹け繝医Ν).
                              ! Vorticity (spectral)
    real(DP), intent(in):: wz_DivB (lmax, 1:kmax)
                              ! $ D (t-\Delta t) $ . 逋ｺ謨｣ (繧ｹ繝壹け繝医Ν).
                              ! Divergence (spectral)
    real(DP), intent(in):: wz_TempB (lmax, 1:kmax)
                              ! $ T (t-\Delta t) $ . 貂ｩ蠎ｦ (繧ｹ繝壹け繝医Ν).
                              ! Temperature (spectral)
    real(DP), intent(in):: w_PiB (lmax)
                              ! $ \pi = \ln p_s (t-\Delta t) $ . 蝨ｰ陦ｨ髱｢豌怜悸 (繧ｹ繝壹け繝医Ν).
                              ! Surface pressure (spectral)

    real(DP), intent(out):: wz_VorA (lmax, 1:kmax)
                              ! $ \zeta (t+\Delta t) $ . 貂ｦ蠎ｦ (繧ｹ繝壹け繝医Ν).
                              ! Vorticity (spectral)
    real(DP), intent(out):: wz_DivA (lmax, 1:kmax)
                              ! $ D (t+\Delta t) $ . 逋ｺ謨｣ (繧ｹ繝壹け繝医Ν).
                              ! Divergence (spectral)
    real(DP), intent(out):: wz_TempA (lmax, 1:kmax)
                              ! $ T (t+\Delta t) $ . 貂ｩ蠎ｦ (繧ｹ繝壹け繝医Ν).
                              ! Temperature (spectral)
    real(DP), intent(out):: w_PiA (lmax)
                              ! $ \pi = \ln p_s (t+\Delta t) $ . 蝨ｰ陦ｨ髱｢豌怜悸 (繧ｹ繝壹け繝医Ν).
                              ! Surface pressure (spectral)

    ! 菴懈･ｭ螟画焚
    ! Work variables
    !
    real(DP):: wz_HDiv (lmax, 1:kmax)
    real(DP):: w_CtDiv (lmax)
    real(DP):: wz_WT   (lmax, 1:kmax)

    real(DP):: wz_siTemp (lmax, 1:kmax)
                              ! 貂ｩ蠎ｦ (繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ豕輔�縺溘ａ縺ｮ菴懈･ｭ螟画焚).
                              ! Temperature (work variable for semi-implicit scheme)
    real(DP):: w_siPi (lmax)
                              ! $ \pi $ (繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ豕輔�縺溘ａ縺ｮ菴懈･ｭ螟画焚).
                              ! $ \pi $ (work variable for semi-implicit scheme)
    real(DP):: wz_siPhi (lmax, 1:kmax)
                              ! $ \Phi = \underline{W} ( 1 - 2 \Delta t \underline{D_H} ) \overline{ \Dvect{T} }^{t}$ .
                              ! (繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ豕輔�縺溘ａ縺ｮ菴懈･ｭ螟画焚).
                              ! (Work variable for semi-implicit scheme)
    real(DP):: wz_siVectF (lmax, 1:kmax)
                              ! $ \Dvect{f} $ .
                              ! 逋ｺ謨｣鬆�↓髢｢縺吶ｋ繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ譁ｹ遞句ｼ上�蜿ｳ霎ｺ.
                              ! Right-hand side of a semi-implicit equation of a divergence term.

    real(DP):: wz_siDivAvrTime (lmax, 1:kmax)
                              ! $ \overline{\Dvect{D}}^{t} $ .
                              ! 譎る俣蟷ｳ蝮�� $ \Dvect{D} $ (繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ豕輔�縺溘ａ縺ｮ菴懈･ｭ螟画焚).
                              ! Time average $ \Dvect{D} $ (a work variable for semi-implicit scheme)

    integer:: k, kk           ! 驩帷峩譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 邨��譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in dimension of constituents

    ! 螳溯｡梧枚 ; Executable statement
    !

    ! 譎る俣遨榊�. 諡｡謨｣縺ｯ蠕梧婿蟾ｮ蛻�
    ! Time integration. Backward difference is applied to diffusion
    !

    ! 貂ｦ蠎ｦ ; Vorticity
    !
    wz_VorA = &
      & ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM ) ) &
      &   * ( wz_VorB + 2.0_DP * DelTime * wz_DVorDtNG )

    ! 逋ｺ謨｣ ; Divergence
    !
    select case ( IDTimeIntegScheme )
    case ( IDTimeIntegSchemeSemiImplicit )

      ! 繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ豕輔〒逕ｨ縺�ｋ陦悟�縺ｮ險育ｮ励↓菴ｿ繧上ｌ繧矩���險育ｮ�
      ! Calculate terms used in making a matrix for semi-implicit method
      !
      wz_siTemp = &
        &   ( 1.0_DP - DelTime * wz_DisCoefH ) * wz_TempB &
        & + DelTime * wz_DTempDtNG

      !     NOTE:
      !     The matrix wz_siMtxDi = ( 1 - 2 \Delta t D_H )^{-1} is diagonal matrix.
      !
      wz_siPhi (:,1) = CpDry * z_HydroAlpha(1) * wz_siMtxDi(:,1) * wz_siTemp(:,1)
      do k = 2, kmax
        wz_siPhi (:,k) = wz_siPhi(:,k-1) &
          & + CpDry * z_HydroAlpha(k  ) * wz_siMtxDi(:,k  ) * wz_siTemp(:,k  ) &
          & + CpDry * z_HydroBeta (k-1) * wz_siMtxDi(:,k-1) * wz_siTemp(:,k-1)
      end do

      w_siPi = w_PiB + DelTime * w_DPiDtNG

      ! 逋ｺ謨｣譁ｹ遞句ｼ上�蜿ｳ霎ｺ縺ｮ險育ｮ�
      ! Calculate right side of divergence equation
      !
      do k = 1, kmax
        wz_siVectF(:,k) =                                                &
          &   ( 1.0_DP - DelTime * wz_DisCoefM (:,k) ) * wz_DivB(:,k)        &
!
          & + DelTime * wz_DDivDtNG(:,k)                                 &
!
          & - DelTime * w_LaplaEigVal(:)                                 &
          &   * ( w_SurfGeoPot + wz_siPhi(:,k) + z_siMtxG(k) * w_siPi )
      end do

      ! 譎る俣蟷ｳ蝮�� $ \Dvect{D} $ 繧� LU 陦悟�縺ｧ隗｣縺�
      ! Solve time average $ \Dvect{D} $ with LU matrix
      !
      wz_siDivAvrTime = LUSolve( wzz_siMtxLU, wz_siMtxPiv, wz_siVectF )

      wz_DivA = 2. * wz_siDivAvrTime - wz_DivB

    case ( IDTimeIntegSchemeExplicit )

      wz_WT = 0.0_DP
      do k = 1, kmax
        do kk = 1, kmax
          wz_WT(:,k) = wz_WT(:,k) + zz_siMtxW(k,kk) * wz_TempN(:,kk)
        end do
      end do

      do k = 1, kmax
        wz_DivA(:,k) =                                          &
          &   ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM(:,k) ) ) &
          & * ( wz_DivB(:,k)                                    &
          &     + 2.0_DP * DelTime                              &
          &       * (   wz_DDivDtNG(:,k)                        &
          &           - w_LaplaEigVal(:)                        &
          &               * (   w_SurfGeoPot                    &
          &                   + wz_WT(:,k)                      &
          &                   + z_siMtxG(k) * w_PiN             &
          &                 )                                   &
          &         )                                           &
          &   )
      end do

    end select

    ! 貂ｩ蠎ｦ ; Temperature
    !
    select case ( IDTimeIntegScheme )
    case ( IDTimeIntegSchemeSemiImplicit )
      wz_HDiv = 0.
      do k = 1, kmax
        do kk = 1, kmax
          wz_HDiv(:,k) = wz_HDiv(:,k) + zz_siMtxH(k,kk) * wz_siDivAvrTime(:,kk)
        end do
      end do
    case ( IDTimeIntegSchemeExplicit )
      wz_HDiv = 0.
      do k = 1, kmax
        do kk = 1, kmax
          wz_HDiv(:,k) = wz_HDiv(:,k) + zz_siMtxH(k,kk) * wz_DivN(:,kk)
        end do
      end do
    end select

    wz_TempA = &
      & ( 1.0_DP  / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefH ) ) &
      &   * ( wz_TempB + 2.0_DP * DelTime * ( wz_DTempDtNG - wz_HDiv ) )

    ! 蝨ｰ陦ｨ髱｢豌怜悸 ; Surface pressure
    !
    select case ( IDTimeIntegScheme )
    case ( IDTimeIntegSchemeSemiImplicit )
      w_CtDiv = 0.0_DP
      do k = 1, kmax
        w_CtDiv = w_CtDiv + z_siMtxC(k) * wz_siDivAvrTime(:,k)
      end do
    case ( IDTimeIntegSchemeExplicit )
      w_CtDiv = 0.0_DP
      do k = 1, kmax
        w_CtDiv = w_CtDiv + z_siMtxC(k) * wz_DivN(:,k)
      end do
    end select

    w_PiA  = w_PiB + 2.0_DP * DelTime * ( w_DPiDtNG - w_CtDiv )

  end subroutine TimeIntegration

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

  subroutine TimeIntegrationQMix( &
    & wzf_DQMixDtN,   & ! (in)
    & wzf_QMixB,      & ! (in)
    & wzf_QMixA       & ! (out)
    & )
    !
    ! 譎る俣遨榊�繧定｡後＞,
    ! 譎ょ綾 $ t $ 縺ｮ迚ｩ逅�㍼縺ｮ譎る俣螟牙喧縺ｨ $ t-\Delta t$ 縺ｮ迚ｩ逅�㍼縺九ｉ
    ! 譎ょ綾 $ t+\Delta t $ 縺ｮ迚ｩ逅�㍼繧定ｨ育ｮ励＠縺ｾ縺�.
    !
    ! 譎る俣遨榊�豕輔↓縺ｯ繝ｪ繝ｼ繝励ヵ繝ｭ繝�げ繧ｹ繧ｭ繝ｼ繝�繧堤畑縺�※縺�∪縺�.
    ! 縺溘□縺玲僑謨｣鬆�↓縺ｯ蠕梧婿蟾ｮ蛻�ｒ逕ｨ縺�※縺�∪縺�.
    ! 繝�ヵ繧ｩ繝ｫ繝医〒縺ｯ, $ \Delta t $ 繧貞､ｧ縺阪￥縺ｨ繧九◆繧√↓, 驥榊鴨豕｢鬆�↓
    ! 繧ｻ繝溘う繝ｳ繝励Μ繧ｷ繝�ヨ豕輔ｒ驕ｩ逕ｨ縺励※縺�∪縺�.
    ! NAMELIST#dynamics_hspl_vas83_nml 縺ｮ TimeIntegScheme 繧貞､画峩縺吶ｋ縺薙→縺ｧ,
    ! 驥榊鴨豕｢鬆�ｒ繧ｨ繧ｯ繧ｹ繝励Μ繧ｷ繝�ヨ豕輔↓繧医▲縺ｦ隗｣縺上％縺ｨ繧ょ庄閭ｽ縺ｧ縺�.
    !
    ! With time integration, physical values at $ t+\Delta t $ is calculated
    ! from tendency at $ t $ and physical values at $ t-\Delta t $ .
    !
    ! Leap-frog scheme is used as time integration scheme.
    ! And backward scheme is used for diffusion terms.
    ! As default setting, semi-implicit scheme is applied to gravitational terms
    ! in order to enlarge the value of $ \Delta t $ .
    ! Explicit scheme can be applied to gravitational terms by changing
    ! "TimeIntegScheme" in "NAMELIST#dynamics_hspl_vas83_nml".
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; USE statements
    !

    ! 譎ょ綾邂｡逅�
    ! Time control
    !
    use timeset, only: DelTime ! $ \Delta t $ [s]

    ! 譬ｼ蟄千せ險ｭ螳�
    ! Grid points settings
    !
    use gridset, only: lmax, & ! 繧ｹ繝壹け繝医Ν繝��繧ｿ縺ｮ驟榊�繧ｵ繧､繧ｺ
                               ! Size of array for spectral data
      &                imax, & ! 邨悟ｺｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in longitude
      &                jmax, & ! 邱ｯ蠎ｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in latitude
      &                kmax    ! 驩帷峩螻､謨ｰ.
                               ! Number of vertical level

    ! 邨��縺ｫ髢｢繧上ｋ驟榊�縺ｮ險ｭ螳�
    ! Settings of array for atmospheric composition
    !
    use composition, only: &
      &                    ncmax
                               ! 謌仙�縺ｮ謨ｰ
                               ! Number of composition

    implicit none

    real(DP), intent(in):: wzf_DQMixDtN(lmax, 1:kmax, 1:ncmax)
                              ! $ \DD{q}{t} (t) $ . 豈疲ｹｿ螟牙喧 (繧ｹ繝壹け繝医Ν).
                              ! Specific humidity tendency (spectral)
    real(DP), intent(in):: wzf_QMixB(lmax, 1:kmax, 1:ncmax)
                              ! $ q (t-\Delta t) $ . 豈疲ｹｿ (繧ｹ繝壹け繝医Ν).
                              ! Specific humidity (spectral)

    real(DP), intent(out):: wzf_QMixA(lmax, 1:kmax, 1:ncmax)
                              ! $ q (t+\Delta t) $ . 豈疲ｹｿ (繧ｹ繝壹け繝医Ν).
                              ! Specific humidity (spectral)

    ! 菴懈･ｭ螟画焚
    ! Work variables
    !
    integer:: k               ! 驩帷峩譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 邨��譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in dimension of constituents

    ! 螳溯｡梧枚 ; Executable statement
    !

    ! 譎る俣遨榊�. 諡｡謨｣縺ｯ蠕梧婿蟾ｮ蛻�
    ! Time integration. Backward difference is applied to diffusion
    !

    ! 豈疲ｹｿ ; Specific humidity
    !
    do n = 1, ncmax
      do k = 1, kmax
        wzf_QMixA(:,k,n) = &
          & ( 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * w_DisCoefQ ) ) &
          &   * ( wzf_QMixB(:,k,n) + 2.0_DP * DelTime * wzf_DQMixDtN(:,k,n) )
      end do
    end do

  end subroutine TimeIntegrationQMix

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

  subroutine DivergenceDamping( &
    & wz_DivA                   & ! (inout)
    & )
    !
    ! 逋ｺ謨｣縺ｮ貂幄｡ｰ鬆��莉伜刈 (蝣ｴ縺後▽繧雁粋縺｣縺ｦ縺�↑縺�ｨ育ｮ怜�譛溘�謫ｾ荵ｱ縺ｮ貂幄｡ｰ繧呈Φ螳�)
    !
    ! Addition of divergence damping term
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; USE statements
    !

    ! 譎ょ綾邂｡逅�
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $ [s]
      & TimeN                 ! 繧ｹ繝�ャ繝� $ t $ 縺ｮ譎ょ綾. Time of step $ t $.

    ! 譬ｼ蟄千せ險ｭ螳�
    ! Grid points settings
    !
    use gridset, only: lmax, & ! 繧ｹ繝壹け繝医Ν繝��繧ｿ縺ｮ驟榊�繧ｵ繧､繧ｺ
                               ! Size of array for spectral data
      &                imax, & ! 邨悟ｺｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in longitude
      &                jmax, & ! 邱ｯ蠎ｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in latitude
      &                kmax    ! 驩帷峩螻､謨ｰ.
                               ! Number of vertical level

    implicit none

    real(DP), intent(inout):: wz_DivA (lmax, 1:kmax)
                              ! $ D (t+\Delta t) $ . 逋ｺ謨｣ (繧ｹ繝壹け繝医Ν).
                              ! Divergence (spectral)

    ! 菴懈･ｭ螟画焚
    ! Work variables
    !
    real(DP) :: DivDampCoef

    ! 螳溯｡梧枚 ; Executable statement
    !

    if ( FlagDivDamp ) then

      DivDampCoef = 1.0_DP / DelTime * ( DivDampPeriod - ( TimeN - DivDampTime0 ) ) / DivDampPeriod
      if ( DivDampCoef < 0.0_DP ) DivDampCoef = 0.0_DP

      wz_DivA = wz_DivA / ( 1.0_DP + 2.0_DP * DelTime * DivDampCoef )

    end if

  end subroutine DivergenceDamping

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

  subroutine OutputDiagnosedVariables(                           &
    & xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB,             & ! (in)
    & xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN,             & ! (in)
    & xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA,             & ! (in)
    & xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy, & ! (in)
    & xyr_SigDot, xy_DPiDt,                                      & ! (in)
    & xyz_PiAdv,                                                 & ! (in)
    & xyz_OMG                                                    & ! (out)
    & )
    !
    ! 險ｺ譁ｭ驥上�蜃ｺ蜉帙ｒ陦後＞縺ｾ縺�.
    !
    ! 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: &
      & DelTime, &            ! $ \Delta t $ [s]
      & TimeN                 ! 繧ｹ繝�ャ繝� $ t $ 縺ｮ譎ょ綾. Time of step $ t $.

    ! 譬ｼ蟄千せ險ｭ螳�
    ! Grid points settings
    !
    use gridset, only: lmax, & ! 繧ｹ繝壹け繝医Ν繝��繧ｿ縺ｮ驟榊�繧ｵ繧､繧ｺ
                               ! Size of array for spectral data
      &                imax, & ! 邨悟ｺｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in longitude
      &                jmax, & ! 邱ｯ蠎ｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in latitude
      &                kmax    ! 驩帷峩螻､謨ｰ.
                               ! Number of vertical level

    ! 邨��縺ｫ髢｢繧上ｋ驟榊�縺ｮ險ｭ螳�
    ! Settings of array for atmospheric composition
    !
    use composition, only: &
      &                    ncmax, &
                               ! 謌仙�縺ｮ謨ｰ
                               ! Number of composition
      &                    a_QMixName, &
      &                    IndexH2OVap

    ! 蠎ｧ讓吶ョ繝ｼ繧ｿ險ｭ螳�
    ! 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)

    ! SPMODEL 繝ｩ繧､繝悶Λ繝ｪ, 逅�擇荳翫�蝠城｡後ｒ逅�擇隱ｿ蜥悟�謨ｰ螟画鋤縺ｫ繧医ｊ隗｣縺�(螟壼ｱ､蟇ｾ蠢�)
    ! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
    !
#ifdef LIB_MPI
#ifdef SJPACK
    use wa_mpi_module_sjpack, only:                                 &
      & w_xy => w_xv,                                               &
      & xy_GradLon_w => xv_GradLon_w, xy_GradLat_w => xv_GradLat_w, &
      & wa_DivLambda_xya => wa_DivLambda_xva,                       &
      & wa_DivMu_xya => wa_DivMu_xva,                               &
      & xya_wa => xva_wa
#else
    use wa_mpi_module, only:                                        &
      & w_xy => w_xv,                                               &
      & xy_GradLon_w => xv_GradLon_w, xy_GradLat_w => xv_GradLat_w, &
      & wa_DivLambda_xya => wa_DivLambda_xva,                       &
      & wa_DivMu_xya => wa_DivMu_xva,                               &
      & xya_wa => xva_wa
#endif
#elif AXISYMMETRY
    use wa_zonal_module, only:      &
      & w_xy,                       &
      & xy_GradLon_w, xy_GradLat_w, &
      & wa_DivLambda_xya,           &
      & wa_DivMu_xya,               &
      & xya_wa
#elif SJPACK
    use wa_module_sjpack, only:     &
      & w_xy,                       &
      & xy_GradLon_w, xy_GradLat_w, &
      & wa_DivLambda_xya,           &
      & wa_DivMu_xya,               &
      & xya_wa
#elif AXISYMMETRY_SJPACK
    use wa_zonal_module_sjpack, only: &
      & w_xy,                         &
      & xy_GradLon_w, xy_GradLat_w,   &
      & wa_DivLambda_xya,             &
      & wa_DivMu_xya,                 &
      & xya_wa
#else
    use wa_module, only:            &
      & w_xy,                       &
      & xy_GradLon_w, xy_GradLat_w, &
      & wa_DivLambda_xya,           &
      & wa_DivMu_xya,               &
      & xya_wa
#endif

    ! 繝偵せ繝医Μ繝��繧ｿ蜃ｺ蜉�
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 螳｣險譁� ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyz_UB    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u (t-\Delta t) $ .   譚ｱ隘ｿ鬚ｨ騾�. Eastward wind
    real(DP), intent(in):: xyz_VB    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v (t-\Delta t) $ .   蜊怜圏鬚ｨ騾�. Northward wind
    real(DP), intent(in):: xyz_TempB (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t-\Delta t) $ .   貂ｩ蠎ｦ. Temperature
    real(DP), intent(in):: xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t-\Delta t) $ .   豈疲ｹｿ. Specific humidity
    real(DP), intent(in):: xy_PsB    (0:imax-1, 1:jmax)
                              ! $ p_s (t-\Delta t) $ . 蝨ｰ陦ｨ髱｢豌怜悸. Surface pressure

    real(DP), intent(in):: xyz_UN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u (t) $ .     譚ｱ隘ｿ鬚ｨ騾�. Eastward wind
    real(DP), intent(in):: xyz_VN   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v (t) $ .     蜊怜圏鬚ｨ騾�. Northward wind
    real(DP), intent(in):: xyz_TempN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t) $ .     貂ｩ蠎ｦ. Temperature
    real(DP), intent(in):: xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t) $ .     豈疲ｹｿ. Specific humidity
    real(DP), intent(in):: xy_PsN    (0:imax-1, 1:jmax)
                              ! $ p_s (t) $ .   蝨ｰ陦ｨ髱｢豌怜悸. Surface pressure

    real(DP), intent(in):: xyz_UA    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u (t+\Delta t) $ .   譚ｱ隘ｿ鬚ｨ騾�. Eastward wind
    real(DP), intent(in):: xyz_VA    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v (t+\Delta t) $ .   蜊怜圏鬚ｨ騾�. Northward wind
    real(DP), intent(in):: xyz_TempA (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t+\Delta t) $ .   貂ｩ蠎ｦ. Temperature
    real(DP), intent(in):: xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t+\Delta t) $ .   豈疲ｹｿ. Specific humidity
    real(DP), intent(in):: xy_PsA    (0:imax-1, 1:jmax)
                              ! $ p_s (t+\Delta t) $ . 蝨ｰ陦ｨ髱｢豌怜悸. Surface pressure

    real(DP), intent(in):: xyz_DUDtPhy    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \left(\DP{u}{t}\right)^{phy} $ .
                              ! 螟門鴨鬆� (迚ｩ逅�℃遞�) 縺ｫ繧医ｋ譚ｱ隘ｿ鬚ｨ騾溷､牙喧.
                              ! Eastward wind tendency by external force terms (physical processes)
    real(DP), intent(in):: xyz_DVDtPhy    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \left(\DP{v}{t}\right)^{phy} $ .
                              ! 螟門鴨鬆� (迚ｩ逅�℃遞�) 縺ｫ繧医ｋ蜊怜圏鬚ｨ騾溷､牙喧.
                              ! Northward wind tendency by external force terms (physical processes)
    real(DP), intent(in):: xyz_DTempDtPhy (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \left(\DP{T}{t}\right)^{phy} $ .
                              ! 螟門鴨鬆� (迚ｩ逅�℃遞�) 縺ｫ繧医ｋ貂ｩ蠎ｦ螟牙喧.
                              ! Temperature tendency by external force terms (physical processes)
    real(DP), intent(in):: xyzf_DQMixDtPhy(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \left(\DP{q}{t}\right)^{phy} $ .
                              ! 螟門鴨鬆� (迚ｩ逅�℃遞�) 縺ｫ繧医ｋ豈疲ｹｿ螟牙喧.
                              ! Temperature tendency by external force terms (physical processes)
    real(DP), intent(in):: xyr_SigDot (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \dot{\sigma} $ .
                              ! 驩帷峩豬�. Vertical flow
    real(DP), intent(in):: xy_DPiDt (0:imax-1, 1:jmax)
                              ! $ Z $ . 蝨ｰ陦ｨ髱｢豌怜悸譎る俣螟牙喧鬆�.
                              ! Surface pressure tendency

    real(DP), intent(in):: xyz_PiAdv (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \Dvect{v} \cdot \nabla \pi $ .
                              ! $ \pi $ 縺ｮ遘ｻ豬�. Advection of $ \pi $
    real(DP), intent(out):: xyz_OMG  (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! OMEGA, DP/Dt

    ! 菴懈･ｭ螟画焚
    ! Work variables
    !
    real(DP):: xyz_DUDtDyn    (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_DVDtDyn    (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_DTempDtDyn (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyzf_DQMixDtDyn(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP):: xy_DPsDtDyn    (0:imax-1, 1:jmax)

    real(DP):: xyz_UCosLat   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ U = u \cos \varphi $ .
    real(DP):: xyz_VCosLat   (0:imax-1, 1:jmax, 1:kmax)
                              ! $ V = v \cos \varphi $ .
    real(DP):: xyz_Vor       (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \zeta $ . 貂ｦ蠎ｦ. Vorticity
    real(DP):: xyz_Div       (0:imax-1, 1:jmax, 1:kmax)
                              ! $ D $ . 逋ｺ謨｣. Divergence

    real(DP):: xyr_PiAdvDivSum(0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xy_Mass (0:imax-1, 1:jmax)
                              ! 雉ｪ驥�.
                              ! Mass
    real(DP):: xyz_KinEngy (0:imax-1, 1:jmax, 1:kmax)
                              ! $ KE $ . 驕句虚繧ｨ繝阪Ν繧ｮ繝ｼ.
                              ! Kinetic energy
    real(DP):: xyz_IntEngy (0:imax-1, 1:jmax, 1:kmax)
                              ! $ IE $ . 蜀�Κ繧ｨ繝阪Ν繧ｮ繝ｼ.
                              ! Internal energy
    real(DP):: xyz_PotEngy (0:imax-1, 1:jmax, 1:kmax)
                              ! $ PE $ . 繝昴ユ繝ｳ繧ｷ繝｣繝ｫ繧ｨ繝阪Ν繧ｮ繝ｼ.
                              ! Potential energy
    real(DP):: xyz_LatEngy (0:imax-1, 1:jmax, 1:kmax)
                              ! $ LE $ . 貎懃�繧ｨ繝阪Ν繧ｮ繝ｼ.
                              ! Latent heat energy
    real(DP):: xyz_TotEngy (0:imax-1, 1:jmax, 1:kmax)
                              ! $ TE $ . 蜈ｨ繧ｨ繝阪Ν繧ｮ繝ｼ.
                              ! Total energy
    real(DP):: xyz_Enstro (0:imax-1, 1:jmax, 1:kmax)
                              ! 繧ｨ繝ｳ繧ｹ繝医Ο繝輔ぅ繝ｼ.
                              ! Enstrophy

    integer:: k               ! 驩帷峩譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in vertical direction
    integer:: n

    ! 螳溯｡梧枚 ; Executable statement
    !

    ! Calculate tendencies
    !
    xyz_DUDtDyn     = ( xyz_UA     - xyz_UB     ) / ( 2.0_DP * DelTime ) &
      &               - xyz_DUDtPhy
    xyz_DVDtDyn     = ( xyz_VA     - xyz_VB     ) / ( 2.0_DP * DelTime ) &
      &               - xyz_DVDtPhy
    xyz_DTempDtDyn  = ( xyz_TempA  - xyz_TempB  ) / ( 2.0_DP * DelTime ) &
      &               - xyz_DTempDtPhy
    xyzf_DQMixDtDyn = ( xyzf_QMixA - xyzf_QMixB ) / ( 2.0_DP * DelTime ) &
      &               - xyzf_DQMixDtPhy
    xy_DPsDtDyn     = ( xy_PsA     - xy_PsB     ) / ( 2.0_DP * DelTime )

    call HistoryAutoPut( TimeN, 'DUDtDyn'   , xyz_DUDtDyn                        )
    call HistoryAutoPut( TimeN, 'DVDtDyn'   , xyz_DVDtDyn                        )
    call HistoryAutoPut( TimeN, 'DTempDtDyn', xyz_DTempDtDyn                     )
    do n = 1, ncmax
      call HistoryAutoPut( TimeN, 'D'//trim(a_QMixName(n))//'DtDyn', xyzf_DQMixDtDyn(:,:,:,n) )
    end do
    call HistoryAutoPut( TimeN, 'DPsDtDyn'  , xy_DPsDtDyn                        )

    ! 驩帷峩豬√→蝨ｰ陦ｨ髱｢豌怜悸譎る俣螟牙喧鬆��蜃ｺ蜉�
    ! Output vertical flow and surface pressure tendency
    !
    call HistoryAutoPut( TimeN, 'SigDot', xyr_SigDot )
    call HistoryAutoPut( TimeN, 'DPiDt',  xy_DPiDt )

    ! 鬚ｨ騾溘°繧画ｸｦ蠎ｦ逋ｺ謨｣縺ｮ險育ｮ�
    ! Calculate vorticity and divergence from wind velocity
    !
    do k = 1, kmax
      xyz_UCosLat(:,:,k) = xyz_UN(:,:,k) * xy_CosLat
      xyz_VCosLat(:,:,k) = xyz_VN(:,:,k) * xy_CosLat
    end do
    xyz_Vor = xya_wa(   wa_DivLambda_xya( xyz_VCosLat )  &
      &               - wa_DivMu_xya(     xyz_UCosLat ) ) / RPlanet
    xyz_Div = xya_wa(   wa_DivLambda_xya( xyz_UCosLat )  &
      &               + wa_DivMu_xya(     xyz_VCosLat ) ) / RPlanet

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

    ! Calculation of Omega (DPressDt)
    ! It should be recognized that the value of xy_Ps here may have been modified
    ! by mass fixer.
    !
    !   Integration from top of the model to k's layer upper interface
    xyr_PiAdvDivSum(:,:,kmax) = 0.0_DP
    do k = kmax-1, 0, -1
      xyr_PiAdvDivSum(:,:,k) = xyr_PiAdvDivSum(:,:,k+1) &
        & + ( xyz_PiAdv(:,:,k+1) + xyz_Div(:,:,k+1) ) * z_DelSigma(k+1)
    end do
    do k = 1, kmax
      xyz_OMG(:,:,k) =                                                                  &
        & xy_PsN * (                                                                    &
        &           z_Sigma(k) * xyz_PiAdv(:,:,k)                                       &
        &         - xyr_PiAdvDivSum(:,:,k)                                              &
        &         - ( xyz_PiAdv(:,:,k) + xyz_Div(:,:,k) ) * ( z_Sigma(k) - r_Sigma(k) ) &
        &          )
    end do
    call HistoryAutoPut( TimeN, 'OMG', xyz_OMG )

    ! 雉ｪ驥上�險育ｮ�
    ! Calculate mass
    !
    xy_Mass = xy_PsN / Grav

    ! 繧ｨ繝阪Ν繧ｮ繝ｼ, 繧ｨ繝ｳ繧ｹ繝医Ο繝輔ぅ繝ｼ縺ｮ險育ｮ�
    ! Calculate energy and enstrophy
    !
    call HydroGrid( xyz_TempN,   & ! (in)
      &             xyz_PotEngy  & ! (out)
      & )

    do k = 1, kmax
      xyz_KinEngy(:,:,k) = ( xyz_UN(:,:,k) ** 2 + xyz_VN(:,:,k) ** 2 ) / 2.0_DP &
        &                    * xy_Mass

      xyz_IntEngy(:,:,k) = CpDry * xyz_TempN(:,:,k) &
        &                    * xy_Mass

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

      xyz_LatEngy(:,:,k) = LatentHeat * xyzf_QMixN(:,:,k,IndexH2OVap) &
        &                    * xy_Mass
    end do

    xyz_TotEngy = xyz_KinEngy + xyz_IntEngy + xyz_PotEngy + xyz_LatEngy

    do k = 1, kmax
      xyz_Enstro(:,:,k) = xyz_Vor(:,:,k) ** 2 &
        &                   * xy_Mass
    end do

    call HistoryAutoPut( TimeN, 'Mass',    xy_Mass     )
    call HistoryAutoPut( TimeN, 'KinEngy', xyz_KinEngy )
    call HistoryAutoPut( TimeN, 'IntEngy', xyz_IntEngy )
    call HistoryAutoPut( TimeN, 'PotEngy', xyz_PotEngy )
    call HistoryAutoPut( TimeN, 'LatEngy', xyz_LatEngy )
    call HistoryAutoPut( TimeN, 'TotEngy', xyz_TotEngy )
    call HistoryAutoPut( TimeN, 'Enstro',  xyz_Enstro  )

  end subroutine OutputDiagnosedVariables

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

  subroutine DynamicsHSplVAS83Init
    !
    ! 險育ｮ励↓蠢�ｦ√↑繝代Λ繝｡繧ｿ縺ｮ險ｭ螳壹ｄ NAMELIST#dynamics_hspl_vas83_nml
    ! 縺ｮ隱ｭ縺ｿ霎ｼ縺ｿ繧定｡後＞縺ｾ縺�.
    !
    ! Configure parameters for calculation,
    ! and load "NAMELIST#dynamics_hspl_vas83_nml"
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; USE statements
    !

    ! 迚ｩ逅�ｮ壽焚險ｭ螳�
    ! Physical constants settings
    !
    use constants, only: &
      & RPlanet, &
                              ! $ 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

    ! 譬ｼ蟄千せ險ｭ螳�
    ! Grid points settings
    !
    use gridset, only: nmax, & ! 譛螟ｧ蜈ｨ豕｢謨ｰ.
                               ! Maximum truncated wavenumber
      &                lmax, & ! 繧ｹ繝壹け繝医Ν繝��繧ｿ縺ｮ驟榊�繧ｵ繧､繧ｺ
                               ! Size of array for spectral data
      &                imax, & ! 邨悟ｺｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in longitude
      &                jmax, & ! 邱ｯ蠎ｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in latitude
      &                kmax    ! 驩帷峩螻､謨ｰ.
                               ! Number of vertical level

    ! 邨��縺ｫ髢｢繧上ｋ驟榊�縺ｮ險ｭ螳�
    ! Settings of array for atmospheric composition
    !
    use composition, only: &
      &                    ncmax, &
                               ! 謌仙�縺ｮ謨ｰ
                               ! Number of composition
      &                    a_QMixName

    ! 蠎ｧ讓吶ョ繝ｼ繧ｿ險ｭ螳�
    ! 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)
      & AxNameX, AxNameY, AxNameZ, AxNameR, AxNameT

    ! SPMODEL 繝ｩ繧､繝悶Λ繝ｪ, 逅�擇荳翫�蝠城｡後ｒ逅�擇隱ｿ蜥悟�謨ｰ螟画鋤縺ｫ繧医ｊ隗｣縺�(螟壼ｱ､蟇ｾ蠢�)
    ! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
    !
#ifdef LIB_MPI
#ifdef SJPACK
    use wa_mpi_module_sjpack, only: &
      & xy_Lat => xv_Lat,    &
      & w_xy   => w_xv,      &
      & rn,                  &
      & nm_l
#else
    use wa_mpi_module, only: &
      & xy_Lat => xv_Lat,    &
      & w_xy   => w_xv,      &
      & rn,                  &
      & nm_l
#endif
#elif AXISYMMETRY
    use wa_zonal_module, only: &
      & xy_Lat,          &
      & w_xy,            &
      & rn,              &
      & nm_l
#elif SJPACK
    use wa_module_sjpack, only: &
      & xy_Lat,          &
      & w_xy,            &
      & rn,              &
      & nm_l
#elif AXISYMMETRY_SJPACK
    use wa_zonal_module_sjpack, only: &
      & xy_Lat,          &
      & w_xy,            &
      & rn,              &
      & nm_l
#else
    use wa_module, only: &
      & xy_Lat,          &
      & w_xy,            &
      & rn,              &
      & nm_l
#endif

    ! NAMELIST 繝輔ぃ繧､繝ｫ蜈･蜉帙↓髢｢縺吶ｋ繝ｦ繝ｼ繝�ぅ繝ｪ繝�ぅ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg

    ! 繝偵せ繝医Μ繝��繧ｿ蜃ｺ蜉�
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! gtool4 繝��繧ｿ蜈･蜉�
    ! Gtool4 data input
    !
    use gtool_history, only: HistoryGet

    ! 繝輔ぃ繧､繝ｫ蜈･蜃ｺ蜉幄｣懷勧
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 遞ｮ蛻･蝙九ヱ繝ｩ繝｡繧ｿ
    ! Kind type parameter
    !
    use dc_types, only: &
      & STDOUT, &             ! 讓呎ｺ門�蜉帙�陬�ｽｮ逡ｪ蜿ｷ. Unit number of standard output
      & STRING                ! 譁�ｭ怜�.       Strings.

    ! 譌･莉倥♀繧医�譎ょ綾縺ｮ蜿悶ｊ謇ｱ縺�
    ! Date and time handler
    !
    use dc_calendar, only: DCCalConvertByUnit

    ! 邨�∩霎ｼ縺ｿ髢｢謨ｰ PRESENT 縺ｮ諡｡蠑ｵ迚磯未謨ｰ
    ! Extended functions of intrinsic function "PRESENT"
    !
    use dc_present, only: present_and_not_empty

    ! 譁�ｭ怜�謫堺ｽ�
    ! Character handling
    !
    use dc_string, only: LChar

    ! 譎ょ綾邂｡逅�
    ! Time control
    !
    use timeset, only: &
      & TimeN                 ! 繧ｹ繝�ャ繝� $ t $ 縺ｮ譎ょ綾. Time of step $ t $.

    ! 雉ｪ驥上�陬懈ｭ｣
    ! Mass fixer
    !
    use mass_fixer, only : MassFixerInit

    ! 繧ｻ繝溘Λ繧ｰ繝ｩ繝ｳ繧ｸ繝･豕輔↓繧医ｋ迚ｩ雉ｪ遘ｻ豬∬ｨ育ｮ�
    ! Semi-Lagrangian method for tracer transport
    use sltt, only: SLTTInit

!!$    ! 迚ｩ雉ｪ遘ｻ豬� (譛蛾剞菴鍋ｩ肴ｳ�)
!!$    ! Tracer Transport (finite volume method)
!!$    !
!!$    use tt_fv, only: TTFVInit

    !
    ! Utility module for advection test
    !
    use adv_test, only: AdvTestInit

    ! 螳｣險譁� ; Declaration statements
    !
    implicit none

    ! 蝓ｺ貅匁ｸｩ蠎ｦ縺ｮ險ｭ螳壹�縺溘ａ縺ｮ菴懈･ｭ螟画焚
    ! Work variable for reference temperature
    !
    character(TOKEN) :: TimeIntegScheme
                              ! 譎る俣遨榊�豕�.
                              ! 莉･荳九�譁ｹ豕輔ｒ驕ｸ謚槫庄閭ｽ.
                              !
                              ! Time integration scheme.
                              ! Available schemes are as follows.
                              !
                              ! * "Semi-implicit"
                              ! * "Explicit"

    real(DP):: RefTemp
                              ! $ \overline{T} $ . 蝓ｺ貅匁ｸｩ蠎ｦ.
                              ! Reference temperature

    ! 豌ｴ蟷ｳ諡｡謨｣, 繧ｹ繝昴Φ繧ｸ螻､縺ｮ縺溘ａ縺ｮ菴懈･ｭ螟画焚
    ! Work variable for horizontal diffusion and sponge layer
    !
    real(DP) :: w_HDifCoefM        (lmax)
                              ! $ {\cal D}_M = -K_{HD} [(-1)^{N_D/2}\nabla^{N_D}- (2/a^2)^{N_D/2}] $ .
                              ! 驕句虚驥乗ｰｴ蟷ｳ諡｡謨｣菫よ焚.
                              ! Coefficient of horizontal momentum diffusion
    real(DP) :: w_HDifCoefH        (lmax)
                              ! $ {\cal D}_H = -(-1)^{N_D/2}K_{HD}\nabla^{N_D} $ .
                              ! 辭ｱ, 豌ｴ豌ｴ蟷ｳ諡｡謨｣菫よ焚.
                              ! Coefficient of horizontal thermal and water diffusion
    real(DP) :: wz_SpongeLayerCoefM(lmax, 1:kmax)
                              ! 繧ｹ繝昴Φ繧ｸ螻､縺ｫ縺翫￠繧矩°蜍暮㍼縺ｮ貂幄｡ｰ菫よ焚
                              ! Damping coefficient for momentum in a sponge layer
    real(DP) :: wz_SpongeLayerCoefH(lmax, 1:kmax)
                              ! 繧ｹ繝昴Φ繧ｸ螻､縺ｫ縺翫￠繧狗�縺ｮ貂幄｡ｰ菫よ焚
                              ! Damping coefficient for temperature zonal perturbation
                              ! in a sponge layer

    integer:: HDOrder           ! 雜�ｲ俶ｧ縺ｮ谺｡謨ｰ.  Order of hyper-viscosity
    real(DP):: VisCoef          ! 雜�ｲ俶ｧ菫よ焚. Hyper-viscosity coefficient
    real(DP):: HDEFoldTimeValue ! 譛螟ｧ豕｢謨ｰ縺ｫ蟇ｾ縺吶ｋ e-folding time.
                                ! 雋�縺ｮ蛟､繧剃ｸ弱∴繧九→, 豌ｴ蟷ｳ諡｡謨｣菫よ焚繧偵ぞ繝ｭ縺ｫ縺励∪縺�.
                                !
                                ! E-folding time for maximum wavenumber.
                                ! If negative value is given,
                                ! coefficients of horizontal diffusion become zero.
    character(TOKEN):: HDEFoldTimeUnit
                                ! 譛螟ｧ豕｢謨ｰ縺ｫ蟇ｾ縺吶ｋ e-folding time 縺ｮ蜊倅ｽ�.
                                ! Unit of e-folding time for maximum wavenumber
    real(DP):: HDEFoldTime      ! 譛螟ｧ豕｢謨ｰ縺ｫ蟇ｾ縺吶ｋ e-folding time [蜊倅ｽ�: 遘綻.
                                ! E-folding time for maximum wavenumber [Unit: sec]

    logical          :: FlagSpongeLayer
    logical          :: FlagSpongeLayerforZonalMean
    logical          :: FlagSpongeLayerforHeat
    real(DP)         :: SLEFoldTimeValue
                                ! 繧ｹ繝昴Φ繧ｸ螻､縺ｮ譛荳雁ｱ､縺ｫ縺翫￠繧区ｸ幄｡ｰ譎ょｮ壽焚
                                ! Damping time constant at the top model layer
                                ! in a sponge layer
    character(TOKEN) :: SLEFoldTimeUnit
                                ! SLEFoldTimeValue 縺ｮ蜊倅ｽ�
                                ! Unit of SLEFoldTimeValue
    real(DP)         :: SLEFoldTime
                                ! 繧ｹ繝昴Φ繧ｸ螻､縺ｮ譛荳雁ｱ､縺ｫ縺翫￠繧区ｸ幄｡ｰ譎ょｮ壽焚 (遘�)
                                ! Damping time constant at the top model layer
                                ! in a sponge layer (sec)
    integer          :: SLOrder
                                ! 繧ｹ繝昴Φ繧ｸ螻､縺ｮ貂幄｡ｰ菫よ焚縺ｮ sigma 萓晏ｭ俶ｧ縺ｮ繧ｪ繝ｼ繝
                                ! Sigma dependence of damping coefficient
                                ! in a sponge layer
    integer          :: SLNumLayer
                                ! 繧ｹ繝昴Φ繧ｸ螻､縺碁←蠢懊＆繧後ｋ繝｢繝�Ν荳顔ｫｯ縺九ｉ縺ｮ螻､縺ｮ謨ｰ
                                ! Number of layers which the sponge layer is applied.
    integer          :: a_DegOrd(2)

    real(DP)         :: DivDampPeriodValue
                                ! Period for divergence damping application
    character(TOKEN) :: DivDampPeriodUnit
                                ! Unit of DivDampPeriodValue

    character(TOKEN):: TracerTransportMethod

    ! NonLinearOnGrid 遲峨〒菴ｿ逕ｨ縺吶ｋ菫よ焚縺ｮ險ｭ螳壹�縺溘ａ縺ｮ菴懈･ｭ螟画焚
    ! Work variable for coefficients for "NonLinearOnGrid", etc.
    !
    real(DP):: Kappa          ! $ \kappa = R/C_p $ .
                              ! 豌嶺ｽ灘ｮ壽焚縺ｮ螳壼悸豈皮�縺ｫ蟇ｾ縺吶ｋ豈�.
                              ! Ratio of gas constant to specific heat

    integer:: k               ! 驩帷峩譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in vertical direction
    integer:: n
    integer:: l               ! 豕｢謨ｰ譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in wavenumber

    integer:: unit_nml        ! NAMELIST 繝輔ぃ繧､繝ｫ繧ｪ繝ｼ繝励Φ逕ｨ陬�ｽｮ逡ｪ蜿ｷ.
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 隱ｭ縺ｿ霎ｼ縺ｿ譎ゅ� IOSTAT.
                              ! IOSTAT of NAMELIST read

    ! NAMELIST 螟画焚鄒､
    ! NAMELIST group name
    !
    namelist /dynamics_hspl_vas83_nml/                                        &
      & TimeIntegScheme,                                                      &
      & HDOrder,                                                              &
      & HDEFoldTimeValue, HDEFoldTimeUnit,                                    &
      & FlagSpongeLayer, FlagSpongeLayerforZonalMean, FlagSpongeLayerforHeat, &
      & SLEFoldTimeValue, SLEFoldTimeUnit, SLOrder, SLNumLayer,               &
      & RefTemp,                                                              &
      & FlagDivDamp, DivDampPeriodValue, DivDampPeriodUnit,                   &
      & DivDampTime0,                                                         &
      & FlagTotMassFix,                                                       &
      & FlagMassHorDifCor,                                                    &
!!$      & TracerTransportMethod,                                                &
      & FlagSLTT,                                                             &
      & FlagAdvTest
          !
          ! 繝�ヵ繧ｩ繝ｫ繝亥､縺ｫ縺､縺�※縺ｯ蛻晄悄蛹匁焔邯� "dynamics_hspl_vas83#DynamicsInit"
          ! 縺ｮ繧ｽ繝ｼ繧ｹ繧ｳ繝ｼ繝峨ｒ蜿ら�縺ｮ縺薙→.
          !
          ! Refer to source codes in the initialization procedure
          ! "dynamics_hspl_vas83#DynamicsInit" for the default values.
          !

    ! 螳溯｡梧枚 ; Executable statement
    !

    if ( dynamics_hspl_vas83_inited ) return

    ! 繝�ヵ繧ｩ繝ｫ繝亥､縺ｮ險ｭ螳�
    ! Default values settings
    !
    TimeIntegScheme             = 'Semi-implicit'

    HDOrder                     = 8
    HDEFoldTimeValue            = 8640.0_DP
    HDEFoldTimeUnit             = 'sec'

    FlagSpongeLayer             = .false.
    FlagSpongeLayerforZonalMean = .false.
    FlagSpongeLayerforHeat      = .false.
    SLEFoldTimeValue            = 86400.0_DP
    SLEFoldTimeUnit             = 'sec'
    SLOrder                     = 1
    SLNumLayer                  = kmax

    RefTemp                     = 300.

    FlagDivDamp                 = .false.
    DivDampPeriodValue          = 2.0_DP
    DivDampPeriodUnit           = 'day'
    DivDampTime0                = 0.0_DP

    FlagTotMassFix              = .true.

    FlagMassHorDifCor           = .false.

!!$    TracerTransportMethod       = "spectrum"
    FlagSLTT                    = .false.
    FlagAdvTest                 = .false.

    ! NAMELIST 縺ｮ隱ｭ縺ｿ霎ｼ縺ｿ
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, &          ! (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')
      IDTimeIntegScheme = IDTimeIntegSchemeSemiImplicit
    case ('explicit')
      IDTimeIntegScheme = IDTimeIntegSchemeExplicit
    case default
      call MessageNotify( 'E', module_name, &
        & '"TimeIntegScheme" must be "Semi-Implicit" or "Explicit".' )
    end select

    !
    ! Check tracer transport method
    !
!!$    select case ( LChar( trim( TracerTransportMethod ) ) )
!!$    case ('spectrum')
!!$      IDTTMethod = IDTTMethodSp
!!$    case ('semi-lagrange')
!!$      IDTTMethod = IDTTMethodSL
!!$    case ('finite volume')
!!$      IDTTMethod = IDTTMethodFV
!!$    case default
!!$      call MessageNotify( 'E', module_name, &
!!$        & '"TracerTransportMethod" must be "spectrum", "semi-lagrange", or "finite volume".' )
!!$    end select

    ! SemiImplMatrix 繧ｵ繝悶Ν繝ｼ繝√Φ逕ｨ縺ｫ $ \Delta t $ 縺ｮ菫晏ｭ伜､縺ｫ蛻晄悄蛟､繧定ｨｭ螳�.
    ! Configure initial value to saved value of $ \Delta t $ for a subroutine "SemiImplMatrix"
    !
    DelTimeSave = - 1.0_DP

    ! NonLinearOnGrid 遲峨〒菴ｿ逕ｨ縺吶ｋ菫よ焚縺ｮ險ｭ螳�
    ! Configure coefficients for "NonLinearOnGrid", etc.
    !

    ! $ \sin \varphi $ 縺ｨ $ \cos \varphi $ 縺ｮ險育ｮ�
    ! Calculate $ \sin \varphi $ and $ \cos \varphi $
    !
    allocate( xy_SinLat (0:imax-1, 1:jmax) )
    allocate( xy_CosLat (0:imax-1, 1:jmax) )
    xy_SinLat = sin( xy_Lat )
    xy_CosLat = cos( xy_Lat )

    ! 繧ｳ繝ｪ繧ｪ繝ｪ繝代Λ繝｡繝ｼ繧ｿ縺ｮ險育ｮ�
    ! Calculate Coriolis parameter
    !
    allocate( xy_Cori (0:imax-1, 1:jmax) )
    xy_Cori = 2.0_DP * Omega * xy_SinLat

    ! 髱呎ｰｴ蝨ｧ縺ｮ蠑上�菫よ焚 $ \alpha $ , $ \beta $ 縺ｮ險育ｮ�
    ! Calculate coefficients $ \alpha $ and $ \beta $ in hydrostatic equation.
    !
    allocate( z_HydroAlpha(1:kmax) )
    allocate( z_HydroBeta (1:kmax) )

    Kappa = GasRDry / CpDry

    do k = 1, kmax
      z_HydroAlpha(k) = &
        & ( r_Sigma(k-1) / z_Sigma(k) ) ** Kappa - 1.0_DP

      z_HydroBeta(k) = &
        & 1.0_DP - ( r_Sigma(k) / z_Sigma(k) ) ** Kappa
    enddo

    ! 貂ｩ蠎ｦ驩帷峩陬憺俣縺ｮ菫よ焚 $ \kappa $, $ a $ , $ b $ 縺ｮ險育ｮ�
    ! Calculate coefficients $ \kappa $, $ a $ , $ b $
    !   for interpolation of temperature
    !
    allocate( z_TInpCoefA (1:kmax) )
    allocate( z_TInpCoefB (1:kmax) )
    allocate( z_TInpCoefK (1:kmax) )

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

    z_TInpCoefA = 0.
    do k = 2, kmax
      z_TInpCoefA(k) = &
        & z_HydroAlpha(k) &
        &   * ( 1.0_DP - (z_Sigma(k) / z_Sigma(k-1)) ** Kappa ) ** ( -1 )
    end do

    z_TInpCoefB = 0.
    do k = 1, kmax - 1
      z_TInpCoefB(k) = &
        & z_HydroBeta(k) &
        &   * ( (z_Sigma(k) / z_Sigma(k+1) ) ** Kappa - 1.0_DP ) ** ( -1 )
    end do

    ! 蝓ｺ貅匁ｸｩ蠎ｦ (蜊頑紛謨ｰ繝ｬ繝吶Ν) 縺ｮ險育ｮ�
    ! Calculate reference temperature on half levels
    !
    allocate( z_RefTemp(1:kmax) )
    allocate( r_RefTemp(0:kmax) )

    z_RefTemp       = RefTemp
    r_RefTemp(0)    = 0.
    r_RefTemp(kmax) = 0.

    do k = 1, kmax - 1
      r_RefTemp(k) = &
        &   z_TInpCoefA(k+1) * z_RefTemp(k+1)   &
        & + z_TInpCoefB( k ) * z_RefTemp( k )
    enddo

    ! 繝ｫ繧ｸ繝｣繝ｳ繝峨Ν髯ｪ髢｢謨ｰ縺ｮ蝗ｺ譛牙､縺ｮ險ｭ螳�
    ! Set eigen value of Associated Legendre Function
    !
    allocate( w_LaplaEigVal(lmax) )

    w_LaplaEigVal(:) = rn(:,1) / RPlanet**2

    ! 豌ｴ蟷ｳ諡｡謨｣菫よ焚縺ｮ險ｭ螳�
    ! Configure coefficient of horizontal diffusion
    !
    ! 邊俶ｧ菫よ焚縺ｮ險育ｮ� (譛螟ｧ豕｢謨ｰ縺ｮ e-folding time 縺� HDEFoldTime 縺ｨ縺ｪ繧九ｈ縺�↓)
    ! Calculate viscosity coefficient
    !
    HDEFoldTime = DCCalConvertByUnit( HDEFoldTimeValue, HDEFoldTimeUnit, 'sec' ) ! (in)

    if ( HDEFoldTimeValue > 0.0_DP ) then
      VisCoef = ( (nmax*(nmax+1)) / RPlanet**2 )**(-HDOrder / 2) &
        &        / HDEFoldTime
    else
      VisCoef = 0.0_DP
    end if

    w_HDifCoefH = - VisCoef * ( ( - w_LaplaEigVal )**( HDOrder / 2 ) )

    w_HDifCoefM = w_HDifCoefH - VisCoef * ( - ( 2.0_DP / RPlanet**2 )**( HDOrder / 2 ) )

    ! 繧ｹ繝昴Φ繧ｸ螻､縺ｮ貂幄｡ｰ菫よ焚縺ｮ險ｭ螳�
    ! Set damping coefficient in a sponge layer
    !
    if ( SLEFoldTimeValue <= 0.0_DP ) then
      call MessageNotify( 'E', module_name, &
        & 'SLEFoldTimeValue must be greater than zero, SLEFoldTimeValue=<%f>.', &
        & d = (/ SLEFoldTimeValue /) )
    end if
    if ( ( SLNumLayer <= 0 ) .or. ( SLNumLayer > kmax ) ) then
      call MessageNotify( 'E', module_name, &
        & 'SLNumLayer must be greater than zero and less than/equal to kmax, SLNumLayer=<%d>.', &
        & i = (/ SLNumLayer /) )
    end if

    ! 繧ｹ繝昴Φ繧ｸ螻､縺ｮ貂幄｡ｰ菫よ焚縺ｮ險育ｮ�
    ! Calculate damping coefficient for momentum in a sponge layer
    !
    SLEFoldTime = DCCalConvertByUnit( SLEFoldTimeValue, SLEFoldTimeUnit, 'sec' ) ! (in)

    if ( FlagSpongeLayer ) then

      ! Sponge layer is not applied for model layers, k = 1, kmax-SLNumLayers.
      !
      do k = 1, kmax-SLNumLayer
        wz_SpongeLayerCoefM(:,k) = 0.0_DP
      end do
      do k = kmax-SLNumLayer+1, kmax
        wz_SpongeLayerCoefM(:,k) = &
          & 1.0_DP / ( SLEFoldTime * ( z_Sigma(k) / z_Sigma(kmax) )**SLOrder )
      end do

      if ( .not. FlagSpongeLayerforZonalMean ) then
        ! Sponge layer is not applied for zonal mean component.
        ! i.e., sponge layer is applied for zonal wave component only.
        !
        do k = kmax-SLNumLayer+1, kmax
          do l = 1, lmax
            a_DegOrd = nm_l( l )
            if ( a_DegOrd(2) == 0 ) then
              wz_SpongeLayerCoefM(l,k) = 0.0_DP
            end if
          end do
        end do
      end if

      if ( FlagSpongeLayerforHeat ) then

        wz_SpongeLayerCoefH = wz_SpongeLayerCoefM

        ! Sponge layer is not applied for zonal mean component.
        ! i.e., sponge layer is applied for zonal wave component only.
        !
        do k = kmax-SLNumLayer+1, kmax
          do l = 1, lmax
            a_DegOrd = nm_l( l )
            if ( a_DegOrd(2) == 0 ) then
              wz_SpongeLayerCoefH(l,k) = 0.0_DP
            end if
          end do
        end do
      else
        wz_SpongeLayerCoefH = 0.0_DP
      end if

    else
      ! Sponge layer is not applied.
      !
      wz_SpongeLayerCoefM = 0.0_DP
      wz_SpongeLayerCoefH = 0.0_DP
    end if

    ! 驕句虚驥上�豌ｴ蟷ｳ諡｡謨｣菫よ焚縺ｨ繧ｹ繝昴Φ繧ｸ螻､貂幄｡ｰ菫よ焚縺ｮ蜥�
    ! Damping coefficients by horizonatl diffusion and a sponge layer
    !
    allocate( wz_DisCoefM(lmax, 1:kmax) )
    allocate( wz_DisCoefH(lmax, 1:kmax) )
    allocate( w_DisCoefQ (lmax) )

    do k = 1, kmax
      wz_DisCoefM(:,k) = w_HDifCoefM - wz_SpongeLayerCoefM(:,k)
      wz_DisCoefH(:,k) = w_HDifCoefH - wz_SpongeLayerCoefH(:,k)
    end do
    w_DisCoefQ = w_HDifCoefH

    ! Calculation of divergence damping period
    !
    DivDampPeriod = DCCalConvertByUnit( DivDampPeriodValue, DivDampPeriodUnit, 'sec' )
    if ( DivDampTime0 > TimeN ) then
      call MessageNotify( 'E', module_name, &
        & 'DivDampTime0, %f, must be less than or equal to TimeN, %f.', &
        & d = (/ DivDampTime0, TimeN /) )
    end if

    ! 繝偵せ繝医Μ繝��繧ｿ蜃ｺ蜉帙�縺溘ａ縺ｮ縺ｸ縺ｮ螟画焚逋ｻ骭ｲ
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'DUDtDyn', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'dynamical tendency of zonal wind', 'm s-2' )
    call HistoryAutoAddVariable( 'DVDtDyn', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'dynamical tendency of meridional wind', 'm s-2' )
    call HistoryAutoAddVariable( 'DTempDtDyn', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'dynamical tendency of temperature', 'K s-1' )
    do n = 1, ncmax
      call HistoryAutoAddVariable( 'D'//trim(a_QMixName(n))//'DtDyn', &
        & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
        & 'dynamical tendency of water vapor', 's-1' )
    end do
    call HistoryAutoAddVariable( 'DPsDtDyn', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'dynamical tendency of surface pressure', 'Ps s-1' )
    call HistoryAutoAddVariable( 'OMG', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'vertical velocity in pressure coordinate (omega, DP/Dt)', 'Pa s-1' )

    call HistoryAutoAddVariable( 'SigDot', &
      & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
      & 'sigma-vertical velocity', '1 s-1' )
    call HistoryAutoAddVariable( 'DPiDt', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'Pi (log Ps) tendency)', 'Pa s-1' )

    call HistoryAutoAddVariable( 'Vor', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'vorticity', 's-1' )
    call HistoryAutoAddVariable( 'Div', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'divergence', 's-1' )

    call HistoryAutoAddVariable( 'Mass', &
      & (/ AxNameX, AxNameY, AxNameT /), &
      & 'mass', 'kg' )
    call HistoryAutoAddVariable( 'KinEngy', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'kinetic energy', 'J m-2' )
    call HistoryAutoAddVariable( 'IntEngy', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'internal energy', 'J m-2' )
    call HistoryAutoAddVariable( 'PotEngy', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'potential energy', 'J m-2' )
    call HistoryAutoAddVariable( 'LatEngy', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'latent energy', 'J m-2' )
    call HistoryAutoAddVariable( 'TotEngy', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'total energy', 'J m-2' )
    call HistoryAutoAddVariable( 'Enstro', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'enstrophy', 'kg m-2 s-2' )

    call HistoryAutoAddVariable( 'AdvTestStreamFunc', &
      & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
      & 'mass stream function', 'kg s-1' )

    ! Initialization of modules used in this module
    !

    ! 雉ｪ驥上�陬懈ｭ｣
    ! Mass fixer
    !
    call MassFixerInit

    if ( FlagAdvTest ) then
      !
      ! Utility module for advection test
      !
      call AdvTestInit
    end if

    if ( FlagSLTT ) then

!!$    select case ( IDTTMethod )
!!$    case ( IDTTMethodSp )
!!$    case ( IDTTMethodSL )
      ! 繧ｻ繝溘Λ繧ｰ繝ｩ繝ｳ繧ｸ繝･豕輔↓繧医ｋ迚ｩ雉ｪ遘ｻ豬�
      ! Semi-Lagrangian method for tracer transport
      !
      call SLTTInit
!!$    case ( IDTTMethodFV )
!!$      ! 迚ｩ雉ｪ遘ｻ豬� (譛蛾剞菴鍋ｩ肴ｳ�)
!!$      ! Tracer Transport (finite volume method)
!!$      !
!!$      call TTFVInit
!!$    end select

    end if

    ! 蜊ｰ蟄� ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  TimeIntegScheme             = %c', c1 = trim( TimeIntegScheme ) )
    call MessageNotify( 'M', module_name, '  HDEFoldTime                 = %f [%c]', &
      & d = (/ HDEFoldTimeValue /), c1 = trim(HDEFoldTimeUnit) )
    call MessageNotify( 'M', module_name, '  HDOrder                     = %d', i = (/ HDOrder /) )
    call MessageNotify( 'M', module_name, '  VisCoef                     = %f', d = (/ VisCoef /) )
    call MessageNotify( 'M', module_name, '  FlagSpongeLayer             = %b', l = (/ FlagSpongeLayer /) )
    call MessageNotify( 'M', module_name, '  FlagSpongeLayerforZonalMean = %b', l = (/ FlagSpongeLayerforZonalMean /) )
    call MessageNotify( 'M', module_name, '  FlagSpongeLayerforHeat      = %b', l = (/ FlagSpongeLayerforHeat /) )
    call MessageNotify( 'M', module_name, '  SLEFoldTime                 = %f [%c]', &
      & d = (/ SLEFoldTimeValue /), c1 = trim(SLEFoldTimeUnit) )
    call MessageNotify( 'M', module_name, '  SLOrder                     = %d', i = (/ SLOrder /) )
    call MessageNotify( 'M', module_name, '  SLNumLayer                  = %d', i = (/ SLNumLayer /) )
    call MessageNotify( 'M', module_name, '  RefTemp                     = %f', d = (/ RefTemp /) )
    call MessageNotify( 'M', module_name, '  FlagDivDamp                 = %b', l = (/ FlagDivDamp /) )
    call MessageNotify( 'M', module_name, '  DivDampPeriod               = %f', d = (/ DivDampPeriod /) )
    call MessageNotify( 'M', module_name, '  DivDampTime0                = %f', d = (/ DivDampTime0 /) )
    call MessageNotify( 'M', module_name, '  FlagTotMassFix              = %b', l = (/ FlagTotMassFix /) )
    call MessageNotify( 'M', module_name, '  FlagMassHorDifCor           = %b', l = (/ FlagMassHorDifCor /) )
    call MessageNotify( 'M', module_name, '  FlagSLTT                    = %b', l = (/ FlagSLTT /) )
!!$    call MessageNotify( 'M', module_name, '  TracerTransportMethod       = %c', c1 = trim( TracerTransportMethod ) )
    call MessageNotify( 'M', module_name, '  FlagAdvTest                 = %b', l = (/ FlagAdvTest /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    dynamics_hspl_vas83_inited = .true.

  end subroutine DynamicsHSplVAS83Init

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

  subroutine DynamicsHSplVAS83Finalize
    !
    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蜀�Κ縺ｮ螟画焚縺ｮ蜑ｲ繧贋ｻ倥￠隗｣髯､繧定｡後＞縺ｾ縺�.
    !
    ! Deallocate variables in this module.
    !

    ! 螳｣險譁� ; Declaration statements
    !
    implicit none

    ! 螳溯｡梧枚 ; Executable statement
    !

    if ( .not. dynamics_hspl_vas83_inited ) return

    ! 繝�ヵ繧ｩ繝ｫ繝亥､縺ｸ謌ｻ縺�
    ! Return to default values
    !
    DelTimeSave = - 1.0_DP

    ! 蜑ｲ繧贋ｻ倥￠隗｣髯､
    ! Deallocation
    !
    if ( allocated( xy_SinLat ) )    deallocate( xy_SinLat )
    if ( allocated( xy_CosLat ) )    deallocate( xy_CosLat )

    if ( allocated( xy_Cori ) )      deallocate( xy_Cori )

    if ( allocated( z_HydroAlpha ) ) deallocate( z_HydroAlpha )
    if ( allocated( z_HydroBeta  ) ) deallocate( z_HydroBeta  )

    if ( allocated( z_TInpCoefA ) )  deallocate( z_TInpCoefA )
    if ( allocated( z_TInpCoefB ) )  deallocate( z_TInpCoefB )
    if ( allocated( z_TInpCoefK ) )  deallocate( z_TInpCoefK )

    if ( allocated( z_RefTemp ) )    deallocate( z_RefTemp )
    if ( allocated( r_RefTemp ) )    deallocate( r_RefTemp )

    if ( allocated( w_LaplaEigVal ) ) deallocate( w_LaplaEigVal )

    if ( allocated( wz_DisCoefM   ) ) deallocate( wz_DisCoefM     )
    if ( allocated( wz_DisCoefH   ) ) deallocate( wz_DisCoefH     )
    if ( allocated( w_DisCoefQ    ) ) deallocate( w_DisCoefQ      )

    if ( allocated( z_siMtxC      ) )  deallocate( z_siMtxC       )
    if ( allocated( z_siMtxG      ) )  deallocate( z_siMtxG       )
    if ( allocated( zz_siMtxH     ) )  deallocate( zz_siMtxH      )
    if ( allocated( zz_siMtxDiH   ) )  deallocate( zz_siMtxDiH    )
    if ( allocated( wzz_siMtxWDiH ) )  deallocate( wzz_siMtxWDiH  )
    if ( allocated( zz_siMtxGCt   ) )  deallocate( zz_siMtxGCt    )

    if ( allocated( zz_siMtxW  ) )  deallocate( zz_siMtxW  )
    if ( allocated( zz_siMtxQ  ) )  deallocate( zz_siMtxQ  )
    if ( allocated( zz_siMtxS  ) )  deallocate( zz_siMtxS  )
    if ( allocated( zz_siMtxR  ) )  deallocate( zz_siMtxR  )

!!$    if ( allocated(nmo           ) )  deallocate( nmo            )
!!$    if ( allocated(wzz_siMtxM    ) )  deallocate( wzz_siMtxM     )
!!$    if ( allocated(z_siMtxPivWork) )  deallocate( z_siMtxPivWork )
    if ( allocated(wzz_siMtxLU   ) )  deallocate( wzz_siMtxLU    )
    if ( allocated(wz_siMtxPiv   ) )  deallocate( wz_siMtxPiv    )

    dynamics_hspl_vas83_inited = .false.

  end subroutine DynamicsHSplVAS83Finalize

  !-------------------------------------------------------------------
  ! Memo: YOT (2009/10/04)
  ! The routine, SemiImplMatrix_old, is old version of SemiImplMatrix.
  ! This is not deleted now, since this may be valuable as a reference in
  ! optimizing the code.
  !
!!$
!!$  subroutine SemiImplMatrix_old
!!$    !
!!$    ! TimeIntegration 縺ｧ菴ｿ逕ｨ縺吶ｋ菫よ焚縺ｮ險ｭ螳壹ｒ陦後＞縺ｾ縺�.
!!$    ! 蛻晏屓縺翫ｈ縺ｳ $ \Delta t $ 縺悟､画峩縺輔ｌ縺溷�ｴ蜷井ｻ･螟悶�,
!!$    ! 蜑榊屓縺ｫ險ｭ螳壹＠縺溷､繧偵◎縺ｮ縺ｾ縺ｾ逕ｨ縺�∪縺�.
!!$    !
!!$    ! Configure coefficients for "TimeIntegration".
!!$    ! Setting values that are set last time are used,
!!$    ! except when first time and ﾎ杯 are changed.
!!$    !
!!$
!!$    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; 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 谺｡譁ｹ遞句ｼ上ｒ隗｣縺上◆繧√�髢｢謨ｰ (spml 蜷梧｢ｱ繝｢繧ｸ繝･繝ｼ繝ｫ)
!!$    ! Functions to solve linear simultaneous equation by LU decomposition
!!$    ! (a module included in spml)
!!$    !
!!$    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.
!!$    !
!!$    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 ( DelTimeSave == DelTime ) return
!!$    DelTimeSave = DelTime
!!$
!!$    call DbgMessage( '%c: %c: (DelTime=%f [sec]) coefficients for "TimeIntegration" is generated. ', &
!!$      & c1 = module_name, c2 = 'SemiImplMatrix', d = (/ DelTime /) )
!!$
!!$    ! TimeIntegration 縺ｧ菴ｿ逕ｨ縺吶ｋ菫よ焚縺ｮ險ｭ螳�
!!$    ! Configure coefficients for "TimeIntegration"
!!$    !
!!$    if ( .not. allocated( z_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_TInpCoefK * 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.
!!$    do k = 1, kmax
!!$      do l = 1, k
!!$        zz_siMtxW(k,l) = CpDry * z_HydroAlpha(l)
!!$      enddo
!!$
!!$      do l = 1, k-1
!!$        zz_siMtxW(k,l) = zz_siMtxW(k,l) + CpDry * z_HydroBeta(l)
!!$      enddo
!!$    enddo
!!$
!!$    zz_siMtxS = 0.
!!$    do k = 1, kmax
!!$      do l = 1, kmax
!!$        zz_siMtxS(k,l) = r_Sigma(k-1) * z_DelSigma(l)
!!$      enddo
!!$      do l = k, kmax
!!$        zz_siMtxS(k,l) = zz_siMtxS(k,l) - z_DelSigma(l)
!!$      enddo
!!$    enddo
!!$
!!$    zz_siMtxQ = 0.
!!$    do k = 1, kmax
!!$      zz_siMtxQ(k,k) = ( r_RefTemp(k-1) - z_RefTemp(k) ) / z_DelSigma(k)
!!$    enddo
!!$    do k = 1, kmax-1
!!$      zz_siMtxQ(k,k+1) = ( z_RefTemp(k) - r_RefTemp(k) ) / z_DelSigma(k)
!!$    enddo
!!$
!!$    zz_siMtxQS = 0.
!!$    zz_siMtxQS = matmul(zz_siMtxQ, zz_siMtxS)
!!$
!!$    zz_siMtxR = 0.
!!$    do k = 1, kmax
!!$      do l = k, kmax
!!$        zz_siMtxR(k,l) = &
!!$          & - z_HydroAlpha(k) / z_DelSigma(k) * z_DelSigma(l) * z_RefTemp(k)
!!$      enddo
!!$      do l = k + 1, kmax
!!$        zz_siMtxR(k,l) = zz_siMtxR(k,l)  &
!!$          & - z_HydroBeta(k) / z_DelSigma(k) * z_DelSigma(l) * z_RefTemp(k)
!!$      enddo
!!$    enddo
!!$
!!$    zz_siMtxH = 0.
!!$    zz_siMtxH = zz_siMtxQS - zz_siMtxR
!!$
!!$    zz_siMtxWH = 0.
!!$    zz_siMtxWH = matmul(zz_siMtxW, zz_siMtxH)
!!$
!!$    if ( .not. allocated(nmo           ) )  allocate( nmo           (1:2, 0:nmax, 0:nmax) )
!!$    if ( .not. allocated(wzz_siMtxM    ) )  allocate( wzz_siMtxM    ((nmax+1)**2, 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
!!$
!!$      ! 繧ｹ繝壹け繝医Ν豺ｻ蟄鈴��ｺ上�蜿悶ｊ蜃ｺ縺�
!!$      ! 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
!!$          wzz_siMtxM ( nm,k,kk ) = &
!!$            & - DelTime**2 * wz_LaplaEigVal(nm,1) &
!!$            &   * (   zz_siMtxWH( k,kk ) &
!!$            &       + zz_siMtxGCt( k,kk ) &
!!$            &         * ( 1. - 2. * DelTime * wz_HDifCoefH(nm,1) )  )
!!$          if ( k == kk ) then
!!$            wzz_siMtxM ( nm,k,kk ) = &
!!$              & wzz_siMtxM ( nm,k,kk ) &
!!$              & +   ( 1. - 2. * DelTime * wz_HDifCoefM(nm,1) ) &
!!$              &   * ( 1. - 2. * DelTime * wz_HDifCoefH(nm,1) )
!!$          endif
!!$        end do
!!$      end do
!!$
!!$      ! LU 陦悟�險育ｮ�
!!$      ! LU matrix calculation
!!$      !
!!$      call LUDecomp( &
!!$        & wzz_siMtxM(nm,:,:), & ! (inout)
!!$        & z_siMtxPivWork )      ! (out)
!!$
!!$      ! 繝繝溘�蛟､縺ｮ莉｣蜈･. (菴咲ｽｮ kmax 縺ｯ縺ｾ縺�譛ｪ螳夂ｾｩ縺ｪ縺溘ａ).
!!$      ! Dummy value is subtituted. (Because position kmax is undefined yet).
!!$      !
!!$      z_siMtxPivWork(kmax) = 0
!!$
!!$      ! 驟榊�縺ｮ隧ｰ繧∵崛縺�
!!$      ! 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 ) = wzz_siMtxM ( nm,k,kk )
!!$              wzz_siMtxLU ( nmo(2,m,l),k,kk ) = wzz_siMtxM ( nm,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 ) = wzz_siMtxM ( nm,k,kk )
!!$        end do
!!$        wz_siMtxPiv ( nm,k ) = z_siMtxPivWork ( k )
!!$      end do
!!$    end do
!!$
!!$  end subroutine SemiImplMatrix_old
!!$

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

  subroutine SemiImplMatrix
    !
    ! TimeIntegration 縺ｧ菴ｿ逕ｨ縺吶ｋ菫よ焚縺ｮ險ｭ螳壹ｒ陦後＞縺ｾ縺�.
    ! 蛻晏屓縺翫ｈ縺ｳ $ \Delta t $ 縺悟､画峩縺輔ｌ縺溷�ｴ蜷井ｻ･螟悶�,
    ! 蜑榊屓縺ｫ險ｭ螳壹＠縺溷､繧偵◎縺ｮ縺ｾ縺ｾ逕ｨ縺�∪縺�.
    !
    ! Configure coefficients for "TimeIntegration".
    ! Setting values that are set last time are used,
    ! except when first time and $ \Delta t $ are changed.
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; USE statements
    !

    ! 迚ｩ逅�ｮ壽焚險ｭ螳�
    ! Physical constants settings
    !
    use constants, only: &
      & RPlanet, &
                              ! $ a $ [m].
                              ! 諠第弌蜊雁ｾ�.
                              ! Radius of planet
      & CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 荵ｾ辯･螟ｧ豌励�螳壼悸豈皮�.
                              ! Specific heat of air at constant pressure

    ! 譬ｼ蟄千せ險ｭ螳�
    ! Grid points settings
    !
    use gridset, only: lmax, & ! 繧ｹ繝壹け繝医Ν繝��繧ｿ縺ｮ驟榊�繧ｵ繧､繧ｺ
                               ! Size of array for spectral data
      &                imax, & ! 邨悟ｺｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in longitude
      &                jmax, & ! 邱ｯ蠎ｦ譬ｼ蟄千せ謨ｰ.
                               ! Number of grid points in latitude
      &                kmax    ! 驩帷峩螻､謨ｰ.
                               ! Number of vertical level

    ! 蠎ｧ讓吶ョ繝ｼ繧ｿ險ｭ螳�
    ! 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)
    !
#ifdef LIB_MPI
#ifdef SJPACK
    use wa_mpi_module_sjpack, only: l_nm
#else
    use wa_mpi_module, only: l_nm
#endif
#elif AXISYMMETRY
    use wa_zonal_module, only: l_nm
#elif SJPACK
    use wa_module_sjpack, only: l_nm
#elif AXISYMMETRY_SJPACK
    use wa_zonal_module_sjpack, only: l_nm
#else
    use wa_module, only: l_nm
#endif

    ! LU 蛻�ｧ｣豕輔↓繧医ｊ騾｣遶� 1 谺｡譁ｹ遞句ｼ上ｒ隗｣縺上◆繧√�髢｢謨ｰ (spml 蜷梧｢ｱ繝｢繧ｸ繝･繝ｼ繝ｫ)
    ! Functions to solve linear simultaneous equation by LU decomposition
    ! (a module included in spml)
    !
    use lumatrix, only: LUDecomp

    ! 繝�ヰ繝�げ逕ｨ繝ｦ繝ｼ繝�ぅ繝ｪ繝�ぅ
    ! Utilities for debug
    !
    use dc_trace, only: DbgMessage, BeginSub, EndSub, Debug

    ! 螳｣險譁� ; Declaration statements
    !
    implicit none

    ! TimeIntegration 遲峨〒菴ｿ逕ｨ縺吶ｋ菫よ焚縺ｮ險ｭ螳壹�縺溘ａ縺ｮ菴懈･ｭ螟画焚
    ! Work variable for coefficients for "TimeIntegration", etc.
    !
    integer:: k, l, kk        ! 驩帷峩譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in vertical direction
    integer:: n

    ! 螳溯｡梧枚 ; Executable statement
    !

    ! $ \Delta t $ [s] 縺ｮ繝√ぉ繝�け繝ｻ菫晏ｭ�
    ! Check and save $ \Delta t $ [s]
    !
    if ( DelTimeSave == DelTime ) return
    DelTimeSave = DelTime

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

    ! TimeIntegration 縺ｧ菴ｿ逕ｨ縺吶ｋ菫よ焚縺ｮ險ｭ螳�
    ! Configure coefficients for "TimeIntegration"
    !
    if ( .not. allocated( z_siMtxC      ) )  allocate( z_siMtxC     (1:kmax) )
    if ( .not. allocated( z_siMtxG      ) )  allocate( z_siMtxG     (1:kmax) )
    if ( .not. allocated( zz_siMtxH     ) )  allocate( zz_siMtxH    (1:kmax, 1:kmax) )
    if ( .not. allocated( wz_siMtxDi    ) )  allocate( wz_siMtxDi   (lmax,1:kmax) )
    if ( .not. allocated( zz_siMtxDiH   ) )  allocate( zz_siMtxDiH  (1:kmax, 1:kmax) )
    if ( .not. allocated( wzz_siMtxWDiH ) )  allocate( wzz_siMtxWDiH(lmax,1:kmax, 1:kmax) )
    if ( .not. allocated( zz_siMtxGCt   ) )  allocate( zz_siMtxGCt  (1:kmax, 1:kmax) )

    if ( .not. allocated( zz_siMtxW  ) )  allocate( zz_siMtxW  (1:kmax, 1:kmax) )
    if ( .not. allocated( zz_siMtxQ  ) )  allocate( zz_siMtxQ  (1:kmax, 1:kmax) )
    if ( .not. allocated( zz_siMtxS  ) )  allocate( zz_siMtxS  (1:kmax, 1:kmax) )
    if ( .not. allocated( zz_siMtxR  ) )  allocate( zz_siMtxR  (1:kmax, 1:kmax) )

    z_siMtxC = z_DelSigma
    z_siMtxG = CpDry * z_TInpCoefK * z_RefTemp

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

    zz_siMtxW = 0.
    do k = 1, kmax
      do l = 1, k
        zz_siMtxW(k,l) = CpDry * z_HydroAlpha(l)
      enddo

      do l = 1, k-1
        zz_siMtxW(k,l) = zz_siMtxW(k,l) + CpDry * z_HydroBeta(l)
      enddo
    enddo

    zz_siMtxS = 0.
    do k = 1, kmax
      do l = 1, kmax
        zz_siMtxS(k,l) = r_Sigma(k-1) * z_DelSigma(l)
      enddo
      do l = k, kmax
        zz_siMtxS(k,l) = zz_siMtxS(k,l) - z_DelSigma(l)
      enddo
    enddo

    zz_siMtxQ = 0.
    do k = 1, kmax
      zz_siMtxQ(k,k) = ( r_RefTemp(k-1) - z_RefTemp(k) ) / z_DelSigma(k)
    enddo
    do k = 1, kmax-1
      zz_siMtxQ(k,k+1) = ( z_RefTemp(k) - r_RefTemp(k) ) / z_DelSigma(k)
    enddo

    zz_siMtxR = 0.
    do k = 1, kmax
      do l = k, kmax
        zz_siMtxR(k,l) = &
          & - z_HydroAlpha(k) / z_DelSigma(k) * z_DelSigma(l) * z_RefTemp(k)
      enddo
      do l = k + 1, kmax
        zz_siMtxR(k,l) = zz_siMtxR(k,l)  &
          & - z_HydroBeta(k) / z_DelSigma(k) * z_DelSigma(l) * z_RefTemp(k)
      enddo
    enddo

    zz_siMtxH = matmul(zz_siMtxQ, zz_siMtxS) - zz_siMtxR

    ! Check of coefficients for horizontal diffusion and sponge layer
    ! A threshold value used below is arbitrary.
    !
    do n = 1, lmax
      do k = 1, kmax
        if ( abs( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM(n,k) ) < 1.0e-10_DP ) then
          call MessageNotify( 'E', module_name, &
            & 'Dissipation coefficient for momentum is inappropriate.' )
        end if
        if ( abs( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefH(n,k) ) < 1.0e-10_DP ) then
          call MessageNotify( 'E', module_name, &
            & 'Dissipation coefficient for heat is inappropriate.' )
        end if
      end do
      if ( abs( 1.0_DP - 2.0_DP * DelTime * w_DisCoefQ(n) ) < 1.0e-10_DP ) then
        call MessageNotify( 'E', module_name, &
          & 'Dissipation coefficient for composition is inappropriate.' )
      end if
    end do

    wz_siMtxDi = 1.0_DP / ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefH )

    do n = 1, lmax
      ! NOTE: wz_siMtiDi is a diagonal matrix.
      !
      do k = 1, kmax
        do kk = 1, kmax
          zz_siMtxDiH(k,kk) = wz_siMtxDi(n,k) * zz_siMtxH(k,kk)
        end do
      end do
      zz_siMtxDiH = matmul( zz_siMtxW, zz_siMtxDiH )
      do k = 1, kmax
        do kk = 1, kmax
          wzz_siMtxWDiH(n,k,kk) = zz_siMtxDiH(k,kk)
        end do
      end do
    end do

    if ( .not. allocated(wzz_siMtxLU) )  allocate( wzz_siMtxLU(lmax, 1:kmax, 1:kmax) )
    if ( .not. allocated(wz_siMtxPiv) )  allocate( wz_siMtxPiv(lmax, 1:kmax) )

    ! 陦悟� $ \underline{M} $ 縺ｮ險育ｮ�
    ! Calculate matrix $ \underline{M} $
    !
    do k = 1, kmax
      do kk = 1, kmax
        wzz_siMtxLU ( :,k,kk ) = &
          & - DelTime**2 * w_LaplaEigVal(:) &
          &   * ( wzz_siMtxWDiH(:,k,kk) + zz_siMtxGCt(k,kk) )
        if ( k == kk ) then
          wzz_siMtxLU ( :,k,kk ) = &
            &   wzz_siMtxLU ( :,k,kk ) &
            & + ( 1.0_DP - 2.0_DP * DelTime * wz_DisCoefM (:,k) )
        endif
      end do
    end do

    ! LU 陦悟�險育ｮ�
    ! LU matrix calculation
    !
    call LUDecomp(    &
      & wzz_siMtxLU,  & ! (inout)
      & wz_siMtxPiv )   ! (out)

  end subroutine SemiImplMatrix

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

end module dynamics_hspl_vas83