!= Module DynamicsHEVI
!
! Authors::   譚牙ｱｱ閠穂ｸ譛�(SUGIYAMA Ko-ichiro)
! Version::   $Id: dynamics_hevi_v3.f90,v 1.4 2014/07/08 00:55:26 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2014. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]

module DynamicsHEVI_v3
  !
  ! 蜉帛ｭｦ繧ｳ繧｢.
  !   繧ｿ繧､繝�繧ｹ繝励Μ繝�ヨ豕輔ｒ蛻ｩ逕ｨ. 髻ｳ豕｢繝｢繝ｼ繝峨→縺昴ｌ莉･螟悶ｒ蛻･縲��譎る俣蛻ｻ縺ｿ縺ｧ隗｣縺�.
  !   遏ｭ縺�凾髢薙せ繝�ャ繝励�險育ｮ励↓縺ｯ HE-VI 豕輔ｒ蛻ｩ逕ｨ.
  !
  ! v3: v1 縺ｨ v2 繧呈紛逅�☆繧九％縺ｨ縺ｧ, 險育ｮ玲婿豕輔ｒ蛻�ｊ譖ｿ縺医ｉ繧後ｋ繧医≧縺ｫ隱ｿ謨ｴ縺励◆迚�
  !
  ! Note:
  !  * 繧ｨ繧ｯ繧ｹ繝翫�髢｢謨ｰ縺ｮ遨ｺ髢捺婿蜷代�髮｢謨｣蛹悶� 2 谺｡邊ｾ蠎ｦ縺ｧ縺ゅｋ縺溘ａ, 豌怜悸蛯ｾ蠎ｦ
  !    蜉幃���險育ｮ励�繝ｭ繧ｰ繝ｩ繝�縺ｫ縺翫＞縺ｦ differentiate_center4 繝｢繧ｸ繝･繝ｼ繝ｫ繧�
  !    謖�ｮ壹☆繧九％縺ｨ縺ｯ縺ｧ縺阪↑縺��縺ｧ豕ｨ諢�.

  !繝｢繧ｸ繝･繝ｼ繝ｫ隱ｭ縺ｿ霎ｼ縺ｿ
  use dc_types,   only : DP

  !證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
  implicit none

  !螻樊ｧ縺ｮ謖�ｮ�
  private

  !螟画焚縺ｮ螳夂ｾｩ
  character(*), parameter:: module_name = 'dynamics_hevi'
                                           ! 繝｢繧ｸ繝･繝ｼ繝ｫ縺ｮ蜷咲ｧｰ.
                                           ! Module name
  real(DP), save, private :: FactorBuoyTemp    = 1.0d0
                                           !豬ｮ蜉� (貂ｩ蠎ｦ縺ｮ蟇�ｸ�) 縺ｮ譛臥┌
                                           !閠��縺励↑縺��ｴ蜷医�蛟､繧偵ぞ繝ｭ縺ｫ縺吶ｋ.
  real(DP), save, private :: FactorBuoyMolWt   = 1.0d0
                                           !豬ｮ蜉� (蛻�ｭ宣㍼蜉ｹ譫�) 縺ｮ譛臥┌
                                           !閠��縺励↑縺��ｴ蜷医�蛟､繧偵ぞ繝ｭ縺ｫ縺吶ｋ.
  real(DP), save, private :: FactorBuoyLoading = 1.0d0
                                           !豬ｮ蜉� (闕ｷ驥榊柑譫�) 縺ｮ譛臥┌
                                           !閠��縺励↑縺��ｴ蜷医�蛟､繧偵ぞ繝ｭ縺ｫ縺吶ｋ.
  real(DP), save, private :: FactorDExnerDtAdv    = 0.0d0
                                           !繧ｨ繧ｯ繧ｹ繝翫�髢｢謨ｰ縺ｮ遘ｻ豬√�譛臥┌
                                           !閠��縺励↑縺��ｴ蜷医�蛟､繧偵ぞ繝ｭ縺ｫ縺吶ｋ.
  real(DP), save, private :: FactorDExnerDtExpnd  = 0.0d0
                                           !繧ｨ繧ｯ繧ｹ繝翫�髢｢謨ｰ縺ｮ閹ｨ蠑ｵ鬆��譛臥┌
                                           !閠��縺励↑縺��ｴ蜷医�蛟､繧偵ぞ繝ｭ縺ｫ縺吶ｋ.
  real(DP), save, private :: FactorExnerFall  = 0.0d0
                                           !關ｽ荳九↓莨ｴ縺�ｻｮ貂ｩ菴榊､牙喧縺ｫ莨ｴ縺�お繧ｯ繧ｹ繝翫�髢｢謨ｰ縺ｮtendency
                                           !閠��縺励↑縺��ｴ蜷医�蛟､繧偵ぞ繝ｭ縺ｫ縺吶ｋ.
  integer, save, private  :: IDAdvection               = 0
                                           ! 遘ｻ豬∬ｨ育ｮ玲婿豕�
  integer, parameter      :: IDAdvectionCenter4_std    = 1
                                           ! 1: 4 谺｡荳ｭ螟ｮ蟾ｮ蛻� (蠕ｮ蛻�ｹｳ蝮�Δ繧ｸ繝･繝ｼ繝ｫ蛻ｩ逕ｨ)
  integer, parameter      :: IDAdvectionCenter4_2D_dry = 2
                                           ! 2: 4 谺｡荳ｭ螟ｮ蟾ｮ蛻� (2D 迚郁ｨ育ｮ励Ν繝ｼ繝√Φ)
  integer, parameter      :: IDAdvectionCenter4_2D     = 3
                                           ! 2: 4 谺｡荳ｭ螟ｮ蟾ｮ蛻� (2D 迚郁ｨ育ｮ励Ν繝ｼ繝√Φ)
  integer, parameter      :: IDAdvectionCenter4_3D_dry = 4
                                           ! 3: 4 谺｡荳ｭ螟ｮ蟾ｮ蛻� (3D 迚郁ｨ育ｮ励Ν繝ｼ繝√Φ)
  integer, parameter      :: IDAdvectionCenter4_3D      = 5
                                           ! 3: 4 谺｡荳ｭ螟ｮ蟾ｮ蛻� (3D 迚郁ｨ育ｮ励Ν繝ｼ繝√Φ)
  integer, save, private  :: IDAcousticmode            = 0
                                           ! 遏ｭ縺�凾髢薙せ繝�ャ繝励�險育ｮ玲婿豕�
  integer, parameter      :: IDAcousticmode_std        = 1
                                           ! 1: 蠕ｮ蛻�ｹｳ蝮�Δ繧ｸ繝･繝ｼ繝ｫ縺ｮ蛻ｩ逕ｨ.
  integer, parameter      :: IDAcousticmode_2D         = 2
                                           ! 2: 2D 迚郁ｨ育ｮ励Ν繝ｼ繝√Φ
  integer, parameter      :: IDAcousticmode_3D         = 3
                                           ! 3: 3D 迚郁ｨ育ｮ励Ν繝ｼ繝√Φ

  ! public
  !
  public Dynamics_Init
  public Dynamics_Long_forcing
  public Dynamics_Short_forcing

contains

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  subroutine Dynamics_Init
    !
    ! 蜉帛ｭｦ繧ｳ繧｢ 3D 迚医�蛻晄悄蛹悶Ν繝ｼ繝√Φ
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ隱ｭ縺ｿ霎ｼ縺ｿ
    !
    use dc_types,   only : DP, STRING
    use dc_iounit,  only : FileOpen
    use dc_message, only : MessageNotify
    use gridset,    only : FlagCalc3D, FlagCalcMoist
    use namelist_util, &
      &             only : namelist_filename
    use advection_center4_std, &
      &             only : advection_center4_std_init
    use advection_center4_2D, &
      &             only : advection_center4_2D_init
    use advection_center4_3D, &
      &             only : advection_center4_3D_init
    use acousticmode_std, &
      &             only : acousticmode_std_init
    use acousticmode_2D, &
      &             only : acousticmode_2D_init
    use acousticmode_3D, &
      &             only : acousticmode_3D_init

    !證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
    !
    implicit none

    ! 螟画焚縺ｮ螳夂ｾｩ
    !
    real(DP)          :: AlphaSound = 5.0d-2    !髻ｳ豕｢貂幄｡ｰ鬆��菫よ焚 (豌苓ｱ｡蠎∵焚蛟､莠亥�ｱ隱ｲ蝣ｱ蜻翫�蛻･蜀�49 繧医ｊ)
    real(DP)          :: AlphaNDiff = 1.0d-3    !4谺｡縺ｮ謨ｰ蛟､諡｡謨｣縺ｮ菫よ焚. CReSS 繝槭ル繝･繧｢繝ｫ繧医ｊ
    real(DP)          :: NDiffRatio = 1.0d0     !騾溷ｺｦ縺ｫ蟇ｾ縺吶ｋ邊俶ｧ繧剃ｸ翫￡繧句�ｴ蜷医�謨ｰ蟄励ｒ 1 莉･荳翫↓縺吶ｋ.
    integer           :: unit                   !陬�ｽｮ逡ｪ蜿ｷ
    character(STRING) :: FlagAdvection = ""     !遘ｻ豬∬ｨ育ｮ玲婿豕輔�驕ｸ謚�
    character(STRING) :: FlagAcousticmode = ""  !髻ｳ豕｢繝｢繝ｼ繝峨�險育ｮ玲婿豕輔�驕ｸ謚�

    !-------------------------------------------------------------------
    ! Namelist 縺九ｉ諠��ｱ繧貞叙蠕励☆繧�
    !
    NAMELIST /Dynamics_nml/                                 &
      & AlphaSound, AlphaNDiff, NDiffRatio,                 &
      & FactorBuoyTemp, FactorBuoyMolWt, FactorBuoyLoading, &
      & FactorDExnerDtAdv, FactorDExnerDtExpnd, FactorExnerFall,  &
      & FlagAdvection, FlagAcousticmode

    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=dynamics_nml)
    close(unit)

    !-------------------------------------------------------------------
    ! 繝輔Λ繧ｰ蜃ｦ逅�
    !
    if ( FlagAdvection == "Center4" .OR. FlagAdvection == "" ) then
      !
      ! 4 谺｡邊ｾ蠎ｦ荳ｭ蠢�ｷｮ蛻�. 蠕ｮ蛻�ｹｳ蝮�Δ繧ｸ繝･繝ｼ繝ｫ繧剃ｽｿ繧上↑縺��ｴ蜷�. (繝�ヵ繧ｩ繝ｫ繝�)
      !
      if ( FlagCalc3D ) then
!        if ( FlagCalcTracer .OR. FlagCalcMoist ) then
        if ( FlagCalcMoist ) then
          IDAdvection = IDAdvectionCenter4_3D
        else
          IDAdvection = IDAdvectionCenter4_3D_dry
        end if
      else
!        if ( FlagCalcTracer .OR. FlagCalcMoist ) then
        if ( FlagCalcMoist ) then
          IDAdvection = IDAdvectionCenter4_2D
        else
          IDAdvection = IDAdvectionCenter4_2D_dry
        end if
      end if
    else
      !
      ! 4 谺｡邊ｾ蠎ｦ荳ｭ蠢�ｷｮ蛻� with 蠕ｮ蛻�ｹｳ蝮�Δ繧ｸ繝･繝ｼ繝ｫ
      !
      IDAdvection = IDAdvectionCenter4_std
    end if

    if ( FlagAcousticmode == "Center2" .OR. FlagAcousticmode == "" ) then
      !
      ! 蠕ｮ蛻�ｹｳ蝮�Δ繧ｸ繝･繝ｼ繝ｫ繧剃ｽｿ繧上↑縺��ｴ蜷�. (繝�ヵ繧ｩ繝ｫ繝�)
      !
      if ( FlagCalc3D ) then
        IDAcousticmode = IDAcousticmode_3D
      else
        IDAcousticmode = IDAcousticmode_2D
      end if
    else
      !
      ! with 蠕ｮ蛻�ｹｳ蝮�Δ繧ｸ繝･繝ｼ繝ｫ
      !
      IDAcousticmode = IDAcousticmode_std
    end if

    !-------------------------------------------------------------------
    ! 遘ｻ豬∬ｨ育ｮ礼畑縺ｮ險育ｮ励Δ繧ｸ繝･繝ｼ繝ｫ縺ｮ蛻晄悄蛹�:
    !
    select case ( IDAdvection )

    case ( IDAdvectionCenter4_std )
      call advection_center4_std_init( AlphaNDiff, NDiffRatio )

    case ( IDAdvectionCenter4_3D_dry, IDAdvectionCenter4_3D )
      call advection_center4_3D_init( AlphaNDiff, NDiffRatio )

    case ( IDAdvectionCenter4_2D_dry, IDAdvectionCenter4_2D )
      call advection_center4_2D_init( AlphaNDiff, NDiffRatio )

    end select

    !-------------------------------------------------------------------
    ! 髻ｳ豕｢繝｢繝ｼ繝峨�險育ｮ励Δ繧ｸ繝･繝ｼ繝ｫ縺ｮ蛻晄悄蛹�:
    !
    select case ( IDAcousticmode )

    case ( IDAcousticmode_std )
      call acousticmode_std_init( AlphaSound )

    case ( IDAcousticmode_3D )
      call acousticmode_3D_init( AlphaSound )

    case ( IDAcousticmode_2D )
      call acousticmode_2D_init( AlphaSound )

    end select

    !-------------------------------------------------------------------
    ! tendency 縺ｮ蜃ｺ蜉�
    !
    call Dynamics_Tendency_output

    !-------------------------------------------------------------------
    ! 蜃ｺ蜉�
    !
    call MessageNotify( "M", module_name, "AlphaSound = %f", d=(/AlphaSound/) )
    call MessageNotify( "M", module_name, "AlphaNDiff = %f", d=(/AlphaNDiff/) )

    call MessageNotify( "M", module_name, "FactorBuoyTemp   = %f", d=(/FactorBuoyTemp/) )
    call MessageNotify( "M", module_name, "FactorBuoyMolWt  = %f", d=(/FactorBuoyMolWt/) )
    call MessageNotify( "M", module_name, "FactorBuoyLoading= %f", d=(/FactorBuoyLoading/) )
    call MessageNotify( "M", module_name, "FactorDExnerDtAdv   = %f", d=(/FactorDExnerDtAdv/) )
    call MessageNotify( "M", module_name, "FactorDExnerDtExpnd = %f", d=(/FactorDExnerDtExpnd/) )

    call MessageNotify( "M", module_name, "IDAdvection = %d", i=(/IDAdvection/) )
    call MessageNotify( "M", module_name, "IDAcousticmode = %d", i=(/IDAcousticmode/) )

  end subroutine Dynamics_Init

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  subroutine Dynamics_Long_forcing(   &
    & pyz_VelXBl,  pyz_VelXNl,        & ! (in)
    & xqz_VelYBl,  xqz_VelYNl,        & ! (in)
    & xyr_VelZBl,  xyr_VelZNl,        & ! (in)
    & xyz_PTempBl, xyz_PTempNl,       & ! (in)
    & xyz_ExnerBl, xyz_ExnerNl,       & ! (in)
    & xyzf_QMixBl, xyzf_QMixNl,       & ! (in)
    & xyz_KmBl,    xyz_KmNl,          & ! (in)
    & pyz_DVelXDtNl,                  & ! (inout)
    & xqz_DVelYDtNl,                  & ! (inout)
    & xyr_DVelZDtNl,                  & ! (inout)
    & xyz_DPTempDtNl,                 & ! (inout)
    & xyz_DExnerDtNl,                 & ! (inout)
    & xyzf_DQMixDtNl,                 & ! (inout)
    & xyz_DKmDtNl                     & ! (inout)
    & )
    !
    ! 蜉帛ｭｦ繧ｳ繧｢: 髟ｷ縺�凾髢薙せ繝�ャ繝励〒隧穂ｾ｡縺吶ｋ鬆��險育ｮ�.
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ隱ｭ縺ｿ霎ｼ縺ｿ
    !
    use dc_types,  only : DP
    use gridset,   only : imin, imax, jmin, jmax, kmin, kmax, ncmax, &
                                                    !驟榊�繧ｵ繧､繧ｺ
      &                   nx, ny, nz,              &!迚ｩ逅��伜沺縺ｮ驟榊�繧ｵ繧､繧ｺ
      &                   FlagCalcMoist !繝輔Λ繧ｰ
    use composition, &
      &             only : SpcWetSymbol
    use timeset,    only : TimeN
    use gtool_historyauto, &
      &             only : HistoryAutoPut
    use advection_center4_std, &
      &             only : advection_center4_std_main
    use advection_center4_2D, &
      &             only : advection_center4_2D_dry, advection_center4_2D_tracer
    use advection_center4_3D, &
      &             only : advection_center4_3D_dry, advection_center4_3D_tracer
    use DExnerDt,   only : xyz_DExnerDt_xyzf

    ! 證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
    !
    implicit none

    ! 螟画焚縺ｮ螳夂ｾｩ
    !
    real(DP), intent(in)    :: pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyz_PTempNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyz_ExnerBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyz_ExnerNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(in)    :: xyzf_QMixNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(in)    :: xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyz_KmNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xyz_DPTempDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xyzf_DQMixDtNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(inout) :: xyz_DKmDtNl(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)                :: pyz_DVelXDtAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xqz_DVelYDtAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyr_DVelZDtAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_DKmDtAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_DExnerDtAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_DPTempDtAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyzf_QMixAdv(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP)                :: xyzf_QMixFall(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP)                :: pyz_VelXnDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xqz_VelYnDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyr_VelZnDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_KmNDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_ExnerNDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_ExnerFall(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_PTempNDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyzf_QMixNDiff(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP)                :: xyr_BuoyT(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyr_BuoyM(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyr_BuoyD(imin:imax,jmin:jmax,kmin:kmax)

    integer                 :: f

    !--------------------------------------------------------------------
    ! 遘ｻ豬����險育ｮ�
    !
    select case ( IDAdvection )

    case ( IDAdvectionCenter4_std )

      call advection_center4_std_main(     &
        & pyz_VelXBl,   pyz_VelXNl,        & ! (in)
        & xqz_VelYBl,   xqz_VelYNl,        & ! (in)
        & xyr_VelZBl,   xyr_VelZNl,        & ! (in)
        & xyz_PTempBl,  xyz_PTempNl,       & ! (in)
        & xyz_ExnerBl,  xyz_ExnerNl,       & ! (in)
        & xyzf_QMixBl,  xyzf_QMixNl,       & ! (in)
        & xyz_KmBl,     xyz_KmNl,          & ! (in)
        & pyz_DVelXDtAdv,  pyz_VelXnDiff,     & !(out)
        & xqz_DVelYDtAdv,  xqz_VelYnDiff,     & !(out)
        & xyr_DVelZDtAdv,  xyr_VelZnDiff,     & !(out)
        & xyz_DPTempDtAdv, xyz_PTempNDiff,    & !(out)
        & xyz_DExnerDtAdv, xyz_ExnerNDiff,    & !(out)
        & xyzf_QMixAdv, xyzf_QMixNDiff,    & !(out)
        & xyz_DKmDtAdv,    xyz_KmNDiff        & !(out)
        & )

    case ( IDAdvectionCenter4_3D_dry )

      call advection_center4_3D_dry(       &
        & pyz_VelXBl,   pyz_VelXNl,        & ! (in)
        & xqz_VelYBl,   xqz_VelYNl,        & ! (in)
        & xyr_VelZBl,   xyr_VelZNl,        & ! (in)
        & xyz_PTempBl,  xyz_PTempNl,       & ! (in)
        & xyz_ExnerBl,  xyz_ExnerNl,       & ! (in)
        & xyz_KmBl,     xyz_KmNl,          & ! (in)
        & pyz_DVelXDtAdv,  pyz_VelXnDiff,     & !(out)
        & xqz_DVelYDtAdv,  xqz_VelYnDiff,     & !(out)
        & xyr_DVelZDtAdv,  xyr_VelZnDiff,     & !(out)
        & xyz_DPTempDtAdv, xyz_PTempNDiff,    & !(out)
        & xyz_DExnerDtAdv, xyz_ExnerNDiff,    & !(out)
        & xyz_DKmDtAdv,    xyz_KmNDiff        & !(out)
        & )

    case ( IDAdvectionCenter4_3D )

      call advection_center4_3D_dry(       &
        & pyz_VelXBl,   pyz_VelXNl,        & ! (in)
        & xqz_VelYBl,   xqz_VelYNl,        & ! (in)
        & xyr_VelZBl,   xyr_VelZNl,        & ! (in)
        & xyz_PTempBl,  xyz_PTempNl,       & ! (in)
        & xyz_ExnerBl,  xyz_ExnerNl,       & ! (in)
        & xyz_KmBl,     xyz_KmNl,          & ! (in)
        & pyz_DVelXDtAdv,  pyz_VelXnDiff,     & !(out)
        & xqz_DVelYDtAdv,  xqz_VelYnDiff,     & !(out)
        & xyr_DVelZDtAdv,  xyr_VelZnDiff,     & !(out)
        & xyz_DPTempDtAdv, xyz_PTempNDiff,    & !(out)
        & xyz_DExnerDtAdv, xyz_ExnerNDiff,    & !(out)
        & xyz_DKmDtAdv,    xyz_KmNDiff        & !(out)
        & )

      call advection_center4_3D_tracer(       &
        & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & ! (in)
        & xyzf_QMixBl,  xyzf_QMixNl,          & ! (in)
        & xyzf_QMixAdv, xyzf_QMixNDiff        & !(out)
        & )

    case ( IDAdvectionCenter4_2D_dry )

      ! Y 譁ｹ蜷代� tendency 縺ｯ髮ｶ.
      !
      xqz_DVelYDtAdv   = 0.0d0
      xqz_VelYnDiff = 0.0d0

      call advection_center4_2D_dry(       &
        & pyz_VelXBl,   pyz_VelXNl,        & ! (in)
        & xyr_VelZBl,   xyr_VelZNl,        & ! (in)
        & xyz_PTempBl,  xyz_PTempNl,       & ! (in)
        & xyz_ExnerBl,  xyz_ExnerNl,       & ! (in)
        & xyz_KmBl,     xyz_KmNl,          & ! (in)
        & pyz_DVelXDtAdv,  pyz_VelXnDiff,     & !(out)
        & xyr_DVelZDtAdv,  xyr_VelZnDiff,     & !(out)
        & xyz_DPTempDtAdv, xyz_PTempNDiff,    & !(out)
        & xyz_DExnerDtAdv, xyz_ExnerNDiff,    & !(out)
        & xyz_DKmDtAdv,    xyz_KmNDiff        & !(out)
        & )

    case ( IDAdvectionCenter4_2D )

      ! Y 譁ｹ蜷代� tendency 縺ｯ髮ｶ.
      !
      xqz_DVelYDtAdv   = 0.0d0
      xqz_VelYnDiff = 0.0d0

      call advection_center4_2D_dry(       &
        & pyz_VelXBl,   pyz_VelXNl,        & ! (in)
        & xyr_VelZBl,   xyr_VelZNl,        & ! (in)
        & xyz_PTempBl,  xyz_PTempNl,       & ! (in)
        & xyz_ExnerBl,  xyz_ExnerNl,       & ! (in)
        & xyz_KmBl,     xyz_KmNl,          & ! (in)
        & pyz_DVelXDtAdv,  pyz_VelXnDiff,     & !(out)
        & xyr_DVelZDtAdv,  xyr_VelZnDiff,     & !(out)
        & xyz_DPTempDtAdv, xyz_PTempNDiff,    & !(out)
        & xyz_DExnerDtAdv, xyz_ExnerNDiff,    & !(out)
        & xyz_DKmDtAdv,    xyz_KmNDiff        & !(out)
        & )

      call advection_center4_2D_tracer(    &
        & pyz_VelXNl, xyr_VelZNl,          & ! (in)
        & xyzf_QMixBl,  xyzf_QMixNl,       & ! (in)
        & xyzf_QMixAdv, xyzf_QMixNDiff     & !(out)
        & )

    end select

    !--------------------------------------------------------------------
    ! tendency 縺ｮ譖ｴ譁ｰ
    !

    ! 諡｡謨｣菫よ焚
    !
    xyz_DKmDtNl = xyz_DKmDtNl + ( xyz_KmNDiff + xyz_DKmDtAdv )

    ! output
    !
    call HistoryAutoPut(TimeN, 'DKmDtAdv',  xyz_DKmDtAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'DKmDtDiff', xyz_KmNDiff(1:nx,1:ny,1:nz) )

    ! 貂ｩ菴�
    !
    xyz_DPTempDtNl = xyz_DPTempDtNl + ( xyz_PTempNDiff + xyz_DPTempDtAdv )

    ! output
    !
    call HistoryAutoPut(TimeN, 'DPTempDtAdv',  xyz_DPTempDtAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'DPTempDtDiff', xyz_PTempNDiff(1:nx,1:ny,1:nz) )

    ! 豺ｷ蜷域ｯ�
    !
    call QMixfall  ! 關ｽ荳矩��

    xyzf_DQMixDtNl = xyzf_DQMixDtNl + ( xyzf_QMixNDiff + xyzf_QMixAdv + xyzf_QMixFall )

    ! output
    !
    do f = 1, ncmax
      call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(f))//'DtAdv',   &
        & xyzf_QMixAdv(1:nx,1:ny,1:nz,f))
      call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(f))//'DtDiff', &
        & xyzf_QMixNDiff(1:nx,1:ny,1:nz,f))
    end do

    ! 繧ｨ繧ｯ繧ｹ繝翫�髢｢謨ｰ
    !
    xyz_ExnerFall   = xyz_DExnerDt_xyzf( xyzf_QMixFall ) * FactorExnerFall
    xyz_ExnerNDiff  = xyz_ExnerNDiff   * FactorDExnerDtAdv
    xyz_DExnerDtAdv = xyz_DExnerDtAdv  * FactorDExnerDtAdv
    xyz_DExnerDtNl  = xyz_DExnerDtNl + ( xyz_ExnerFall+ xyz_ExnerNDiff + xyz_DExnerDtAdv )

    ! output
    !
    call HistoryAutoPut(TimeN, 'DExnerDtAdv',  xyz_DExnerDtAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'DExnerDtDiff', xyz_ExnerNDiff(1:nx,1:ny,1:nz) )

    ! VelX
    !
    pyz_DVelXDtNl = pyz_DVelXDtNl + ( pyz_VelXNDiff + pyz_DVelXDtAdv )

    call HistoryAutoPut(TimeN, 'DVelXDtAdv',  pyz_DVelXDtAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'DVelXDtDiff', pyz_VelXnDiff(1:nx,1:ny,1:nz) )

    ! VelY
    !
    xqz_DVelYDtNl = xqz_DVelYDtNl + ( xqz_VelYNDiff + xqz_DVelYDtAdv )

    call HistoryAutoPut(TimeN, 'DVelYDtAdv',  xqz_DVelYDtAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'DVelYDtDiff', xqz_VelYnDiff(1:nx,1:ny,1:nz) )

    ! VelZ
    !
    call BuoyancyLong_xyr       ! 豬ｮ蜉幃��

    xyr_DVelZDtNl =                            &
      &   xyr_DVelZDtNl                        &
      & + (                                    &
      &      xyr_VelZnDiff                     &
      &    + xyr_DVelZDtAdv                       &
      &    + (                                 &
      &        + xyr_BuoyM * FactorBuoyMolWt   &
      &        + xyr_BuoyD * FactorBuoyLoading &
      &        + xyr_BuoyT * FactorBuoyTemp    &
      &      )                                 &
      &   )

    call HistoryAutoPut(TimeN, 'DVelZDtAdv',   xyr_DVelZDtAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'DVelZDtDiff',  xyr_VelZnDiff(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'DVelZDtBuoyT', xyr_BuoyT(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'DVelZDtBuoyM', xyr_BuoyM(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'DVelZDtBuoyD', xyr_BuoyD(1:nx,1:ny,1:nz) )

  contains

    subroutine BuoyancyLong_xyr
      !
      ! 豬ｮ蜉幃���險育ｮ�
      !

      use composition,  only: GasNum,       &!
        &                     IdxG,         &!
        &                     MolWtWet       ! 貉ｿ貎､謌仙�縺ｮ蛻�ｭ宣㍼
      use constants,   only : MolWtDry,     &! 荵ｾ辯･謌仙�縺ｮ蛻�ｭ宣㍼
        &                     Grav           ! 驥榊鴨蜉�騾溷ｺｦ
      use basicset,    only : xyz_PTempBZ,  &! 蝓ｺ譛ｬ蝣ｴ縺ｮ貂ｩ菴�
        &                     xyr_QMixBZ,   &! 蝓ｺ譛ｬ蝣ｴ縺ｮ豺ｷ蜷域ｯ�
        &                     xyr_QMixBZPerMolWt

      ! 證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
      !
      implicit none

      ! 菴懈･ｭ螟画焚
      !
      real(DP)    :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, 1:GasNum)
      real(DP)    :: xyz_QMixPerMolWtSum(imin:imax,jmin:jmax,kmin:kmax)
      real(DP)    :: xyz_QMixNlSumGas(imin:imax,jmin:jmax,kmin:kmax)
      real(DP)    :: xyz_QMixNlSumCnd(imin:imax,jmin:jmax,kmin:kmax)
      real(DP)    :: xyz_BuoyT(imin:imax,jmin:jmax,kmin:kmax)
      integer     :: i, j, k, f, n

      ! Buoyancy due to temperature disturbunce
      !
      xyz_BuoyT = Grav * xyz_PTempNl / xyz_PTempBZ

      ! 譬ｼ蟄蝉ｽ咲ｽｮ縺ｮ螟画鋤
      !
      do k = kmin, kmax-1
        xyr_BuoyT(:,:,k) =                           &
          & (                                        &
          &   xyz_BuoyT(:,:,k+1) + xyz_BuoyT(:,:,k)  &
          & ) * 5.0d-1
      end do

      !遨ｴ蝓九ａ
      xyr_BuoyT(:,:,kmax) = 0.0d0

      ! 荵ｾ辯･縺ｮ蝣ｴ蜷医� BuoyD, BuoyM 縺ｯ髮ｶ.
      !
      if ( .NOT. FlagCalcMoist ) then

        xyr_BuoyM = 0.0d0
        xyr_BuoyD = 00d0

        ! 繧ｵ繝悶Ν繝ｼ繝√Φ繧呈栢縺代ｋ
        !
        return
      end if

      ! Buoyancy due to molecular weight
      !
      do f = 1, GasNum
        n = IdxG(f)
        xyzf_QMixPerMolWt(:,:,:,f) = xyzf_QMixNl(:,:,:,n) / MolWtWet(n)
      end do

      xyz_QMixPerMolWtSum = sum(xyzf_QMixPerMolWt, 4)
      xyz_QMixNlSumGas    = sum(xyzf_QMixNl(:,:,:,1:GasNum), 4)

      do k = kmin, kmax-1
        do j = jmin, jmax
          do i = imin, imax

            xyr_BuoyM(i,j,k) =                                         &
              & + Grav                                                 &
              &   * (                                                  &
              &         xyz_QMixPerMolWtSum(i,j,k+1)                   &
              &       + xyz_QMixPerMolWtSum(i,j,k)                     &
              &     ) * 5.0d-1                                         &
              &   / ( 1.0d0 / MolWtDry + xyr_QMixBZPerMolWt(i,j,k) )   &
              & - Grav                                                 &
              &   * (                                                  &
              &         xyz_QMixNlSumGas(i,j,k+1)                      &
              &       + xyz_QMixNlSumGas(i,j,k)                        &
              &     ) * 5.0d-1                                         &
              &   / ( 1.0d0 + xyr_QMixBZ(i,j,k) )

          end do
        end do
      end do

      ! 遨ｴ蝓九ａ
      !
      xyr_BuoyM(:,:,kmax) = 0.0d0

      ! Buoyancy due to loading
      !
      xyz_QMixNlSumCnd = sum( xyzf_QMixNl(:,:,:,GasNum+1:ncmax), 4 )

      do k = kmin, kmax-1
        do j = jmin, jmax
          do i = imin, imax

            xyr_BuoyD(i,j,k) =                       &
              & - Grav                               &
              &   * (                                &
              &         xyz_QMixNlSumCnd(i,j,k+1)    &
              &       + xyz_QMixNlSumCnd(i,j,k)      &
              &     ) * 5.0d-1                       &
              &   / ( 1.0d0 + xyr_QMixBZ(i,j,k) )

          end do
        end do
      end do

      ! 遨ｴ蝓九ａ
      !
      xyr_BuoyD(:,:,kmax) = 0.0d0

    end subroutine BuoyancyLong_xyr

    subroutine QmixFall
      !
      ! 髮ｨ邊偵�關ｽ荳九↓繧医ｋ遘ｻ豬√ｒ豎ゅａ繧�.
      !

      ! 繝｢繧ｸ繝･繝ｼ繝ｫ蜻ｼ縺ｳ蜃ｺ縺�
      !
      use average,          only : xyr_xyz
      use differentiate_center4, &
        &                   only : xyz_dz_xyr
      use constants,        only : FactorJ
      use composition,      only : IdxR, RainNum
      use basicset,         only : xyzf_QMixBZ,  &!蝓ｺ譛ｬ蝣ｴ縺ｮ豺ｷ蜷域ｯ�
        &                          xyz_DensBZ     !蝓ｺ譛ｬ蝣ｴ縺ｮ蟇�ｺｦ

      !證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
      !
      implicit none

      !螟画焚螳夂ｾｩ
      !
      real(DP)  :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax,ncmax)
                                                 !闥ｸ豌玲ｷｷ蜷域ｯ�(謫ｾ荵ｱ + 蟷ｳ蝮��ｴ)
      real(DP)  :: xyrf_QMixFallFlux(imin:imax,jmin:jmax,kmin:kmax,ncmax)
                                                 !髮ｨ邊偵�關ｽ荳句柑譫�
      real(DP)  :: xyz_VelZFall(imin:imax,jmin:jmax,kmin:kmax)
                                                 !髮ｨ邊定誠荳矩溷ｺｦ
      integer   :: s, iR

      ! 蛻晄悄蛹�
      !
      xyzf_QMixAll = max( 0.0d0, xyzf_QMixBl + xyzf_QMixBZ )
      xyrf_QMixFallFlux = 0.0d0
      xyzf_QMixFall = 0.0d0
      xyz_VelZFall = 0.0d0

      ! 關ｽ荳矩���險育ｮ�
      !
      do s = 1, RainNum

        iR = IdxR(s)

        !髮ｨ邊堤ｵらｫｯ騾溷ｺｦ
        xyz_VelZFall = - 12.2d0 * FactorJ * ( xyzf_QMixAll(:,:,:,iR) ** 0.125d0 )

        ! 繝輔Λ繝�け繧ｹ縺ｮ險育ｮ�
        !
        xyrf_QMixFallFlux(:,:,:,iR) =                                &
          &  xyr_xyz (                                               &
          &    xyz_DensBZ * xyzf_QMixAll(:,:,:,iR) * xyz_VelZFall    &
          &  )

        ! 荳顔ｫｯ縺ｮ繝輔Λ繝�け繧ｹ縺ｯ繧ｼ繝ｭ
        !
        xyrf_QMixFallFlux(:,:,nz,iR) = 0.0d0

        ! 髮ｨ邊定誠荳九↓繧医ｋ譎る俣螟牙喧
        !
        xyzf_QMixFall(:,:,:,iR) =                                      &
          &  - xyz_dz_xyr( xyrf_QMixFallFlux(:,:,:,iR) ) / xyz_DensBZ

        call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(iR))//'DtFall', &
          & xyzf_QMixFall(1:nx,1:ny,1:nz,iR))

      end do

    end subroutine QmixFall

  end subroutine Dynamics_Long_forcing

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

  subroutine Dynamics_Short_forcing(  &
        &  pyz_VelXNs,          & ! (in)
        &  xqz_VelYNs,          & ! (in)
        &  xyr_VelZNs,          & ! (in)
        &  xyz_ExnerNs,         & ! (in)
        &  pyz_DVelXDtNl,       & ! (in)
        &  xqz_DVelYDtNl,       & ! (in)
        &  xyr_DVelZDtNl,       & ! (in)
        &  xyz_DExnerDtNl,      & ! (in)
        &  xyz_DExnerDtNs,      & ! (in)
        &  pyz_VelXAs,          & ! (out)
        &  xqz_VelYAs,          & ! (out)
        &  xyr_VelZAs,          & ! (out)
        &  xyz_ExnerAs          & ! (out)
        & )
    !
    ! 遏ｭ縺�凾髢薙せ繝�ャ繝励〒隧穂ｾ｡縺吶ｋ鬆��險育ｮ�.
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ隱ｭ縺ｿ霎ｼ縺ｿ
    !
    use dc_types,  only : DP
    use gridset,   only : imin, imax, jmin, jmax, kmin, kmax, &
                                             !驟榊�繧ｵ繧､繧ｺ
      &                  nx, ny, nz          !迚ｩ逅��伜沺縺ｮ螟ｧ縺阪＆
    use timeset,  only : DelTimeShort, TimeN
    use setmargin,only : SetMargin_xyzf, SetMargin_xyz, &
      &                  SetMargin_pyz, SetMargin_xqz, SetMargin_xyr
    use gtool_historyauto, &
      &           only : HistoryAutoPut
    use acousticmode_std, &
      &           only : acousticmode_std_exp, acousticmode_std_imp
    use acousticmode_2D, &
      &           only : acousticmode_2D_exp, acousticmode_2D_imp
    use acousticmode_3D, &
      &           only : acousticmode_3D_exp, acousticmode_3D_imp
    use constants,only : CvDry,    &! 荵ｾ辯･謌仙�縺ｮ豈皮�
      &                  GasRDry

    ! 證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
    !
    implicit none

    ! 螟画焚縺ｮ螳夂ｾｩ
    !
    real(DP), intent(in)    :: pyz_VelXNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xqz_VelYNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyr_VelZNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyz_ExnerNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xyz_DExnerDtNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out)   :: pyz_VelXAs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out)   :: xqz_VelYAs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out)   :: xyr_VelZAs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out)   :: xyz_ExnerAs(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: pyz_DVelXDtNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_DVelYDtNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyr_DVelZDtNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_VelDivNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_Expnd(imin:imax,jmin:jmax,kmin:kmax)

    ! subroutine 縺ｮ蠑墓焚縺ｨ縺ｪ繧矩�蛻励�, 迚ｩ逅��伜沺縺ｮ螟ｧ縺阪＆縺ｫ.
    !
    real(DP) :: pyz_PGrad(1:nx,1:ny,1:nz)
    real(DP) :: xqz_PGrad(1:nx,1:ny,1:nz)
    real(DP) :: xyr_PGrad(1:nx,1:ny,1:nz)
    real(DP) :: pyz_SWF(1:nx,1:ny,1:nz)
    real(DP) :: xqz_SWF(1:nx,1:ny,1:nz)
    real(DP) :: xyr_SWF(1:nx,1:ny,1:nz)

    !------------------------------------------------------------
    ! 蛻晄悄蛹�
    !
    pyz_DVelXDtNs = 0.0d0
    xqz_DVelYDtNs = 0.0d0
    xyr_DVelZDtNs = 0.0d0

    !------------------------------------------------------------
    ! 豌ｴ蟷ｳ譁ｹ蜷�: 髯ｽ隗｣豕�
    !
    select case ( IDAcousticmode )

    case ( IDAcousticmode_std )
      call acousticmode_std_exp(              &
        & pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, & !(IN)
        & xyz_ExnerNs,                        & !(IN)
        & xyz_VelDivNs,                       & !(OUT)
        & pyz_PGrad, pyz_SWF,                 & !(OUT)
        & xqz_PGrad, xqz_SWF                  & !(OUT)
        & )

    case ( IDAcousticmode_3D )
      call acousticmode_3d_exp(               &
        & pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, & !(IN)
        & xyz_ExnerNs,                        & !(IN)
        & xyz_VelDivNs,                       & !(IN)
        & pyz_PGrad, pyz_SWF,                 & !(OUT)
        & xqz_PGrad, xqz_SWF                  & !(OUT)
        & )

    case ( IDAcousticmode_2D )
      xqz_PGrad = 0.0d0
      xqz_SWF   = 0.0d0

      call acousticmode_2d_exp(               &
        & pyz_VelXNs, xyr_VelZNs,             & !(IN)
        & xyz_ExnerNs,                        & !(IN)
        & xyz_VelDivNs,                       & !(IN)
        & pyz_PGrad, pyz_SWF                  & !(OUT)
        & )

    end select

    ! tendency
    !
    pyz_DVelXDtNs(1:nx,1:ny,1:nz) = pyz_PGrad + pyz_SWF
    xqz_DVelYDtNs(1:nx,1:ny,1:nz) = xqz_PGrad + xqz_SWF

    ! 蛟､縺ｮ菫晉ｮ｡
    !
!    call HistoryAutoPut(TimeN, 'VelDiv', xyz_VelDivNs(1:nx,1:ny,1:nz))
!    call HistoryAutoPut(TimeN, 'VelDiv', pyz_VelXNs(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelDiv', xyr_VelZNs(1:nx,1:ny,1:nz))

    call HistoryAutoPut(TimeN, 'DVelXDtPGrad', pyz_PGrad )
    call HistoryAutoPut(TimeN, 'DVelXDtSWF',   pyz_SWF   )
    call HistoryAutoPut(TimeN, 'DVelYDtPGrad', xqz_PGrad )
    call HistoryAutoPut(TimeN, 'DVelYDtSWF',   xqz_SWF   )

    ! Time integration
    !
    pyz_VelXAs = pyz_VelXNs + DelTimeShort * ( pyz_DVelXDtNl + pyz_DVelXDtNs )
    xqz_VelYAs = xqz_VelYNs + DelTimeShort * ( xqz_DVelYDtNl + xqz_DVelYDtNs )

    ! Set Margin
    !
    call SetMargin_pyz( pyz_VelXAs ) ! (inout)
    call SetMargin_xqz( xqz_VelYAs ) ! (inout)

    !------------------------------------------------------------
    ! 驩帷峩譁ｹ蜷�: 髯ｰ隗｣豕�
    !  蜃晉ｵ宣��ｒ遏ｭ縺�凾髢薙せ繝�ャ繝励〒隧穂ｾ｡縺吶ｋ蝣ｴ蜷医ｒ閠�∴縺ｦ, 髟ｷ縺�凾髢薙せ繝�ャ繝励→
    !  遏ｭ縺�凾髢薙せ繝�ャ繝励〒隧穂ｾ｡縺輔ｌ縺� tendency 繧定ｶｳ縺怜粋繧上○縺溘ｂ縺ｮ繧貞ｼ墓焚縺ｫ.
    !

    ! 遏ｭ縺�凾髢薙せ繝�ャ繝励〒隧穂ｾ｡縺吶ｋ蝨ｧ蜉帙�蠑上� tendency
    !
    xyz_Expnd = GasRDry * xyz_ExnerNs * xyz_VelDivNs / CvDry   &
      &         * FactorDExnerDtExpnd

    ! tendency 縺ｮ蜷郁ｨ�
    !
    xyz_DExnerDtNs = xyz_DExnerDtNs + xyz_Expnd

    ! 蛟､縺ｮ菫晉ｮ｡
    !
    call HistoryAutoPut(TimeN, 'DExnerDtExpnd', xyz_Expnd(1:nx,1:ny,1:nz))

    select case ( IDAcousticmode )

    case ( IDAcousticmode_std )

      call acousticmode_std_imp(                             &
        & pyz_VelXAs, xqz_VelYAs, xyr_VelZNs, xyz_VelDivNs,  & !(IN)
        & xyz_ExnerNs,                                       & !(IN)
        & xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs,     & !(IN)
        & xyz_ExnerAs,                                       & !(OUT)
        & xyr_PGrad, xyr_SWF                                 & !(OUT)
        & )

    case ( IDAcousticmode_3D )

      call acousticmode_3D_imp(                              &
        & pyz_VelXAs, xqz_VelYAs, xyr_VelZNs, xyz_VelDivNs,  & !(IN)
        & xyz_ExnerNs,                                       & !(IN)
        & xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs,     & !(IN)
        & xyz_ExnerAs,                                       & !(OUT)
        & xyr_PGrad, xyr_SWF                                 & !(OUT)
        & )

    case ( IDAcousticmode_2D )

      call acousticmode_2D_imp(                              &
        & pyz_VelXAs, xyr_VelZNs, xyz_VelDivNs,              & !(IN)
        & xyz_ExnerNs,                                       & !(IN)
        & xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs,     & !(IN)
        & xyz_ExnerAs,                                       & !(OUT)
        & xyr_PGrad, xyr_SWF                                 & !(OUT)
        & )

    end select

    ! tendency
    !
    xyr_DVelZDtNs(1:nx,1:ny,1:nz) = xyr_PGrad + xyr_SWF

    ! Time integration
    !
    xyr_VelZAs = xyr_VelZNs + DelTimeShort * ( xyr_DVelZDtNl + xyr_DVelZDtNs )

    ! Set Margin
    !
    call SetMargin_xyz( xyz_ExnerAs ) ! (inout)
    call SetMargin_xyr( xyr_VelZAs )  ! (inout)

    ! 蛟､縺ｮ菫晉ｮ｡
    !
    call HistoryAutoPut(TimeN, 'DVelZDtPGrad',  xyr_PGrad )
    call HistoryAutoPut(TimeN, 'DVelZDtSWF',    xyr_SWF   )

  end subroutine Dynamics_Short_forcing

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

  subroutine Dynamics_Tendency_Output
    !
    ! 蜃ｺ蜉帙�險ｭ螳�
    !

    use gtool_historyauto, only : HistoryAutoAddVariable
    use gridset,           only : ncmax             ! 迚ｩ雉ｪ謨ｰ
    use composition,       only : SpcWetSymbol, &
      &                           IdxR, RainNum

    implicit none

    integer :: f, iR

    call HistoryAutoAddVariable(                              &
      & varname='DPTempDtAdv',                                   &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Advection term of potential temperature',   &
      & units='K.s-1',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                                      &
      & varname='DPTempDtDiff',                                          &
      & dims=(/'x','y','z','t'/),                                     &
      & longname='Numerical diffusion term of potential temperature', &
      & units='K.s-1',                                                &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DExnerDtAdv',                                   &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Advection term of exner function',          &
      & units='s-1',                                          &
      & xtype='double')

    call HistoryAutoAddVariable(                               &
      & varname='DExnerDtDiff',                                   &
      & dims=(/'x','y','z','t'/),                              &
      & longname='Numerical diffusion term of exner function', &
      & units='s-1',                                           &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DExnerDtExpnd',                                 &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Expanding term of exner function',          &
      & units='s-1',                                          &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='CDensAdv',                                   &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Advection term of cloud density',           &
      & units='kg.m-3.s-1',                                   &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='CDensDiff',                                  &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Numerical diffusion term of cloud density', &
      & units='kg.m-3.s-1',                                   &
      & xtype='double')

    do f = 1, ncmax
      call HistoryAutoAddVariable(                            &
        & varname='D'//trim(SpcWetSymbol(f))//'DtAdv',              &
        & dims=(/'x','y','z','t'/),                           &
        & longname='Advection term of '                       &
        &           //trim(SpcWetSymbol(f))//' mixing ratio', &
        & units='kg.kg-1.s-1',                                &
        & xtype='double')

      call HistoryAutoAddVariable(                            &
        & varname='D'//trim(SpcWetSymbol(f))//'DtDiff',             &
        & dims=(/'x','y','z','t'/),                           &
        & longname='Diffusion term of '                       &
        &           //trim(SpcWetSymbol(f))//' mixing ratio', &
        & units='kg.kg-1.s-1',                                &
        & xtype='double')
    end do

    do f = 1, RainNum
      iR = IdxR(f)
      call HistoryAutoAddVariable(                             &
        & varname='D'//trim(SpcWetSymbol(iR))//'DtFall',             &
        & dims=(/'x','y','z','t'/),                            &
        & longname='Fall term of '                             &
        &           //trim(SpcWetSymbol(iR))//' mixing ratio', &
        & units='kg.kg-1.s-1',                                 &
        & xtype='double')

    end do

    call HistoryAutoAddVariable(                              &
      & varname='DVelXDtAdv',                                    &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Advection term of velocity (x)',            &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelXDtDiff',                                   &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Numerical diffusion term of velocity (x)',  &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelXDtPGrad',                                  &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Pressure gradient term of velocity (x)',    &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelXDtSWF',                                    &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Filter for acoustic mode (x)',              &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelYDtAdv',                                    &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Advection term of velocity (y)',            &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelYDtDiff',                                   &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Numerical diffusion term of velocity (y)',  &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelYDtPGrad',                                  &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Pressure gradient term of velocity (y)',    &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelYDtSWF',                                    &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Filter for acoustic mode (y)',              &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelZDtAdv',                                    &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Advection term of velocity (z)',            &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelZDtDiff',                                   &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Numerical diffusion term of Velocity (z)',  &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelZDtBuoyT',                                  &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Buoyancy (Temperature)',                    &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelZDtBuoyM',                                  &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Buoyancy (MolWt)',                          &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelZDtBuoyD',                                  &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Buoyancy (Drag)',                           &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelZDtPGrad',                                  &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Pressure gradient term of velocity (z)',    &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DVelZDtSWF',                                    &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Filter for acoustic mode (z)',              &
      & units='m.s-2',                                        &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='VelDiv',                                     &
      & dims=(/'x','y','z','t'/),                             &
      & longname='velocity divergence',                       &
      & units='s-2',                                          &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DKmDtAdv',                                      &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Advection of Km',                           &
      & units='s-1',                                          &
      & xtype='double')

    call HistoryAutoAddVariable(                              &
      & varname='DKmDtDiff',                                     &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Diffusion term of Km',                      &
      & units='s-1',                                          &
      & xtype='double')

!    call HistoryAutoAddVariable(  &
!      & varname='D', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='matrix D',  &
!      & units='s-1',    &
!      & xtype='double')

!    call HistoryAutoAddVariable(  &
!      & varname='E', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='matrix E',  &
!      & units='s-1',    &
!      & xtype='double')

!    call HistoryAutoAddVariable(  &
!      & varname='F', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='matrix F',  &
!      & units='s-1',    &
!      & xtype='double')

!    call HistoryAutoAddVariable(  &
!      & varname='ENs', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='DExnerDtNs',  &
!      & units='s-1',    &
!      & xtype='double')

!    call HistoryAutoAddVariable(  &
!      & varname='ENl', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='DExnerDtNl',  &
!      & units='s-1',    &
!      & xtype='double')

  end subroutine Dynamics_Tendency_Output

end module DynamicsHEVI_v3