!= åŠ›å¦éŽç¨‹ (ã‚¹ãƒšã‚¯ãƒˆãƒ«æ³•, 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 Î”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
!!$
!!$    ! åº§æ¨™ãƒ‡ãƒ¼ã‚¿è¨å®š
!!$    ! 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