!= 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: arare5-20150225 $
! Copyright:: Copyright (C) GFD Dennou Club, 2014. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]

module DynamicsHEVI
  !
  ! 力学コア. 
  !   タイムスプリット法を利用. 音波モードとそれ以外を別々の時間刻みで解く. 
  !   短い時間ステップの計算には HE-VI 法を利用.
  !

  !モジュール読み込み
  !
  use dc_types,   only : DP

  !暗黙の型宣言禁止
  !
  implicit none

  !属性の指定
  !
  private

  ! 変数の定義
  !
  character(*), parameter:: module_name = 'dynamics_hevi'
                                                  ! モジュールの名称.
                                                  ! Module name
  real(DP), save :: FactorBuoyTemp    = 1.0d0     !浮力 (温度の寄与) の有無
                                                  !考慮しない場合は値をゼロにする.
  real(DP), save :: FactorBuoyMolWt   = 1.0d0     !浮力 (分子量効果) の有無
                                                  !考慮しない場合は値をゼロにする.
  real(DP), save :: FactorBuoyLoading = 1.0d0     !浮力 (荷重効果) の有無
                                                  !考慮しない場合は値をゼロにする.
  real(DP), save :: FactorExnerAdv = 1.0d0        !エクスナー関数の移流の有無
                                                  !考慮しない場合は値をゼロにする.
  real(DP), save :: FactorExnerExpnd = 1.0d0      !エクスナー関数の膨張項の有無
                                                  !考慮しない場合は値をゼロにする.

  integer, save      :: 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      :: 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_integrate
  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 mpi_wrapper,only : myrank
    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, &
      & FactorExnerAdv, FactorExnerExpnd,                    &
      & 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

    !-------------------------------------------------------------------
    ! 出力
    !
    if (myrank == 0) then 
      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, "FactorExnerAdv   = %f", d=(/FactorExnerAdv/) )
      call MessageNotify( "M", module_name, "FactorExnerExpnd = %f", d=(/FactorExnerExpnd/) )

      call MessageNotify( "M", module_name, "IDAdvection = %d", i=(/IDAdvection/) )
      call MessageNotify( "M", module_name, "IDAcousticmode = %d", i=(/IDAcousticmode/) )
    end if
    
  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 : FlagCalcMoist,   &
      &                  imin,            &! x 方向の配列の下限
      &                  imax,            &! x 方向の配列の上限
      &                  jmin,            &! y 方向の配列の下限
      &                  jmax,            &! y 方向の配列の上限
      &                  kmin,            &! z 方向の配列の下限
      &                  kmax,            &! z 方向の配列の上限
      &                  nx, ny, nz,      &! 物理領域の配列の上限
      &                  ncmax             ! 物質数
    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

    ! 暗黙の型宣言禁止
    !
    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_VelXAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xqz_VelYAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyr_VelZAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_KmAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_ExnerAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_PTempAdv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyzf_QMixAdv(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_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_VelXAdv,  pyz_VelXnDiff,     & !(out)
        & xqz_VelYAdv,  xqz_VelYnDiff,     & !(out)
        & xyr_VelZAdv,  xyr_VelZnDiff,     & !(out)
        & xyz_PTempAdv, xyz_PTempNDiff,    & !(out)
        & xyz_ExnerAdv, xyz_ExnerNDiff,    & !(out)
        & xyzf_QMixAdv, xyzf_QMixNDiff,    & !(out)
        & xyz_KmAdv,    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_VelXAdv,  pyz_VelXnDiff,     & !(out)
        & xqz_VelYAdv,  xqz_VelYnDiff,     & !(out)
        & xyr_VelZAdv,  xyr_VelZnDiff,     & !(out)
        & xyz_PTempAdv, xyz_PTempNDiff,    & !(out)
        & xyz_ExnerAdv, xyz_ExnerNDiff,    & !(out)
        & xyz_KmAdv,    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_VelXAdv,  pyz_VelXnDiff,     & !(out)
        & xqz_VelYAdv,  xqz_VelYnDiff,     & !(out)
        & xyr_VelZAdv,  xyr_VelZnDiff,     & !(out)
        & xyz_PTempAdv, xyz_PTempNDiff,    & !(out)
        & xyz_ExnerAdv, xyz_ExnerNDiff,    & !(out)
        & xyz_KmAdv,    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_VelYAdv   = 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_VelXAdv,  pyz_VelXnDiff,     & !(out)
        & xyr_VelZAdv,  xyr_VelZnDiff,     & !(out)
        & xyz_PTempAdv, xyz_PTempNDiff,    & !(out)
        & xyz_ExnerAdv, xyz_ExnerNDiff,    & !(out)
        & xyz_KmAdv,    xyz_KmNDiff        & !(out)
        & )

    case ( IDAdvectionCenter4_2D ) 
      
      ! Y 方向は tendency は零. 
      !
      xqz_VelYAdv   = 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_VelXAdv,  pyz_VelXnDiff,     & !(out)
        & xyr_VelZAdv,  xyr_VelZnDiff,     & !(out)
        & xyz_PTempAdv, xyz_PTempNDiff,    & !(out)
        & xyz_ExnerAdv, xyz_ExnerNDiff,    & !(out)
        & xyz_KmAdv,    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_KmAdv + xyz_KmNDiff
    
    ! output
    !
    call HistoryAutoPut(TimeN, 'KmAdv',   xyz_KmAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'KmNDiff', xyz_KmNDiff(1:nx,1:ny,1:nz) )
    
    ! エクスナー関数
    !   移流項を考慮しない場合は, FactorExnerAdv = 0 とする. 
    !
    xyz_DExnerDtNl =                                         &
      &    xyz_DExnerDtNl                                    &
      &  + ( xyz_ExnerAdv + xyz_ExnerNDiff ) * FactorExnerAdv
    
    ! output
    !
    call HistoryAutoPut(TimeN, 'ExnerAdv',   xyz_ExnerAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'ExnerNDiff', xyz_ExnerNDiff(1:nx,1:ny,1:nz) )
    
    ! 温位
    !
    xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_PTempAdv + xyz_PTempNDiff
    
    ! output
    !
    call HistoryAutoPut(TimeN, 'PTempAdv',   xyz_PTempAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'PTempNDiff', xyz_PTempNDiff(1:nx,1:ny,1:nz) )

    ! 混合比
    !
    xyzf_DQMixDtNl = xyzf_DQMixDtNl + xyzf_QMixAdv + xyzf_QMixNDiff
    
    ! output
    !
    do f = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_Adv',   &
        & xyzf_QMixAdv(1:nx,1:ny,1:nz,f))
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_NDiff', &
        & xyzf_QMixNDiff(1:nx,1:ny,1:nz,f))
    end do
    
    ! VelX 
    !
    pyz_DVelXDtNl = pyz_DVelXDtNl + pyz_VelXAdv + pyz_VelXNDiff
    
    call HistoryAutoPut(TimeN, 'VelXAdv',   pyz_VelXAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'VelXNDiff', pyz_VelXnDiff(1:nx,1:ny,1:nz) )
    
    ! VelY
    !
    xqz_DVelYDtNl = xqz_DVelYDtNl + xqz_VelYAdv + xqz_VelYNDiff
    
    call HistoryAutoPut(TimeN, 'VelYAdv',   xqz_VelYAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'VelYNDiff', xqz_VelYnDiff(1:nx,1:ny,1:nz) )
    
    ! VelZ
    !
    call BuoyancyLong_xyr       ! 浮力項
    
    xyr_DVelZDtNl =                        &
      &   xyr_DVelZDtNl                    &
      & + xyr_VelZAdv                      &
      & + xyr_VelZnDiff                    & 
      & + xyr_BuoyT * FactorBuoyTemp       &
      & + xyr_BuoyM * FactorBuoyMolWt      &
      & + xyr_BuoyD * FactorBuoyLoading
    
    call HistoryAutoPut(TimeN, 'VelZAdv',   xyr_VelZAdv(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'VelZNDiff', xyr_VelZnDiff(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'VelZBuoyT', xyr_BuoyT(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'VelZBuoyM', xyr_BuoyM(1:nx,1:ny,1:nz) )
    call HistoryAutoPut(TimeN, 'VelZBuoyD', 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 :   xyr_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)
      integer     :: i, j, k, f, n
      
      ! Buoyancy due to temperature disturbunce
      !
      do k = kmin, kmax-1
        do j = jmin, jmax
          do i = imin, imax
            
            xyr_BuoyT(i,j,k) =              &
              & Grav                        &
              & * (                         &
              &       xyz_PTempNl(i,j,k+1)  &
              &     + xyz_PTempNl(i,j,k)    &
              &   ) * 5.0d-1                &
              & / xyr_PTempBZ(i,j,k) 
            
          end do
        end do
      end do

      ! 穴埋め
      !
      xyr_BuoyT(:,:,kmax) = 0.0d0

      if ( .NOT. FlagCalcMoist ) then 
        ! 乾燥の場合は BuoyD, BuoyM は零. 
        !
        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
    
  end subroutine Dynamics_Long_forcing
  
!!!-------------------------------------------------------------!!!

  subroutine Dynamics_Long_integrate( &
    & xyz_VarBl,  xyz_DVarDtNl,                & !(IN)
    & xyzf_VarBl, xyzf_DVarDtNl,               & !(IN)
    & xyz_VarAl,                               & !(OUT)
    & xyzf_VarAl                               & !(OUT)
    )

    ! モジュール読み込み
    !
    use dc_types, only : DP
    use gridset,  only : imin,            &! x 方向の配列の下限
      &                  imax,            &! x 方向の配列の上限
      &                  jmin,            &! y 方向の配列の下限
      &                  jmax,            &! y 方向の配列の上限
      &                  kmin,            &! z 方向の配列の下限
      &                  kmax,            &! z 方向の配列の上限
      &                  ncmax             ! 物質数
    use timeset,  only : DelTimeLong
    use setmargin,only : SetMargin_xyz,   &
      &                  SetMargin_xyzf
    use fillnegative, &
      &           only : FillNegativeQMix

    ! 暗黙の型宣言禁止
    !
    implicit none

    ! 変数の定義
    !
    real(DP), intent(in)   :: xyz_VarBl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(in)   :: xyz_DVarDtNl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(in)   :: xyzf_VarBl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP), intent(in)   :: xyzf_DVarDtNl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP), intent(out)  :: xyz_VarAl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(out)  :: xyzf_VarAl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)

    !------------------------------------------------------------------
    ! スカラー量の積分 
    ! Integration
    !
    xyz_VarAl = xyz_VarBl + (2.0d0 * DelTimeLong) * xyz_DVarDtNl

    ! Set Margin
    !
    call SetMargin_xyz(xyz_VarAl)

    !------------------------------------------------------------------
    ! 混合比の積分 
    ! Integration
    !
    xyzf_VarAl = xyzf_VarBl + (2.0d0 * DelTimeLong) * xyzf_DVarDtNl

    call SetMargin_xyzf(xyzf_VarAl)

    ! QMix については, 周囲を削って負の値を埋める
    ! 穴埋めするためには, 先に setmargin しておく必要がある. 
    !
    call FillNegativeQMix(xyzf_VarAl)

    ! のり代の値の設定. Set Margin
    ! 
    call SetMargin_xyzf(xyzf_VarAl)

  end subroutine Dynamics_Long_integrate
  
!!!---------------------------------------------------------------!!!
  
  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,            &! x 方向の配列の下限
      &                  imax,            &! x 方向の配列の上限
      &                  jmin,            &! y 方向の配列の下限
      &                  jmax,            &! y 方向の配列の上限
      &                  kmin,            &! z 方向の配列の下限
      &                  kmax,            &! z 方向の配列の上限
      &                  nx, ny, nz        ! 物理領域の大きさ
    use timeset,  only : DelTimeShort, TimeN
    use setmargin,only : SetMargin_xyzf, SetMargin_xyz, &
      &                  SetMargin_pyz, SetMargin_xqz, SetMargin_xyr
    use fillnegative, &
      &           only : FillNegativeQMix
    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,                       & !(OUT)
        & 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,                       & !(OUT)
        & 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, 'VelXPGrad', pyz_PGrad )
    call HistoryAutoPut(TimeN, 'VelXSWF',   pyz_SWF   )
    call HistoryAutoPut(TimeN, 'VelYPGrad', xqz_PGrad )
    call HistoryAutoPut(TimeN, 'VelYSWF',   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   &
      &         * FactorExnerExpnd

    ! tendency の合計
    !
    xyz_DExnerDtNs = xyz_DExnerDtNs + xyz_Expnd

    ! 値の保管
    !
    call HistoryAutoPut(TimeN, 'ExnerExpnd', 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, 'VelZPGrad',  xyr_PGrad )
    call HistoryAutoPut(TimeN, 'VelZSWF',    xyr_SWF   )
    
  end subroutine Dynamics_Short_forcing
  
!!!------------------------------------------------------------------------!!!

  subroutine Dynamics_Tendency_Output
    
    use gtool_historyauto, only : HistoryAutoAddVariable
    use gridset,           only : ncmax             ! 物質数
    use composition,       only : SpcWetSymbol

    implicit none
   
    integer :: f

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

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

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

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

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

    end do

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

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

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

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

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

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

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

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

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

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

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

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

    call HistoryAutoAddVariable(                              &
      & varname='KmAdv',                                      &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Advection of Km',                           &
      & units='s-1',                                          &
      & xtype='float')
    
    call HistoryAutoAddVariable(                              &
      & varname='KmNDiff',                                    &
      & dims=(/'x','y','z','t'/),                             &
      & longname='Diffusion term of Km',                      &
      & units='s-1',                                          &
      & xtype='float')
    
  end subroutine Dynamics_Tendency_Output
  
  
end module DynamicsHEVI

