!= Module DynamicsHEVI
!
! Authors::   杉山耕一朗(SUGIYAMA Ko-ichiro), ODAKA Masatsugu 
! Version::   $Id: dynamics_hevi_v2.test.f90,v 1.1 2013/06/15 15:31:41 sugiyama Exp $ 
! Tag Name::  $Name: arare5-20150220 $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
!モデルの力学過程を計算するために必要となる関数群を束ねたモジュール
!具体的には以下の項を計算するための関数を格納する.  
!  * 移流項
!  * 浮力項
!  * 気圧傾度力項
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * エクスナー関数の空間方向の離散化において, 2 次精度の離散化を陽に利
!    用しているため, 気圧傾度力項の計算プログラムにおいて
!    differentiate_center4 モジュールを指定することはできないので注意.
!
!== Future Plans
!

module DynamicsHEVI
  !
  !陽解法を用いた力学過程の各項の計算モジュール. 
  !具体的には以下の項を計算するための関数を格納する.  
  !  * 移流項
  !  * 浮力項
  !  * 気圧傾度力項
  !

  !モジュール読み込み
  use dc_types,   only : DP, STRING
  use dc_iounit,  only : FileOpen
  use dc_message, only : MessageNotify
  use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut 

  use axesset
  use xyz_deriv_module
  use xyz_deriv_c4_module


  use mpi_wrapper,only: myrank
  use gridset, only: &
    &                 imin,            &! x 方向の配列の下限
    &                 imax,            &! x 方向の配列の上限
    &                 jmin,            &! y 方向の配列の下限
    &                 jmax,            &! y 方向の配列の上限
    &                 kmin,            &! z 方向の配列の下限
    &                 kmax,            &! z 方向の配列の上限
    &                 nx,              &! x 方向の物理領域の上限
    &                 ny,              &! x 方向の物理領域の上限
    &                 nz,              &! y 方向の物理領域の上限
    &                 ncmax,           &! 物質数
    &                 FlagCalc3D
  use constants,only: CpDry             ! 乾燥成分の比熱
  use composition, only: SpcWetSymbol
  use timeset, only:  DelTimeShort, DelTimeLong, TimeN
!  use axesset, only:  dx, dy, dz       ! 格子間隔
  use basicset, only: xyz_VelSW,       &!基本場の音速 
    &                 xyz_VelSoundBZ,  &!基本場の音速 
    &                 xyz_DensBZ,      &!基本場の密度
    &                 xyz_PTempBZ,     &!基本場の温位
    &                 xyz_VPTempBZ,    &!基本場の温位
    &                 pyz_VPTempBZ,    &!基本場の温位
    &                 xqz_VPTempBZ,    &!基本場の温位
    &                 xyr_VPTempBZ,    &!基本場の温位
    &                 xyz_ExnerBZ,     &
    &                 xyzf_QMixBZ,     &
    &                 xyr_QMixBZPerMolWt, &
    &                 xyr_QMixBZ, xyz_EffMolWtBZ   
  use setmargin,only: SetMargin_xyzf, SetMargin_xyz, &
    &                 SetMargin_pyz, SetMargin_xqz, SetMargin_xyr
  use fillnegative, only: FillNegativeQMix
  use namelist_util, only: namelist_filename

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private
  
  real(DP), save         :: beta  = 1.0d0         !クランクニコルソン法なら 0.5
                                                  !完全陰解法なら 1

  integer, save          :: N = 10                !係数行列/改行列の次数, 整合寸法
  integer, save          :: M = 10                !方程式の組数
  integer, save          :: NUD = 1               !係数行列の上三角部分の帯幅
  integer, save          :: NLD = 1               !係数行列の下三角部分の帯幅
  integer, save          :: NAL = 1               !LU 分解の結果 L の整合寸法
  integer, save          :: NA = 3                !NUD + NLD + 1

  real(DP), allocatable, save :: xyz_F3BZ(:,:,:)            !係数行列の計算に用いる配列
  real(DP), allocatable, save :: xyz_F1BZ(:,:,:)            !係数行列の計算に用いる配列
  real(DP), allocatable, save :: xyr_F2BZ(:,:,:)            !係数行列の計算に用いる配列
  real(DP), allocatable, save :: xyr_CpVPTempBZ(:,:,:)      !係数行列の計算に用いる配列
  real(DP), allocatable, save :: xyr_CpDensVPTemp2BZ(:,:,:) !係数行列の計算に用いる配列
  real(DP), allocatable, save :: xyr_DensVPTempBZ(:,:,:)    !係数行列の計算に用いる配列

  real(DP), allocatable, save :: A(:)             !係数行列の対角成分
  real(DP), allocatable, save :: B(:)             !係数行列の上三角部分
  real(DP), allocatable, save :: C(:)             !係数行列の下三角部分
  real(DP), allocatable, save :: AL1(:)           !LU 分解の結果 L (1 次元配列)
  integer,  allocatable, save :: IP(:)            !部分ピボット交換の情報を格納

  real(DP), save :: Alpha = 0.0d0                 !音波減衰項の減衰係数
  real(DP), save :: NuHh  = 0.0d0                 !数値粘性の係数 (水平方向)
  real(DP), save :: NuVh  = 0.0d0                 !数値粘性の係数 (鉛直方向)
  real(DP), save :: NuHm  = 0.0d0                 !数値粘性の係数 (水平方向)
  real(DP), save :: NuVm  = 0.0d0                 !数値粘性の係数 (鉛直方向)
  character(*), parameter:: module_name = 'DynamicHEVI'
                                                  ! モジュールの名称.
                                                  ! Module name
  real(DP), save :: FactorBuoyTemp    = 1.0d0     !浮力 (温度の寄与) の有無
                                                  !考慮しない場合は値をゼロにする.
  real(DP), save :: FactorBuoyMolWt   = 1.0d0     !浮力 (分子量効果) の有無
                                                  !考慮しない場合は値をゼロにする.
  real(DP), save :: FactorBuoyLoading = 1.0d0     !浮力 (荷重効果) の有無
                                                  !考慮しない場合は値をゼロにする.

  !public 
  public Dynamics_Init
  public Dynamics_Long_forcing
  public Dynamics_Short_forcing
  public Dynamics_Km_forcing

contains

  subroutine Dynamics_Init

    !暗黙の型宣言禁止
    implicit none
    
    real(DP)  :: DelXMin, DelYMin, DelZMin
    real(DP)  :: AlphaSound = 5.0d-2  !音波減衰項の係数 (気象庁数値予報課報告・別冊49 より)
    real(DP)  :: AlphaNDiff = 1.0d-3  !4次の数値拡散の係数. CReSS マニュアルより
    real(DP)  :: NDiffRatio = 1.0d0   !速度に対する粘性を上げる場合は数字を 1 以上にする. 
    integer   :: unit                 !装置番号

    NAMELIST /Dynamics_nml/                                    &
         & AlphaSound, AlphaNDiff, NDiffRatio, beta,           &
         & FactorBuoyTemp, FactorBuoyMolWt, FactorBuoyLoading

    !ファイルオープン. 情報取得. 
    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=dynamics_nml)
    close(unit)
   
    !-------------------------------------------------------------------
    ! 音波減衰項の減衰率   Min(DelX, DelZ) ** 2.0 に比例
    ! 2 次元計算の場合には DelY に依存しないようにするために if 文を利用.
    !
    DelXMin = dx
    DelYMin = dy
    DelZMin = dz
    if ( FlagCalc3D ) then 
      Alpha = AlphaSound * ( Min(DelXMin, DelYMin, DelZMin) ** 2.0d0 ) / DelTimeShort
    else
      Alpha = AlphaSound * ( Min(DelXMin, DelZMin) ** 2.0d0 ) / DelTimeShort
    end if

    if ( FlagCalc3D ) then 
      NuHh = AlphaNDiff * ( SQRT(DelXMin*DelYMin) ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    else
      NuHh = AlphaNDiff * ( DelXMin ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    end if

    NuVh = AlphaNDiff * ( DelZMin ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    NuHm = NuHh * NDiffRatio
    NuVm = NuVh * NDiffRatio

    if (myrank == 0) then 
      call MessageNotify( "M", module_name, "Alpha = %f", d=(/Alpha/) )
      call MessageNotify( "M", module_name, "NuHh = %f", d=(/NuHh/) )
      call MessageNotify( "M", module_name, "NuVh = %f", d=(/NuVh/) )
      call MessageNotify( "M", module_name, "NuHm = %f", d=(/NuHm/) )
      call MessageNotify( "M", module_name, "NuVm = %f", d=(/NuVm/) )
      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/) )
    end if

    ! 鉛直陰解法を用いるために, 行列の準備を行う. 
    !
    call Dynamics_VI_init

    ! tendency の出力
    !
    call Dynamics_Tendency_output

  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)
    & xyzf_QMixBl, xyzf_QMixNl,       & ! (in)
    & pyz_DVelXDtNl,                  & ! (inout)
    & xqz_DVelYDtNl,                  & ! (inout)
    & xyr_DVelZDtNl,                  & ! (inout)
    & xyz_DPTempDtNl,                 & ! (inout)
    & xyzf_DQMixDtNl,                 & ! (inout)
    & xyz_PTempAl,                    & ! (out)
    & xyzf_QMixAl                     & ! (out)
    & )

    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)    :: 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(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) :: xyzf_DQMixDtNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(out)   :: xyz_PTempAl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(out)   :: xyzf_QMixAl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)

    real(DP)             :: pyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_nDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_nDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_nDiff(imin:imax,jmin:jmax,kmin:kmax)
    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)
    real(DP)             :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyzf_Adv(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_nDiff(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)

    integer              :: f

    real(DP)                :: xyz_Adv2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)                :: xyzf_Adv2(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP)                :: xyzf_nDiff2(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)

    real(DP)                :: pyz_Adv2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: pyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)                :: xqz_Adv2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xqz_nDiff2(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)                :: xyr_Adv2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyr_nDiff2(imin:imax,jmin:jmax,kmin:kmax)

    character(15) :: info
    integer :: k

    !--------------------------------------------------------------------
    ! 温位

    ! 移流については, 基本場も考慮する.
    !
    xyz_PTempAll = xyz_PTempNl + xyz_PTempBZ

    call AdvC4_nDiff_xyz(                    &
      &  xyz_PTempBl, xyz_PTempAll,          & !(IN)
      &  pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN)
      &  xyz_Adv, xyz_nDiff )                  !(OUT)

    xyz_Adv2 = &
      & - xyz_avr_pyz(pyz_VelXNl * pyz_c4dx_xyz(xyz_PTempAll))  &
      & - xyz_avr_xqz(xqz_VelYNl * xqz_c4dy_xyz(xyz_PTempAll))  &
      & - xyz_avr_xyr(xyr_VelZNl * xyr_c4dz_xyz(xyz_PTempAll))  

    xyz_nDiff2 = &
      &  - NuHh * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))))) &
      &  - NuHh * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))))) &
      &  - NuVh * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_PTempBl ))))) 

    info = "PTemp, Adv: "
    call CHK_Val( xyz_Adv(imin:imax,jmin:jmax,kmin:kmax), xyz_Adv2(imin:imax,jmin:jmax,kmin:kmax), info )
    info = "PTemp, nDiff: "
    call CHK_Val( xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax), xyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax), info )

    ! tendency の更新
    !
    xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_Adv + xyz_nDiff
    
    ! output
    !
    call HistoryAutoPut(TimeN, 'PTempAdv',   xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'PTempNDiff', xyz_nDiff(1:nx,1:ny,1:nz))

    
    !--------------------------------------------------------------------
    ! 混合比
    ! 

    ! 移流については, 基本場も考慮する.
    !
    xyzf_QMixAll = xyzf_QMixNl + xyzf_QMixBZ

    call AdvC4_nDiff_xyzf(                   & 
      &  xyzf_QMixBl, xyzf_QMixAll,          & !(IN)
      &  pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN)
      &  xyzf_Adv, xyzf_nDiff )                !(OUT)

    do f = 1, ncmax
      xyzf_Adv2(:,:,:,f) = &
        & - xyz_avr_pyz(pyz_VelXNl * pyz_c4dx_xyz(xyzf_QMixAll(:,:,:,f))) &
        & - xyz_avr_xqz(xqz_VelYNl * xqz_c4dy_xyz(xyzf_QMixAll(:,:,:,f))) &
        & - xyz_avr_xyr(xyr_VelZNl * xyr_c4dz_xyz(xyzf_QMixAll(:,:,:,f)))    

      xyzf_nDiff2(:,:,:,f) = &
        &  - NuHh * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))))) &
        &  - NuHh * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))))) &
        &  - NuVh * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyzf_QMixBl(:,:,:,f) )))))
    end do

    do f = 1, ncmax
      info = "QMix, Adv: "
      call CHK_Val( xyzf_Adv(imin:imax,jmin:jmax,kmin:kmax,f), xyzf_Adv2(imin:imax,jmin:jmax,kmin:kmax,f), info )
      info = "QMix, nDiff: "
      call CHK_Val( xyzf_nDiff(imin:imax,jmin:jmax,kmin:kmax,f), xyzf_nDiff2(imin:imax,jmin:jmax,kmin:kmax,f), info )
    end do

    ! tendency の更新    
    !
    xyzf_DQMixDtNl = xyzf_DQMixDtNl + xyzf_Adv + xyzf_nDiff
    
    ! output
    !
    do f = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_Adv',   xyzf_Adv(1:nx,1:ny,1:nz,f))
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_NDiff', xyzf_NDiff(1:nx,1:ny,1:nz,f))
    end do

    !------------------------------------------------------------------
    ! VelX, VelY, VelZ
    ! 

    ! 移流項・数値拡散項をまとめて計算
    !
    call AdvC4_nDiff_pyz_xqz_xyr(            & 
      &  pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, & !(IN)
      &  pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN)
      &  pyz_Adv,    xqz_Adv,    xyr_Adv,    & !(OUT)
      &  pyz_nDiff,  xqz_nDiff,  xyr_nDiff )   !(OUT)

    pyz_Adv2  = &
      & - pyz_VelXNl * pyz_avr_xyz( xyz_c4dx_pyz( pyz_VelXNl ) )                &
      & - pyz_avr_pqz( pqz_avr_xqz( xqz_VelYNl ) * pqz_c4dy_pyz( pyz_VelXNl ) ) &
      & - pyz_avr_pyr( pyr_avr_xyr( xyr_VelZNl ) * pyr_c4dz_pyz( pyz_VelXNl ) )

    pyz_nDiff2 = &
      & - NuHm * ( pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz( pyz_VelXBl ))))) &
      & - NuHm * ( pyz_dy_pqz(pqz_dy_pyz(pyz_dy_pqz(pqz_dy_pyz( pyz_VelXBl ))))) &
      & - NuVm * ( pyz_dz_pyr(pyr_dz_pyz(pyz_dz_pyr(pyr_dz_pyz( pyz_VelXBl )))))

    xqz_Adv2  = &
      & - xqz_avr_pqz( pqz_avr_pyz( pyz_VelXNl ) * pqz_c4dx_xqz( xqz_VelYNl ) ) &
      & - xqz_VelYNl * xqz_avr_xyz( xyz_c4dy_xqz( xqz_VelYNl ) )                &
      & - xqz_avr_xqr( xqr_avr_xyr( xyr_VelZNl ) * xqr_c4dz_xqz( xqz_VelYNl ) )

    xqz_nDiff2 = &
      & - NuHm * ( xqz_dx_pqz(pqz_dx_xqz(xqz_dx_pqz(pqz_dx_xqz( xqz_VelYBl ))))) &
      & - NuHm * ( xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz( xqz_VelYBl ))))) &
      & - NuVm * ( xqz_dz_xqr(xqr_dz_xqz(xqz_dz_xqr(xqr_dz_xqz( xqz_VelYBl )))))

    xyr_Adv2  = &
      & - xyr_avr_pyr( pyr_avr_pyz( pyz_VelXNl ) * pyr_c4dx_xyr( xyr_VelZNl ) ) &
      & - xyr_avr_xqr( xqr_avr_xqz( xqz_VelYNl ) * xqr_c4dy_xyr( xyr_VelZNl ) ) &
      & - xyr_VelZNl * xyr_avr_xyz( xyz_c4dz_xyr( xyr_VelZNl ) )

    xyr_nDiff2 = &
      & - NuHm * ( xyr_dx_pyr(pyr_dx_xyr(xyr_dx_pyr(pyr_dx_xyr( xyr_VelZBl ))))) &
      & - NuHm * ( xyr_dy_xqr(xqr_dy_xyr(xyr_dy_xqr(xqr_dy_xyr( xyr_VelZBl ))))) & 
      & - NuVm * ( xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr( xyr_VelZBl )))))    

    info = "VelX, Adv: "
    call CHK_Val( pyz_Adv(imin:imax,jmin:jmax,kmin:kmax), pyz_Adv2(imin:imax,jmin:jmax,kmin:kmax), info )
    info = "VelX, nDiff: "
    call CHK_Val( pyz_nDiff(imin:imax,jmin:jmax,kmin:kmax), pyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax), info )

    info = "VelY, Adv: "
    call CHK_Val( xqz_Adv(imin:imax,jmin:jmax,kmin:kmax), xqz_Adv2(imin:imax,jmin:jmax,kmin:kmax), info )
    info = "VelY, nDiff: "
    call CHK_Val( xqz_nDiff(imin:imax,jmin:jmax,kmin:kmax), xqz_nDiff2(imin:imax,jmin:jmax,kmin:kmax), info )

    info = "VelZ, Adv: "
    call CHK_Val( xyr_Adv(imin:imax,jmin:jmax,kmin:kmax), xyr_Adv2(imin:imax,jmin:jmax,kmin:kmax), info )
    info = "VelZ, nDiff: "
    call CHK_Val( xyr_nDiff(imin:imax,jmin:jmax,kmin:kmax), xyr_nDiff2(imin:imax,jmin:jmax,kmin:kmax), info )


    ! tendency of VelX
    !
    pyz_DVelXDtNl = pyz_DVelXDtNl + pyz_Adv + pyz_NDiff
    
    call HistoryAutoPut(TimeN, 'VelXAdv',   pyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXNDiff', pyz_nDiff(1:nx,1:ny,1:nz))

    ! tendency of VelY
    !
    xqz_DVelYDtNl = xqz_DVelYDtNl + xqz_Adv + xqz_NDiff
    
    call HistoryAutoPut(TimeN, 'VelYAdv',   xqz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYNDiff', xqz_nDiff(1:nx,1:ny,1:nz))

    ! Buoyancy 
    ! 
    call BuoyancyLong_xyr(                 &
      &  xyz_PTempNl, xyzf_QMixNl,         & !(IN)
      &  xyr_BuoyT, xyr_BuoyM, xyr_BuoyD )   !(OUT)

    ! tendency of VelZ
    !
    xyr_DVelZDtNl = xyr_DVelZDtNl + xyr_Adv + xyr_NDiff          &
         &          + xyr_BuoyT * FactorBuoyTemp                 &
         &          + xyr_BuoyM * FactorBuoyMolWt                &
         &          + xyr_BuoyD * FactorBuoyLoading

    call HistoryAutoPut(TimeN, 'VelZAdv',   xyr_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelZNDiff', xyr_NDiff(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))    


    !------------------------------------------------------------------
    ! 長い時間ステップでの積分 
    ! Integration
    !
    xyz_PTempAl = xyz_PTempBl + (2.0d0 * DelTimeLong) * xyz_DPTempDtNl
    xyzf_QMixAl = xyzf_QMixBl + (2.0d0 * DelTimeLong) * xyzf_DQMixDtNl

    ! Set Margin
    !
    call SetMargin_xyz(xyz_PTempAl)
    call SetMargin_xyzf(xyzf_QMixAl)

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

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

  end subroutine Dynamics_Long_forcing

  
  subroutine Dynamics_Long_Mars_forcing( &
    & pyz_VelXBl,  pyz_VelXNl,        & ! (in)
    & xqz_VelYBl,  xqz_VelYNl,        & ! (in)
    & xyr_VelZBl,  xyr_VelZNl,        & ! (in)
    & xyz_PTempBl, xyz_PTempNl,       & ! (in)
    & xyz_CDensBl, xyz_CDensNl,       & ! (in)
    & pyz_DVelXDtNl,                  & ! (inout)
    & xqz_DVelYDtNl,                  & ! (inout)
    & xyr_DVelZDtNl,                  & ! (inout)
    & xyz_DPTempDtNl,                 & ! (inout)
    & xyz_DCDensDtNl,                 & ! (inout)
    & xyz_PTempAl                     & ! (out)
    & )

    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_CDensBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyz_CDensNl(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_DCDensDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out), optional :: xyz_PTempAl(imin:imax,jmin:jmax,kmin:kmax) 

    real(DP)             :: pyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_nDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_nDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_nDiff(imin:imax,jmin:jmax,kmin:kmax)
    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)
    real(DP)             :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax)

    !--------------------------------------------------------------------
    ! 温位

    ! 移流については, 基本場も考慮する.
    !
    xyz_PTempAll = xyz_PTempNl + xyz_PTempBZ

    call AdvC4_nDiff_xyz(                    &
      &  xyz_PTempBl, xyz_PTempAll,          & !(IN)
      &  pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN)
      &  xyz_Adv, xyz_nDiff )                  !(OUT)
    
    ! tendency の更新
    !   基本場の移流の分を加えている. 
    !
    xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_Adv + xyz_nDiff
    
    ! output
    !
    call HistoryAutoPut(TimeN, 'PTempAdv',   xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'PTempNDiff', xyz_nDiff(1:nx,1:ny,1:nz))


    !--------------------------------------------------------------------
    ! 雲密度
    !
    call AdvC4_nDiff_xyz(                    &
      &  xyz_CDensBl, xyz_CDensNl,           & !(IN)
      &  pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN)
      &  xyz_Adv, xyz_nDiff )                  !(OUT)
    
    ! tendency の更新
    !
    xyz_DCDensDtNl = xyz_DCDensDtNl + xyz_Adv + xyz_nDiff

    ! 値の保管
    !
    call HistoryAutoPut(TimeN, 'CDensAdv',   xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'CDensNDiff', xyz_nDiff(1:nx,1:ny,1:nz))
    
    !------------------------------------------------------------------
    ! VelX, VelY, VelZ
    ! 

    ! 移流項・数値拡散項をまとめて計算
    !
    call AdvC4_nDiff_pyz_xqz_xyr(            & 
      &  pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, & !(IN)
      &  pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN)
      &  pyz_Adv,    xqz_Adv,    xyr_Adv,    & !(OUT)
      &  pyz_nDiff,  xqz_nDiff,  xyr_nDiff )   !(OUT)
    
    ! tendency of VelX
    !
    pyz_DVelXDtNl = pyz_DVelXDtNl + pyz_Adv + pyz_NDiff
    
    call HistoryAutoPut(TimeN, 'VelXAdv',   pyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXNDiff', pyz_nDiff(1:nx,1:ny,1:nz))

    ! tendency of VelY
    !
    xqz_DVelYDtNl = xqz_DVelYDtNl + xqz_Adv + xqz_NDiff
    
    call HistoryAutoPut(TimeN, 'VelYAdv',   xqz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYNDiff', xqz_nDiff(1:nx,1:ny,1:nz))

    ! Buoyancy 
    !  火星版の浮力は? 要チェック. CHECK!!
    ! 
!    call BuoyancyLong_xyr(                 &
!      &  xyz_PTempNl, xyzf_QMixNl,         & !(IN)
!      &  xyr_BuoyT, xyr_BuoyM, xyr_BuoyD )   !(OUT)
    xyr_BuoyT = 0.0d0; xyr_BuoyM = 0.0d0; xyr_BuoyD = 0.0d0
    
    ! tendency of VelZ
    !
    xyr_DVelZDtNl = xyr_DVelZDtNl + xyr_Adv + xyr_NDiff          &
         &          + xyr_BuoyT * FactorBuoyTemp                 &
         &          + xyr_BuoyM * FactorBuoyMolWt                &
         &          + xyr_BuoyD * FactorBuoyLoading

    call HistoryAutoPut(TimeN, 'VelZAdv',   xyr_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelZNDiff', xyr_NDiff(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))    


    !------------------------------------------------------------------
    ! 長い時間ステップでの積分
    ! Integration
    !
    if (present(xyz_PTempAl)) then 
      ! 積分
      !
      xyz_PTempAl = xyz_PTempBl + (2.0d0 * DelTimeLong) * xyz_DPTempDtNl

      ! Set Margin
      !
      call SetMargin_xyz(xyz_PTempAl)
    end if

  end subroutine Dynamics_Long_Mars_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)
        & )

    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(inout)  :: xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax) !test
    real(DP), intent(inout)  :: xyz_DExnerDtNs(imin:imax,jmin:jmax,kmin:kmax) !test
    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) :: pyz_PGrad(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_PGrad(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyr_PGrad(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: pyz_SWF(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_SWF(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyr_SWF(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: pyz_PGrad2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_PGrad2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyr_PGrad2(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: pyz_SWF2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_SWF2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyr_SWF2(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: xyz_VelDivNs2(imin:imax,jmin:jmax,kmin:kmax)

    integer  :: k
    character(15) :: info

    ! initialize: Divergence of velocity
    !
    call VelDivC2(                           &
      &  pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, & !(IN)
      &  xyz_VelDivNs )                        !(OUT)
    
    xyz_VelDivNs2 = &
      &   xyz_dx_pyz( pyz_VelXNs ) &
      & + xyz_dy_xqz( xqz_VelYNs ) &
      & + xyz_dz_xyr( xyr_VelZNs )

    info = "Vel, Div: "
    call CHK_Val( xyz_VelDivNs(imin:imax,jmin:jmax,kmin:kmax), xyz_VelDivNs2(imin:imax,jmin:jmax,kmin:kmax), info )

    !------------------------------------------------------------
    ! VelX, VelY
    !  水平方向は陽解法で解く. 
    !
    call PGrad_HE(                  &
      &  xyz_VelDivNs, xyz_ExnerNs, & !(IN)
      &  pyz_PGrad, pyz_SWF,        & !(OUT)
      &  xqz_PGrad, xqz_SWF )         !(OUT)

    pyz_PGrad2 = - pyz_avr_xyz( CpDry * xyz_VPTempBZ ) * pyz_dx_xyz( xyz_ExnerNs ) 
    pyz_SWF2   =   Alpha * pyz_dx_xyz( xyz_VelDivNs ) 
    xqz_PGrad2 = - xqz_avr_xyz( CpDry * xyz_VPTempBZ ) * xqz_dy_xyz( xyz_ExnerNs ) 
    xqz_SWF2   =   Alpha * xqz_dy_xyz( xyz_VelDivNs ) 

    info = "VelX, PGrad: "
    call CHK_Val( pyz_PGrad(imin:imax,jmin:jmax,kmin:kmax), pyz_PGrad2(imin:imax,jmin:jmax,kmin:kmax), info )
    info = "VelX, SWF: "
    call CHK_Val( pyz_SWF(imin:imax,jmin:jmax,kmin:kmax), pyz_SWF2(imin:imax,jmin:jmax,kmin:kmax), info )

    info = "VelY, PGrad: "
    call CHK_Val( xqz_PGrad(imin:imax,jmin:jmax,kmin:kmax), xqz_PGrad2(imin:imax,jmin:jmax,kmin:kmax), info )
    info = "VelY, SWF: "
    call CHK_Val( xqz_SWF(imin:imax,jmin:jmax,kmin:kmax), xqz_SWF2(imin:imax,jmin:jmax,kmin:kmax), info )

    ! tendency 
    !
    pyz_DVelXDtNs = pyz_PGrad + pyz_SWF
    xqz_DVelYDtNs = xqz_PGrad + xqz_SWF

    ! 値の保管
    !
    call HistoryAutoPut(TimeN, 'VelXPGrad', pyz_PGrad(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXSWF',   pyz_SWF(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYPGrad', xqz_PGrad(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYSWF',   xqz_SWF(1:nx,1:ny,1:nz))

    ! 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)
    
    !------------------------------------------------------------
    ! Exner function
    !  積分値を返すことに注意. 
    !
    call Exner_HEVI(                               & 
      &  pyz_VelXAs,  xqz_VelYAs,  xyr_VelZNs,     & !(IN)
      &  xyz_VelDivNs, xyz_ExnerNs,                & !(IN)
      &  xyr_DVelZDtNl, xyz_DExnerDtNl,            & !(IN)
      &  xyz_ExnerAs                               & !(OUT)
      & )

    !------------------------------------------------------------
    ! VelZ
    !
    call PGrad_VI(                               &
      &  Alpha, beta, dz, CpDry,  xyr_VPTempBZ,  &
      &  imin, imax, jmin, jmax, kmin, kmax,     &  
      &  xyz_VelDivNs, xyz_ExnerNs, xyz_ExnerAs, &!(IN)
      &  xyr_PGrad, xyr_SWF )                     !(OUT)

    xyr_SWF2 =  Alpha * xyr_dz_xyz( xyz_VelDivNs ) 
    xyr_PGrad2 =                                             &
      & - xyr_avr_xyz(CpDry * xyz_VPTempBZ )                 &
      &   * (                                                &
      &       beta * xyr_dz_xyz( xyz_ExnerAs )               &
      &       + (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs )   &
      &     )

    info = "VelZ, PGrad: "
    call CHK_Val( xyr_PGrad(imin:imax,jmin:jmax,kmin:kmax), xyr_PGrad2(imin:imax,jmin:jmax,kmin:kmax), info )
    info = "VelZ, SWF: "
    call CHK_Val( xyr_SWF(imin:imax,jmin:jmax,kmin:kmax), xyr_SWF2(imin:imax,jmin:jmax,kmin:kmax), info )

    ! tendency
    !
    xyr_DVelZDtNs = xyr_PGrad + xyr_SWF

    ! 値の保管
    !
    call HistoryAutoPut(TimeN, 'VelZPGrad', xyr_PGrad(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelZSWF',   xyr_SWF(1:nx,1:ny,1:nz))

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

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

  end subroutine Dynamics_Short_forcing

  
!!!--------------------------------------------------------------------!!!
  subroutine Dynamics_VI_init()
    !
    !エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素を決め, 
    !LU 分解を行う. 
    !

    !暗黙の型宣言禁止
    implicit none

    real(DP)  :: dt      ! 短い時間格子
    integer   :: INFO  !解のコンディションチェック
    integer   :: i, j, k

    real(8) :: A2(1:nz) 
    real(8) :: B2(2:nz) 
    real(8) :: C2(1:nz-1) 
    real(8) :: xyz_F1BZ2(imin:imax, jmin:jmax, kmin:kmax)
    real(8) :: xyr_F2BZ2(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)  :: dts      ! 短い時間格子

    character(15) :: info2

    !----------------------------------------------------------------
    ! 初期化

    ! 変数名が長すぎたので, 名前を置き換える
    !
    dt = DelTimeShort
    dts = DelTimeShort

    ! 配列の割り付け
    !
    allocate( A(1:nz) )
    allocate( B(2:nz) )
    allocate( C(1:nz-1) )
    allocate( xyz_F3BZ(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_F1BZ(1:nx,1:ny,1:nz) )
    allocate( xyr_F2BZ(1:nx,1:ny,kmin:kmax) )
    allocate( xyr_CpDensVPTemp2BZ(1:nx,1:ny,kmin:kmax) )
    allocate( xyr_DensVPTempBZ(1:nx,1:ny,kmin:kmax) )
    allocate( xyr_CpVPTempBZ(1:nx,1:ny,kmin:kmax) )

    !----------------------------------------------------------------
    ! 係数行列および共通して利用される配列の値を決める
    !   A, B, C を求める際, 基本場 (BZ) の量は X 方向に一様なので. 
    !   nx, ny の値で代表させることとした. 
    !
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx
          xyz_F1BZ(i,j,k) =                                                    &
            &  ( xyz_VelSW(i,j,k) ** 2.0d0 )                                   &
            &  / ( CpDry * xyz_DensBZ(i,j,k) * (xyz_VPTempBZ(i,j,k) ** 2.0d0) )
        end do
      end do
    end do

    xyz_F1BZ2 =                                               &
      &  ( xyz_VelSoundBZ ** 2.0d0 )                          &
      &   / (CpDry * xyz_DensBZ * (xyz_VPTempBZ ** 2.0d0))    

    xyz_F3BZ = xyz_F1BZ2

    info2 = "VI_init, F1BZ: "
    call CHK2_Val( xyz_F1BZ(1:nx,1:ny,1:nz), xyz_F1BZ2(1:nx,1:ny,1:nz), info2 )

    do k = kmin, kmax-1
      do j = 1, ny        
        do i = 1, nx
          xyr_CpDensVPTemp2BZ(i,j,k)=                    &
            &  CpDry                                     &
            &  * (                                       &
            &    + xyz_DensBZ(i,j,k+1)                   &
            &      * ( xyz_VPTempBZ(i,j,k+1) ** 2.0d0 )  &
            &    + xyz_DensBZ(i,j,k)                     &
            &      * ( xyz_VPTempBZ(i,j,k) ** 2.0d0 )    &
            &    ) * 5.0d-1
        end do
      end do
    end do
    xyr_CpDensVPTemp2BZ(:,:,kmax) = 0.0d0  !穴埋め
    
    xyr_F2BZ2 =                                                &
      &  xyr_avr_xyz(                                         &
      &    CpDry * xyz_DensBZ * ( xyz_VPTempBZ ** 2.0d0 )     &
      &   )

    info2 = "VI_init, xyr_CpDensVPTemp2BZ: "
    call CHK2_Val( xyr_CpDensVPTemp2BZ(1:nx,1:ny,1:nz) , xyr_F2BZ2(1:nx,1:ny,1:nz) , info2 )

    xyr_F2BZ = xyr_CpDensVPTemp2BZ

    do k = kmin, kmax-1
      do j = 1, ny
        do i = 1, nx
          xyr_DensVPTempBZ(i,j,k) =                             &
            & + (                                               &
            &   + xyz_DensBZ(i,j,k+1) * xyz_VPTempBZ(i,j,k+1)   &
            &   + xyz_DensBZ(i,j,k)   * xyz_VPTempBZ(i,j,k)     &
            &   ) * 5.0d-1    
        end do
      end do
    end do
    xyr_DensVPTempBZ(:,:,kmax) = 0.0d0  !穴埋め

    do k = kmin, kmax-1
      do j = 1, ny
        do i = 1, nx
          xyr_CpVPTempBZ(i,j,k) =          &
            &   CpDry                      &
            &   * (                        &
            &     + xyz_VPTempBZ(i,j,k+1)  &
            &     + xyz_VPTempBZ(i,j,k)    &
            &     ) * 5.0d-1
        end do
      end do
    end do
    xyr_CpVPTempBZ(:,:,kmax) = 0.0d0
          
    do k = 2, nz-1
      A(k) =                                        &
        & + 1.0d0                                   &
        & + ( beta ** 2.0d0 )                       &
        &    * xyz_F1BZ(nx,ny,k) * ( dt * dt )      &
        &    * (                                    &
        &         xyr_CpDensVPTemp2BZ(nx,ny,k)      &
        &       + xyr_CpDensVPTemp2BZ(nx,ny,k-1)    &
        &       )                                   &
        &    / ( dz * dz )
    end do

    A(1) =                                   &
      & + 1.0d0                              &
      & + ( beta ** 2.0d0 )                  &
      &   * xyz_F1BZ(nx,ny,1) * ( dt * dt )  &
      &   * xyr_CpDensVPTemp2BZ(nx,ny,1)     &
      &   / ( dz * dz ) 

    A(nz) =                                  &
      & + 1.0d0                              &
      & + ( beta ** 2.0d0 )                  &
      &   * xyz_F1BZ(nx,ny,nz) * ( dt * dt ) &
      &   * xyr_CpDensVPTemp2BZ(nx,ny,nz-1)  &
      &   / ( dz * dz )  

    do k = 2, nz
      B(k) =                                     &
        & - ( beta ** 2.0d0 )                    &
        &   * xyz_F1BZ(nx,ny,k-1) * ( dt * dt )  &
        &   * xyr_CpDensVPTemp2BZ(nx,ny,k-1)     &
        &   / ( dz * dz )
    end do

    do k = 1, nz-1
      C(k) =                                     &
        & - ( beta ** 2.0d0 )                    &
        &   * xyz_F1BZ(nx,ny,k+1) * (dt * dt )   &
        &   * xyr_CpDensVPTemp2BZ(nx,ny,k)       &
        &   / ( dz * dz )
    end do

    A2(1+1: nz-1) =                                &
      & (beta ** 2.0d0)                                      &
      &    * xyz_F1BZ2(nx,ny,1+1: nz-1)  &  
      &    * (DTS ** 2.0d0)                                  &
      &    * (                                               &
      &          xyr_F2BZ2(nx,ny,1+1: nz-1)  &
      &            / r_dz(1+1: nz-1)              &
      &        + xyr_F2BZ2(nx,ny,1  : nz-2)  &
      &            / r_dz(1: nz-2)                &
      &       )                                              &
      &    / z_dz(1+1: nz-1)                      &
      & + 1.0d0

    A2(1) =                                        &
      & (beta ** 2.0d0)                               &
      &   * xyz_F1BZ2(nx,ny,1)       &
      &   * xyr_F2BZ2(nx,ny,1)       &
      &     / r_dz(1)                             &
      &   * (DTS ** 2.0d0)                              &
      &   / z_dz(1)                               &
      & + 1.0d0                                         

    A2(nz) =                                        &
      & (beta ** 2.0d0)                               &
      &   * xyz_F1BZ2(nx,ny,nz)       &
      &   * xyr_F2BZ2(nx,ny,nz-1)     &
      &     / r_dz(nz-1)                           &
      &   * (DTS ** 2.0d0)                              &
      &   / z_dz(nz)                               &
      & + 1.0d0                                         
    
    B2(1+1:nz) =                                 &
      & - (beta ** 2.0d0)                                  &
      &   * xyz_F1BZ2(nx,ny,1:nz-1)    &
      &   * xyr_F2BZ2(nx,ny,1:nz-1)    &
      &   * (DTS ** 2.0d0)                                 &
      &   / ( r_dz(1:nz-1) * z_dz(1:nz-1) )
    
    C2(1: nz-1) =                           &
      & - ( beta ** 2.0d0 )                           &
      &   * xyz_F1BZ2(nx,ny,1+1:nz)    &
      &   * xyr_F2BZ2(nx,ny,1  :nz-1)  &
      &   * (DTS ** 2.0d0) &
      &   / ( r_dz(1:nz-1) * z_dz(1+1:nz) )

    write(*,*) "Dynamics_VI_init, A:   ", &
      & minval( A(1:nz)  - A2(1:nz) ), &
      & maxval( A(1:nz)  - A2(1:nz) ), &
      & minval( (A(1:nz) - A2(1:nz)) / A2(1:nz)), &
      & maxval( (A(1:nz) - A2(1:nz)) / A2(1:nz))

    write(*,*) "Dynamics_VI_init, B:   ", &
      & minval( B(2:nz)  - B2(2:nz) ), &
      & maxval( B(2:nz)  - B2(2:nz) ), &
      & minval( (B(2:nz) - B2(2:nz)) / B2(2:nz)), &
      & maxval( (B(2:nz) - B2(2:nz)) / B2(2:nz))

    write(*,*) "Dynamics_VI_init, C:   ", &
      & minval( C(1:nz-1)  - C2(1:nz-1) ), &
      & maxval( C(1:nz-1)  - C2(1:nz-1) ), &
      & minval( (C(1:nz-1) - C2(1:nz-1)) / C2(1:nz-1)), &
      & maxval( (C(1:nz-1) - C2(1:nz-1)) / C2(1:nz-1))

!    write(*,*) A(1), A2(1), C(1), C2(1)
!    do k = 1, nz-1
!      write(*,*) A(k), A2(k), B(k), B2(k), C(k), C2(k)
!    end do
!    write(*,*) A(nz), A2(nz), B(nz), B2(nz)



    !----------------------------------------------------------------
    ! 係数行列を LU 分解
    !
    !配列の大きさを保管
    N    = nz                   !係数行列/改行列の次数, 整合寸法
    M    = nx * ny              !方程式の組数
    NUD  = 1                    !係数行列の上三角部分の帯幅
    NLD  = 1                    !係数行列の下三角部分の帯幅
    NAL  = NLD                  !LU 分解の結果 L の整合寸法
    NA   = NUD + NLD + 1
    INFO = 0

    ! 配列の割り当て
    !
    allocate( AL1(N), IP(N) )

    ! 解行列の計算. LAPACK を使用. 
    !
    call DGTTRF(N, C, A, B, AL1, IP, INFO)
    
    ! 解のコンディションをチェック. 
    !
    if (INFO /= 0) then
      call MessageNotify("Error", "lapack_linear", "INFO is not 0")
      stop
    end if
    
  end subroutine Dynamics_VI_init
  

!!!--------------------------------------------------------------------!!!
  subroutine Exner_HEVI(                   &
    &  pyz_VelXAs, xqz_VelYAs, xyr_VelZNs, &
    &  xyz_VelDivNs, xyz_ExnerNs,          &
    &  xyr_DVelZDtNl, xyz_DExnerDtNl,      &
    &  xyz_ExnerAs                         &
    & )
    !
    !陰解法を用いたエクスナー関数の計算. 
    !

    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数
    real(DP), intent(in)   :: pyz_VelXAs(imin:imax,jmin:jmax,kmin:kmax) 
                                                           !速度 u [τ+Δτ]
    real(DP), intent(in)   :: xqz_VelYAs(imin:imax,jmin:jmax,kmin:kmax) 
                                                           !速度 v [τ+Δτ]
    real(DP), intent(in)   :: xyr_VelZNs(imin:imax,jmin:jmax,kmin:kmax) 
                                                           !速度 w [τ]
    real(DP), intent(in)   :: xyz_VelDivNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)   :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax) 
                                                           !Z 方向の外力項[t]
    real(DP), intent(in)   :: xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax) 
                                                           !Z 方向の外力項[t]
    real(DP), intent(in)   :: xyz_ExnerNs(imin:imax,jmin:jmax,kmin:kmax) 
                                                           !無次元圧力
    real(DP), intent(out)  :: xyz_ExnerAs(imin:imax,jmin:jmax,kmin:kmax) 
                                                           !無次元圧力[τ+Δτ]
    !作業変数定義
    real(DP)               :: D1(1:nx,1:ny,1:nz)  
    real(DP)               :: D(nx*ny,nz)
    real(DP)               :: E(1:nx,1:ny,0:nz)
!    real(DP)               :: E(imin:imax,jmin:jmax,kmin:kmax)    
    real(DP)               :: F(1:nx,1:ny,1:nz)
    real(DP)               :: F0(1:nx,1:ny,kmin:kmax-1)  
!    real(DP)               :: xyz_DivVelNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: dt ! 短い時間格子間隔
    integer                :: i, j, k

    real(DP)               :: X(M, N)     !定数/解行列
    real(DP)               :: TX(N, M)    !解行列を転置したもの
    integer                :: NRHS        
    integer                :: INFO
    integer                :: LDB
    character(1),parameter :: TRANS = 'N'

    real(DP)  :: dts
    real(DP)  :: nD1(imin:imax,jmin:jmax,kmin:kmax)  
    real(DP)  :: nD2(1:nx,1:ny,1:nz)  
!    real(DP)  :: nD(nx*ny,nz)
    real(DP)  :: nE(imin:imax,jmin:jmax,kmin:kmax)  
    real(DP)  :: nF(imin:imax,jmin:jmax,kmin:kmax)  
    character(15) :: info2

    real(DP)  :: tmpA(imin:imax,jmin:jmax,kmin:kmax)  
    real(DP)  :: tmpB(imin:imax,jmin:jmax,kmin:kmax)  
    real(DP)  :: tmpC(imin:imax,jmin:jmax,kmin:kmax)  

    ! Initialize
    !
    dt = DelTimeShort
    dts = DelTimeShort

    !---------------------------------------------------------------
    !行列計算のための係数

    !  添字の範囲は, 1:nx, 1:ny, 0:nz
    !  D(:,:,1) を求める時に D(:,:,0) の値が必要になるため. 
    !
    tmpA = xyr_dz_xyz( xyz_ExnerNs )   
    tmpB = xyr_dz_xyz( xyz_VelDivNs )
    tmpC = xyr_avr_xyz( CpDry * xyz_VPTempBZ )
    
    do k = kmin, kmax-1
      do j = jmin, jmax-1
        do i = imin, imax-1
          tmpC(i,j,k) =                  &
            & (                          &
            &  + xyz_VelDivNs(i,j,k+1)   & 
            &  - xyz_VelDivNs(i,j,k)     & 
            &  ) / dz         
        end do
      end do
    end do
    tmpC(imax,:,:) = 1.0d10 
    tmpC(:,jmax,:) = 1.0d10 
    tmpC(:,:,kmax) = 1.0d10 

!    i = 25
!    j = 23
!    write(*,*) "***", kmin
!    do k = kmin, kmax-1
!      write(*,*) k, xyz_VelDivNs(i,j,k+1), xyz_VelDivNs(i,j,k), tmpC(i,j,k), tmpB(i,j,k), tmpC(i,j,k) - tmpB(i,j,k) 
!    end do

    do k = 0, nz
      do j = 1, ny
        do i = 1, nx
          
          E(i,j,k) =                            &
            & - ( 1.0d0 - beta )                &
            &   * (                             &
            &     + xyz_ExnerNs(i,j,k+1)        & !
            &     - xyz_ExnerNs(i,j,k)          & ! xyz => xyr
            &     ) / dz                        &
            & + (                               &
            &   + Alpha                         &
            &     * (                           &
            &       + xyz_VelDivNs(i,j,k+1)     & !
            &       - xyz_VelDivNs(i,j,k)       & ! xyz => xyr
            &       ) / dz                      &
            &   + xyr_DVelZDtNl(i,j,k)          &
            &   ) / xyr_CpVPTempBZ(i,j,k)  

        end do
      end do
    end do

    nE =   &
      & - ( 1.0d0 - beta ) * xyr_dz_xyz( xyz_ExnerNs )  &
      & + ( Alpha * xyr_dz_xyz( xyz_VelDivNs ) + xyr_DVelZDtNl )    &
      &   / xyr_avr_xyz( CpDry * xyz_VPTempBZ ) 


    info2 = "VI, E: "
    call CHK2_Val( E(1:nx,1:ny,1:nz), nE(1:nx,1:ny,1:nz), info2 )


    ! 被微分関数
    !   配列 F0 の添字の範囲は, 1:nx, 1:ny, kmin:kmax
    !   配列 F を求める際に F0 を z 方向に微分するため. 
    !
    do k = kmin, kmax-1
      do j = 1, ny
        do i = 1, nx

          F0(i,j,k)  =                            &
            & + xyr_DensVPTempBZ(i,j,k)           &
            &   * (                               &
            &     + xyr_VelZNs(i,j,k)             & 
            &     - xyr_CpVPTempBZ(i,j,k)         &
            &       * (1.0d0 - beta)              &
            &       * (                           &
            &           xyz_ExnerNs(i,j,k+1)      & !
            &         - xyz_ExnerNs(i,j,k)        & ! xyz => xyr
            &         ) / dz * dt                 &
            &     + Alpha                         &
            &       * (                           &
            &           xyz_VelDivNs(i,j,k+1)     &
            &         - xyz_VelDivNs(i,j,k)       &
            &         ) / dz * dt                 &
            &     + xyr_DVelZDtNl(i,j,k) * dt     &
            &     )
        end do
      end do
    end do

    !行列計算のための係数
    !  添字の範囲は, 1:nx, 1:ny, 1:nz
    !
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx

          F(i,j,k) = &
            & - beta * xyz_F1BZ(i,j,k) * dt   &
            &   * (                           &
            &       F0(i,j,k)                 &
            &     - F0(i,j,k-1)               &
            &     ) / dz                      &
            & + xyz_DExnerDtNl(i,j,k) * dt 

        end do
      end do
    end do

    !行列計算のための係数
    !  添字の範囲は, 1:nx, 1:ny, 1:nz
    !
    do k = 1, nz 
      do j = 1, ny
        do i = 1, nx

          D1(i,j,k) =                                 &
            & + xyz_ExnerNs(i,j,k)                    &
            & - (1.0d0 - beta)                        &
            &   * xyz_F1BZ(i,j,k) * dt                &
            &   * (                                   &
            &       xyr_DensVPTempBZ(i,j,k)           &
            &       * xyr_VelZNs(i,j,k)               &
            &     - xyr_DensVPTempBZ(i,j,k-1)         &
            &       * xyr_VelZNs(i,j,k-1)             &
            &     ) / dz                              &
            & - (xyz_VelSW(i,j,k) ** 2.0d0) * dt      &
            &   / (CpDry * xyz_VPTempBZ(i,j,k))       &
            &   * (                                   &
            &     + (                                 &
            &         pyz_VelXAs(i,j,k)               &
            &       - pyz_VelXAs(i-1,j,k)             &
            &       ) / dx                            &
            &     + (                                 &
            &         xqz_VelYAs(i,j,k)               &
            &       - xqz_VelYAs(i,j-1,k)             &
            &       ) / dy                            &
            &     )                                   &
            & + F(i,j,k)
        end do
      end do
    end do

    ! 行列計算のための係数
    !
    do j = 1, ny
      do i = 1, nx

        D1(i,j,1) =                                    &
          & + D1(i,j,1)                                &
          & - beta * xyz_F1BZ(i,j,1) * (dt ** 2.0d0)   &
          &   * xyr_CpDensVPTemp2BZ(i,j,0)             &
          &   * E(i,j,0)                               &
          &   / dz
    
        D1(i,j,nz) =                                   &
          & + D1(i,j,nz)                               &
          & + beta * xyz_F1BZ(i,j,nz) * (dt ** 2.0d0)  &
          &   * xyr_CpDensVPTemp2BZ(i,j,nz)            &
          &   * E(i,j,nz)                              &
          &   / dz
      end do
    end do


    nF = - beta * xyz_F3BZ * DTS                                    &
      & * xyz_dz_xyr(                                               &
      &    xyr_avr_xyz( xyz_DensBZ * xyz_VPTempBZ)                  &
      &    * (                                                      &
      &         xyr_VelZNs                                          &
      &       - xyr_avr_xyz(CpDry * xyz_VPTempBZ) * DTS             &
      &         * (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs )        &
      &       + Alpha * xyr_dz_xyz( xyz_VelDivNs ) * DTS           &
      &       + xyr_DVelZDtNl * DTS                                 &
      &      )                                                      &
      &   )                                                         &
      & + ( xyz_DExnerDtNl ) * DTS

    nD1 = xyz_ExnerNs                                                &
      & - (1.0d0 - beta)                                            &
      &   * xyz_F3BZ * DTS                                          &
      &   * xyz_dz_xyr(                                             &
      &       xyr_avr_xyz(xyz_DensBZ * xyz_VPTempBZ) * xyr_VelZNs   &
      &     )                                                       &
      & - (xyz_VelSoundBZ ** 2.0d0) * DTS                           &
      &   / (CpDry * xyz_VPTempBZ) * xyz_dx_pyz( pyz_VelXAs )       &
      & - (xyz_VelSoundBZ ** 2.0d0) * DTS                           & 
      &   / (CpDry * xyz_VPTempBZ) * xyz_dy_xqz( xqz_VelYAs )       &
      & + nF

    nD1(1:nx,1:ny,1) = nD1(1:nx,1:ny,1)                           &
      & - beta * xyz_F3BZ(1:nx,1:ny,1) * (DTS ** 2.0d0)         &
      &   * xyr_F2BZ(1:nx,1:ny, 1-1) * nE(1:nx,1:ny,1-1)         &
      &   / z_dz(1)
    
    nD1(1:nx,1:ny,nz) = nD1(1:nx,1:ny,nz)                           &
      & + beta * xyz_F3BZ(1:nx,1:ny,nz) * (DTS ** 2.0d0)         &
      &   * xyr_F2BZ(1:nx,1:ny,nz) * nE(1:nx,1:ny,nz)              &
      &   / z_dz(nz)     

    nD2 = nD1(1:nx,1:ny,1:nz)


    info2 = "VI, F: "
    call CHK2_Val( F(1:nx,1:ny,1:nz), nF(1:nx,1:ny,1:nz), info2 )

    info2 = "VI, D1: "
    call CHK2_Val( D1(1:nx,1:ny,1:nz), nD2(1:nx,1:ny,1:nz), info2 )


    !-----------------------------------------------------------
    !連立一次方程式の解を求める

    !変数の初期化
    !
    NRHS = M
    INFO = 0
    LDB  = N

    ! LAPACK の仕様に合わせて変形 
    !
    do k = 1, nz
      do j = 1, ny 	
        do i = 1, nx 	
	  D(i + nx * (j - 1), k) =  D1(i,j,k)
        end do 
      end do
    end do	

    TX = transpose( D )
    
    !解行列の計算. LAPACK を使用. 
    !
    call DGTTRS(TRANS, N, NRHS, C, A, B, AL1, IP, TX, LDB, INFO)

    !解のコンディションをチェック. 
    !
    if (INFO /= 0) then
      call MessageNotify("Error", "lapack_linear", "INFO is not 0")
      stop
    end if

    !戻り値を出力
    !
    X = transpose( TX )
     
    do k = 1, nz
      do j = 1, ny 	
        do i = 1, nx 	
          xyz_ExnerAs(i,j,k) = X(i + nx * (j - 1 ), k)
        end do 
      end do
    end do	
    xyz_ExnerAs(imin:0,:,:) = 0.0d0
    xyz_ExnerAs(:,jmin:0,:) = 0.0d0
    xyz_ExnerAs(:,:,kmin:0) = 0.0d0
    xyz_ExnerAs(nx+1:imax,:,:) = 0.0d0
    xyz_ExnerAs(:,ny+1:jmax,:) = 0.0d0
    xyz_ExnerAs(:,:,nz+1:kmax) = 0.0d0

    ! のり代に値を入れる
    !
    call SetMargin_xyz( xyz_ExnerAs ) ! (inout)
  
  end subroutine Exner_HEVI

  
  
  subroutine AdvC4_nDiff_xyz(              &
    &  xyz_VarBl,  xyz_VarNl,              & !(IN)
    &  pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN)
    &  xyz_Adv, xyz_nDiff )                  !(OUT)

    implicit none

    real(DP), intent(in)  :: xyz_VarBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyz_VarNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)              :: fct1, fct2
    integer               :: i, j, k

    ! 微分に用いる係数を予め計算
    !
    fct1 = 9.0d0 / 8.0d0
    fct2 = 1.0d0 / 24.0d0

    ! 移流項の計算. 移流項: 4 次精度中心差分
    ! 
    do k = kmin + 2, kmax - 2
      do j = jmin + 2, jmax - 2
        do i = imin + 2, imax - 2

          xyz_Adv(i,j,k) =                                                  &
            & - (                                                           &
            &      pyz_VelXNl(i,j,k)                                        &
            &        * (                                                    &
            &            fct1 * ( xyz_VarNl(i+1,j,k) - xyz_VarNl(i,j,k)   ) &
            &          - fct2 * ( xyz_VarNl(i+2,j,k) - xyz_VarNl(i-1,j,k) ) &
            &          )                                                    &
            &    + pyz_VelXNl(i-1,j,k)                                      &
            &        * (                                                    &
            &            fct1 * ( xyz_VarNl(i,j,k)   - xyz_VarNl(i-1,j,k) ) &
            &          - fct2 * ( xyz_VarNl(i+1,j,k) - xyz_VarNl(i-2,j,k) ) &
            &          )                                                    &
            &   ) * 5.0d-1 / dx                                             &
            & - (                                                           &
            &      xqz_VelYNl(i,j,k)                                        &
            &        * (                                                    &
            &            fct1 * ( xyz_VarNl(i,j+1,k) - xyz_VarNl(i,j,k)   ) &
            &          - fct2 * ( xyz_VarNl(i,j+2,k) - xyz_VarNl(i,j-1,k) ) &
            &          )                                                    &
            &    + xqz_VelYNl(i,j-1,k)                                      &
            &        * (                                                    &
            &            fct1 * ( xyz_VarNl(i,j,k)   - xyz_VarNl(i,j-1,k) ) &
            &          - fct2 * ( xyz_VarNl(i,j+1,k) - xyz_VarNl(i,j-2,k) ) &
            &          )                                                    &
            &   ) * 5.0d-1 / dy                                             &
            & - (                                                           &
            &      xyr_VelZNl(i,j,k)                                        &
            &        * (                                                    &
            &            fct1 * ( xyz_VarNl(i,j,k+1) - xyz_VarNl(i,j,k)   ) &
            &          - fct2 * ( xyz_VarNl(i,j,k+2) - xyz_VarNl(i,j,k-1) ) &
            &          )                                                    &
            &    + xyr_VelZNl(i,j,k-1)                                      &
            &        * (                                                    &
            &            fct1 * ( xyz_VarNl(i,j,k)   - xyz_VarNl(i,j,k-1) ) &
            &          - fct2 * ( xyz_VarNl(i,j,k+1) - xyz_VarNl(i,j,k-2) ) &
            &          )                                                    &
            &   ) * 5.0d-1 / dz
        end do
      end do
    end do
    
    xyz_Adv(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0
    xyz_Adv(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0
    xyz_Adv(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0
    xyz_Adv(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0
    xyz_Adv(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0
    xyz_Adv(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0

    ! 4 次の数値拡散: 2 次精度中心差分
    ! 
    do k = kmin + 2, kmax - 2
      do j = jmin + 2, jmax - 2
        do i = imin + 2, imax - 2
          
          xyz_nDiff(i,j,k) =                            &
            & - (                                       &
            &     + xyz_VarBl(i+2,j,k)                  &
            &     + xyz_VarBl(i-2,j,k)                  &
            &     - xyz_VarBl(i+1,j,k) * 4.0d0          &
            &     - xyz_VarBl(i-1,j,k) * 4.0d0          &
            &     + xyz_VarBl(i  ,j,k) * 6.0d0          &
            &   ) * NuHh / ( dx ** 4.0d0 )              &
            & - (                                       &
            &       xyz_VarBl(i,j+2,k)                  &
            &     + xyz_VarBl(i,j-2,k)                  &
            &     - xyz_VarBl(i,j+1,k) * 4.0d0          &
            &     - xyz_VarBl(i,j-1,k) * 4.0d0          &
            &     + xyz_VarBl(i,j  ,k) * 6.0d0          &
            &   ) * NuHh / ( dy ** 4.0d0 )              &
            & - (                                       &
            &       xyz_VarBl(i,j,k+2)                  &
            &     + xyz_VarBl(i,j,k-2)                  &
            &     - xyz_VarBl(i,j,k+1) * 4.0d0          &
            &     - xyz_VarBl(i,j,k-1) * 4.0d0          &
            &     + xyz_VarBl(i,j,k  ) * 6.0d0          &
            &   ) * NuVh / ( dz ** 4.0d0 )
        end do
      end do
    end do

    xyz_nDiff(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0
    xyz_nDiff(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0
    xyz_nDiff(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0
    xyz_nDiff(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0
    xyz_nDiff(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0
    xyz_nDiff(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0

  end subroutine AdvC4_nDiff_xyz


  subroutine AdvC4_nDiff_xyzf(               &
      &  xyzf_VarBl, xyzf_VarNl,             & !(IN)
      &  pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN)
      &  xyzf_Adv, xyzf_nDiff )                !(OUT)

    implicit none

    real(DP), intent(in)  :: xyzf_VarBl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP), intent(in)  :: xyzf_VarNl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP), intent(in)  :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyzf_Adv(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP), intent(out) :: xyzf_nDiff(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)

    real(DP)              :: fct1, fct2
    integer              :: s, i, j, k

    ! 微分に用いる係数を予め計算
    !
    fct1 = 9.0d0 / 8.0d0
    fct2 = 1.0d0 / 24.0d0

    ! 移流項の計算. 移流項: 4 次精度中心差分
    ! 
    do s = 1, ncmax
      do k = kmin + 2, kmax - 2
        do j = jmin + 2, jmax - 2
          do i = imin + 2, imax - 2
            
            xyzf_Adv(i,j,k,s) =                                                     &
              & - (                                                                 &
              &      pyz_VelXNl(i,j,k)                                              &
              &        * (                                                          &
              &            fct1 * ( xyzf_VarNl(i+1,j,k,s) - xyzf_VarNl(i,j,k,s)   ) &
              &          - fct2 * ( xyzf_VarNl(i+2,j,k,s) - xyzf_VarNl(i-1,j,k,s) ) &
              &          )                                                          &
              &    + pyz_VelXNl(i-1,j,k)                                            &
              &        * (                                                          &
              &            fct1 * ( xyzf_VarNl(i,j,k,s)   - xyzf_VarNl(i-1,j,k,s) ) &
              &          - fct2 * ( xyzf_VarNl(i+1,j,k,s) - xyzf_VarNl(i-2,j,k,s) ) &
              &          )                                                          &
              &   ) * 5.0d-1 / dx                                                   &
              & - (                                                                 &
              &      xqz_VelYNl(i,j,k)                                              &
              &        * (                                                          &
              &            fct1 * ( xyzf_VarNl(i,j+1,k,s) - xyzf_VarNl(i,j,k,s)   ) &
              &          - fct2 * ( xyzf_VarNl(i,j+2,k,s) - xyzf_VarNl(i,j-1,k,s) ) &
              &          )                                                          &
              &    + xqz_VelYNl(i,j-1,k)                                            &
              &        * (                                                          &
              &            fct1 * ( xyzf_VarNl(i,j,k,s)   - xyzf_VarNl(i,j-1,k,s) ) &
              &          - fct2 * ( xyzf_VarNl(i,j+1,k,s) - xyzf_VarNl(i,j-2,k,s) ) &
              &          )                                                          &
              &   ) * 5.0d-1 / dy                                                   &
              & - (                                                                 &
              &      xyr_VelZNl(i,j,k)                                              &
              &        * (                                                          &
              &            fct1 * ( xyzf_VarNl(i,j,k+1,s) - xyzf_VarNl(i,j,k,s)   ) &
              &          - fct2 * ( xyzf_VarNl(i,j,k+2,s) - xyzf_VarNl(i,j,k-1,s) ) &
              &          )                                                          &
              &    + xyr_VelZNl(i,j,k-1)                                            &
              &        * (                                                          &
              &            fct1 * ( xyzf_VarNl(i,j,k,s)   - xyzf_VarNl(i,j,k-1,s) ) &
              &          - fct2 * ( xyzf_VarNl(i,j,k+1,s) - xyzf_VarNl(i,j,k-2,s) ) &
              &          )                                                          &
              &   ) * 5.0d-1 / dz
          end do
        end do
      end do
    end do

    xyzf_Adv(imin:imin+1,jmin:jmax,kmin:kmax,1:ncmax) = 0.0d0
    xyzf_Adv(imax-1:imax,jmin:jmax,kmin:kmax,1:ncmax) = 0.0d0
    xyzf_Adv(imin:imax,jmin:jmin+1,kmin:kmax,1:ncmax) = 0.0d0
    xyzf_Adv(imin:imax,jmax-1:jmax,kmin:kmax,1:ncmax) = 0.0d0
    xyzf_Adv(imin:imax,jmin:jmax,kmin:kmin+1,1:ncmax) = 0.0d0
    xyzf_Adv(imin:imax,jmin:jmax,kmax-1:kmax,1:ncmax) = 0.0d0

    ! 数値拡散: 2 次精度中心差分
    ! 
    do s = 1, ncmax
      do k = kmin + 2, kmax - 2
        do j = jmin + 2, jmax - 2
          do i = imin + 2, imax - 2

            xyzf_nDiff(i,j,k,s) =                            &
              & - (                                          &
              &       xyzf_VarBl(i+2,j,k,s)                  &
              &     - xyzf_VarBl(i+1,j,k,s) * 4.0d0          &
              &     + xyzf_VarBl(i  ,j,k,s) * 6.0d0          &
              &     - xyzf_VarBl(i-1,j,k,s) * 4.0d0          &
              &     + xyzf_VarBl(i-2,j,k,s)                  &
              &   ) * NuHh / ( dx ** 4.0d0 )                 &
              & - (                                          &
              &       xyzf_VarBl(i,j+2,k,s)                  &
              &     - xyzf_VarBl(i,j+1,k,s) * 4.0d0          &
              &     + xyzf_VarBl(i,j  ,k,s) * 6.0d0          &
              &     - xyzf_VarBl(i,j-1,k,s) * 4.0d0          &
              &     + xyzf_VarBl(i,j-2,k,s)                  &
              &   ) * NuHh / ( dy ** 4.0d0 )                 &
              & - (                                          &
              &       xyzf_VarBl(i,j,k+2,s)                  &
              &     - xyzf_VarBl(i,j,k+1,s) * 4.0d0          &
              &     + xyzf_VarBl(i,j,k  ,s) * 6.0d0          &
              &     - xyzf_VarBl(i,j,k-1,s) * 4.0d0          &
              &     + xyzf_VarBl(i,j,k-2,s)                  &
              &   ) * NuVh / ( dz ** 4.0d0 )
          end do
        end do
      end do
    end do

    xyzf_nDiff(imin:imin+1,jmin:jmax,kmin:kmax,1:ncmax) = 0.0d0
    xyzf_nDiff(imax-1:imax,jmin:jmax,kmin:kmax,1:ncmax) = 0.0d0
    xyzf_nDiff(imin:imax,jmin:jmin+1,kmin:kmax,1:ncmax) = 0.0d0
    xyzf_nDiff(imin:imax,jmax-1:jmax,kmin:kmax,1:ncmax) = 0.0d0
    xyzf_nDiff(imin:imax,jmin:jmax,kmin:kmin+1,1:ncmax) = 0.0d0
    xyzf_nDiff(imin:imax,jmin:jmax,kmax-1:kmax,1:ncmax) = 0.0d0

  end subroutine AdvC4_nDiff_xyzf


  subroutine AdvC4_nDiff_pyz_xqz_xyr(        & 
      &  pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, & !(IN)
      &  pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN)
      &  pyz_Adv,    xqz_Adv,    xyr_Adv,    & !(OUT)
      &  pyz_nDiff,  xqz_nDiff,  xyr_nDiff )   !(OUT)

    implicit none
    
    real(DP), intent(in)  :: pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: pyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xqz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyr_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: pyz_nDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xqz_nDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyr_nDiff(imin:imax,jmin:jmax,kmin:kmax)    

    real(DP)              :: fct1, fct2
    integer               :: i, j, k

    ! 微分に用いる係数を予め計算
    !
    fct1 = 9.0d0 / 8.0d0
    fct2 = 1.0d0 / 24.0d0

    ! 移流項の計算. 移流項: 4 次精度中心差分
    ! 
    do k = kmin + 2, kmax - 2
      do j = jmin + 2, jmax - 2
        do i = imin + 2, imax - 2
          
          pyz_Adv(i,j,k) =                                                   &
            & - pyz_VelXNl(i,j,k)                                            &
            &   * (                                                          &
            &       fct1 * (   pyz_VelXNl(i+1,j,k) - pyz_VelXNl(i-1,j,k) )   &
            &     - fct2 * (   pyz_VelXNl(i+2,j,k) + pyz_VelXNl(i+1,j,k)     &
            &                - pyz_VelXNl(i-1,j,k) - pyz_VelXNl(i-2,j,k) )   &
            &     ) * 5.0d-1 / dx                                            &
            & - (                                                            &
            &   + ( xqz_VelYNl(i+1,j,k) + xqz_VelYNl(i,j,k) )                &  
            &     * (                                                        &  
            &         fct1 * ( pyz_VelXNl(i,j+1,k) - pyz_VelXNl(i,j,k)   )   &
            &       - fct2 * ( pyz_VelXNl(i,j+2,k) - pyz_VelXNl(i,j-1,k) )   &
            &       )                                                        &
            &   + ( xqz_VelYNl(i+1,j-1,k) + xqz_VelYNl(i,j-1,k) )            & 
            &     * (                                                        &
            &         fct1 * ( pyz_VelXNl(i,j,k)   - pyz_VelXNl(i,j-1,k) )   &
            &       - fct2 * ( pyz_VelXNl(i,j+1,k) - pyz_VelXNl(i,j-2,k) )   &
            &       )                                                        &
            &   ) * 2.5d-1 / dy                                              &
            & - (                                                            &
            &   + ( xyr_VelZNl(i+1,j,k) + xyr_VelZNl(i,j,k) )                & 
            &     * (                                                        &
            &         fct1 * ( pyz_VelXNl(i,j,k+1) - pyz_VelXNl(i,j,k)   )   &
            &       - fct2 * ( pyz_VelXNl(i,j,k+2) - pyz_VelXNl(i,j,k-1) )   &
            &       )                                                        &
            &   + ( xyr_VelZNl(i+1,j,k-1) + xyr_VelZNl(i,j,k-1) )            & 
            &     * (                                                        &
            &         fct1 * ( pyz_VelXNl(i,j,k)   - pyz_VelXNl(i,j,k-1) )   &
            &       - fct2 * ( pyz_VelXNl(i,j,k+1) - pyz_VelXNl(i,j,k-2) )   &
            &       )                                                        &
            &   ) * 2.5d-1 / dz

        end do
      end do
    end do

    pyz_Adv(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0
    pyz_Adv(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0
    pyz_Adv(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0
    pyz_Adv(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0
    pyz_Adv(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0
    pyz_Adv(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0

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

          xqz_Adv(i,j,k) =                                                   &
            & - (                                                            &
            &   + ( pyz_VelXNl(i,j+1,k) + pyz_VelXNl(i,j,k) )                & 
            &     * (                                                        &
            &         fct1 * ( xqz_VelYNl(i+1,j,k) - xqz_VelYNl(i,j,k)   )   &
            &       - fct2 * ( xqz_VelYNl(i+2,j,k) - xqz_VelYNl(i-1,j,k) )   &
            &       )                                                        &
            &   + ( pyz_VelXNl(i-1,j+1,k) + pyz_VelXNl(i-1,j,k) )            & 
            &     * (                                                        &
            &         fct1 * ( xqz_VelYNl(i,j,k)   - xqz_VelYNl(i-1,j,k) )   &
            &       - fct2 * ( xqz_VelYNl(i+1,j,k) - xqz_VelYNl(i-2,j,k) )   &
            &       )                                                        &
            &   ) * 2.5d-1 / dx                                              &
            & - xqz_VelYNl(i,j,k)                                            &
            &   * (                                                          &
            &       fct1 * (   xqz_VelYNl(i,j+1,k) - xqz_VelYNl(i,j-1,k) )   &
            &     - fct2 * (   xqz_VelYNl(i,j+2,k) + xqz_VelYNl(i,j+1,k)     &
            &                - xqz_VelYNl(i,j-1,k) - xqz_VelYNl(i,j-2,k) )   &
            &     ) * 5.0d-1 / dy                                            &
            & - (                                                            &
            &   + ( xyr_VelZNl(i,j+1,k) + xyr_VelZNl(i,j,k) )                & 
            &     * (                                                        &
            &         fct1 * ( xqz_VelYNl(i,j,k+1) - xqz_VelYNl(i,j,k)   )   &
            &       - fct2 * ( xqz_VelYNl(i,j,k+2) - xqz_VelYNl(i,j,k-1) )   &
            &       )                                                        &
            &   + ( xyr_VelZNl(i,j+1,k-1) + xyr_VelZNl(i,j,k-1) )            & 
            &     * (                                                        &
            &         fct1 * ( xqz_VelYNl(i,j,k)   - xqz_VelYNl(i,j,k-1) )   &
            &       - fct2 * ( xqz_VelYNl(i,j,k+1) - xqz_VelYNl(i,j,k-2) )   &
            &       )                                                        &
            &   ) * 2.5d-1 / dz
        end do
      end do
    end do

    xqz_Adv(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0
    xqz_Adv(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0
    xqz_Adv(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0
    xqz_Adv(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0
    xqz_Adv(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0
    xqz_Adv(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0

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

          xyr_Adv(i,j,k) =                                                   &
            & - (                                                            &
            &   + ( pyz_VelXNl(i,j,k+1) + pyz_VelXNl(i,j,k) )                & 
            &     * (                                                        &
            &         fct1 * ( xyr_VelZNl(i+1,j,k) - xyr_VelZNl(i,j,k)   )   &
            &       - fct2 * ( xyr_VelZNl(i+2,j,k) - xyr_VelZNl(i-1,j,k) )   &
            &       )                                                        &
            &   + ( pyz_VelXNl(i-1,j,k+1) + pyz_VelXNl(i-1,j,k) )            & 
            &     * (                                                        &
            &         fct1 * ( xyr_VelZNl(i,j,k)   - xyr_VelZNl(i-1,j,k) )   &
            &       - fct2 * ( xyr_VelZNl(i+1,j,k) - xyr_VelZNl(i-2,j,k) )   &
            &       )                                                        &
            &   ) * 2.5d-1 / dx                                              &
            & - (                                                            &
            &   + ( xqz_VelYNl(i,j,k+1) + xqz_VelYNl(i,j,k) )                & 
            &     * (                                                        &
            &         fct1 * ( xyr_VelZNl(i,j+1,k) - xyr_VelZNl(i,j,k)   )   &
            &       - fct2 * ( xyr_VelZNl(i,j+2,k) - xyr_VelZNl(i,j-1,k) )   &
            &       )                                                        &
            &   + ( xqz_VelYNl(i,j-1,k+1) + xqz_VelYNl(i,j-1,k) )            & 
            &     * (                                                        &
            &         fct1 * ( xyr_VelZNl(i,j,k)   - xyr_VelZNl(i,j-1,k) )   &
            &       - fct2 * ( xyr_VelZNl(i,j+1,k) - xyr_VelZNl(i,j-2,k) )   &
            &       )                                                        &
            &   ) * 2.5d-1 / dy                                              &
            & - xyr_VelZNl(i,j,k)                                            &
            &   * (                                                          &
            &       fct1 * (   xyr_VelZNl(i,j,k+1) - xyr_VelZNl(i,j,k-1) )   &
            &     - fct2 * (   xyr_VelZNl(i,j,k+2) + xyr_VelZNl(i,j,k+1)     &
            &                - xyr_VelZNl(i,j,k-1) - xyr_VelZNl(i,j,k-2) )   &
            &     ) * 5.0d-1 / dz
        end do
      end do
    end do

    xyr_Adv(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0
    xyr_Adv(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0
    xyr_Adv(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0
    xyr_Adv(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0
    xyr_Adv(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0
    xyr_Adv(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0

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

          pyz_nDiff(i,j,k) =                             &
            & - (                                        &
            &     + pyz_VelXBl(i+2,j,k)                  &
            &     + pyz_VelXBl(i-2,j,k)                  &
            &     - pyz_VelXBl(i+1,j,k) * 4.0d0          &
            &     - pyz_VelXBl(i-1,j,k) * 4.0d0          &
            &     + pyz_VelXBl(i  ,j,k) * 6.0d0          &
            &   ) * NuHm / ( dx ** 4.0d0 )               &
            & - (                                        &
            &     + pyz_VelXBl(i,j+2,k)                  &
            &     + pyz_VelXBl(i,j-2,k)                  &
            &     - pyz_VelXBl(i,j+1,k) * 4.0d0          &
            &     - pyz_VelXBl(i,j-1,k) * 4.0d0          &
            &     + pyz_VelXBl(i,j  ,k) * 6.0d0          &
            &   ) * NuHm / ( dy ** 4.0d0 )               &
            & - (                                        & 
            &     + pyz_VelXBl(i,j,k+2)                  &
            &     + pyz_VelXBl(i,j,k-2)                  &
            &     - pyz_VelXBl(i,j,k+1) * 4.0d0          &
            &     - pyz_VelXBl(i,j,k-1) * 4.0d0          &
            &     + pyz_VelXBl(i,j,k  ) * 6.0d0          &
            &   ) * NuVm / ( dz ** 4.0d0 )

        end do
      end do
    end do

    pyz_nDiff(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0
    pyz_nDiff(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0
    pyz_nDiff(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0
    pyz_nDiff(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0
    pyz_nDiff(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0
    pyz_nDiff(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0

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

          xqz_nDiff(i,j,k) =                             &
            & - (                                        &
            &     + xqz_VelYBl(i+2,j,k)                  &
            &     + xqz_VelYBl(i-2,j,k)                  &
            &     - xqz_VelYBl(i+1,j,k) * 4.0d0          &
            &     - xqz_VelYBl(i-1,j,k) * 4.0d0          &
            &     + xqz_VelYBl(i  ,j,k) * 6.0d0          &
            &   ) * NuHm / ( dx ** 4.0d0 )               &
            & - (                                        &
            &     + xqz_VelYBl(i,j+2,k)                  &
            &     + xqz_VelYBl(i,j-2,k)                  &
            &     - xqz_VelYBl(i,j+1,k) * 4.0d0          &
            &     - xqz_VelYBl(i,j-1,k) * 4.0d0          &
            &     + xqz_VelYBl(i,j  ,k) * 6.0d0          &
            &   ) * NuHm / ( dy ** 4.0d0 )               &
            & - (                                        & 
            &     + xqz_VelYBl(i,j,k+2)                  &
            &     + xqz_VelYBl(i,j,k-2)                  &
            &     - xqz_VelYBl(i,j,k+1) * 4.0d0          &
            &     - xqz_VelYBl(i,j,k-1) * 4.0d0          &
            &     + xqz_VelYBl(i,j,k  ) * 6.0d0          &
            &   ) * NuVm / ( dz ** 4.0d0 )
        end do
      end do
    end do
    
    xqz_nDiff(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0
    xqz_nDiff(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0
    xqz_nDiff(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0
    xqz_nDiff(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0
    xqz_nDiff(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0
    xqz_nDiff(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0

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

          xyr_nDiff(i,j,k) =                             &
            & - (                                        &
            &     + xyr_VelZBl(i+2,j,k)                  &
            &     + xyr_VelZBl(i-2,j,k)                  &
            &     - xyr_VelZBl(i+1,j,k) * 4.0d0          &
            &     - xyr_VelZBl(i-1,j,k) * 4.0d0          &
            &     + xyr_VelZBl(i  ,j,k) * 6.0d0          &
            &   ) * NuHm / ( dx ** 4.0d0 )               &
            & - (                                        &
            &     + xyr_VelZBl(i,j+2,k)                  &
            &     + xyr_VelZBl(i,j-2,k)                  &
            &     - xyr_VelZBl(i,j+1,k) * 4.0d0          &
            &     - xyr_VelZBl(i,j-1,k) * 4.0d0          &
            &     + xyr_VelZBl(i,j  ,k) * 6.0d0          &
            &   ) * NuHm / ( dy ** 4.0d0 )               &
            & - (                                        &
            &     + xyr_VelZBl(i,j,k+2)                  &
            &     + xyr_VelZBl(i,j,k-2)                  &
            &     - xyr_VelZBl(i,j,k+1) * 4.0d0          &
            &     - xyr_VelZBl(i,j,k-1) * 4.0d0          &
            &     + xyr_VelZBl(i,j,k  ) * 6.0d0          &
            &   ) * NuVm / ( dz ** 4.0d0 )
        end do
      end do
    end do

    xyr_nDiff(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0
    xyr_nDiff(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0
    xyr_nDiff(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0
    xyr_nDiff(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0
    xyr_nDiff(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0
    xyr_nDiff(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0

  end subroutine AdvC4_nDiff_pyz_xqz_xyr
  

  subroutine BuoyancyLong_xyr(           &
    &  xyz_PTempNl, xyzf_QMixNl,         & !(IN)
    &  xyr_BuoyT, xyr_BuoyM, xyr_BuoyD )   !(OUT)

    use composition, only: GasNum,       &! 
      &                    IdxG,         &!
      &                    MolWtWet       ! 湿潤成分の分子量
    use constants,only:    MolWtDry,     &! 乾燥成分の分子量
      &                    Grav           ! 重力加速度
    
    implicit none

    real(DP), intent(in)  :: xyz_PTempNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyzf_QMixNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(out) :: xyr_BuoyT(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyr_BuoyM(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyr_BuoyD(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)              :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, 1:GasNum)
    real(DP)              :: tmp1(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)              :: tmp2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)              :: tmp3(imin:imax,jmin:jmax,kmin:kmax)
    integer               :: i, j, k, f, n

    real(DP)             :: xyr_BuoyT2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_BuoyM2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_BuoyD2(imin:imax,jmin:jmax,kmin:kmax)

    character(15)        :: info2


    do f = 1, GasNum
      n = IdxG(f)
      xyzf_QMixPerMolWt(:,:,:,f) = xyzf_QMixNl(:,:,:,n) / MolWtWet(n)
    end do

    ! 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_PTempBZ(i,j,k+1) &
            &   + xyz_PTempNl(i,j,k)   / xyz_PTempBZ(i,j,k)   &
            &   ) * 5.0d-1

        end do
      end do
    end do

    xyr_BuoyT(:,:,kmax) = 0.0d0

    ! Buoyancy due to molecular weight
    !
    tmp1 = sum(xyzf_QMixPerMolWt, 4) 
    tmp2 = 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                                               &
            &   * ( tmp1(i,j,k+1) + tmp1(i,j,k) ) * 5.0d-1         &
            &   / ( 1.0d0 / MolWtDry + xyr_QMixBZPerMolWt(i,j,k) ) &
            & - Grav                                               &
            &   * ( tmp2(i,j,k+1) + tmp2(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
    !
    tmp3 = 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                                        &
            &   * ( tmp3(i,j,k+1) + tmp3(i,j,k) ) * 5.0d-1  &
            &   / ( 1.0d0 + xyr_QMixBZ(i,j,k) )

        end do
      end do
    end do

    xyr_BuoyD(:,:,kmax) = 0.0d0

    xyr_BuoyT2 = Grav * xyr_avr_xyz( xyz_PTempNl / xyz_PTempBZ)

    xyr_BuoyM2 =                                                     &
      & + Grav * xyr_avr_xyz( sum(xyzf_QMixPerMolWt, 4) )           &
      &    / ( 1.0d0 / MolWtDry + xyr_QMixBZPerMolWt )              &
      & - Grav * xyr_avr_xyz( sum(xyzf_QMixNl(:,:,:,1:GasNum), 4) ) &
      &    / ( 1.0d0 + xyr_QmixBZ ) 

    xyr_BuoyD2 =                                                           &
      & - Grav * xyr_avr_xyz( sum(xyzf_QMixNl(:,:,:,GasNum+1:ncmax), 4) ) &
      &    / ( 1.0d0 + xyr_QMixBZ )
    
    info2 = "VelZ, BuoyT: "
    call CHK_Val( xyr_BuoyT(imin:imax,jmin:jmax,kmin:kmax), xyr_BuoyT2(imin:imax,jmin:jmax,kmin:kmax), info2 )

    info2 = "VelZ, BuoyM: "
    call CHK_Val( xyr_BuoyM(imin:imax,jmin:jmax,kmin:kmax), xyr_BuoyM2(imin:imax,jmin:jmax,kmin:kmax), info2 )

    info2 = "VelZ, BuoyT: "
    call CHK_Val( xyr_BuoyD(imin:imax,jmin:jmax,kmin:kmax), xyr_BuoyD2(imin:imax,jmin:jmax,kmin:kmax), info2 )
    
  end subroutine BuoyancyLong_xyr


  subroutine VelDivC2(               &
    &  pyz_VelX, xqz_VelY, xyr_VelZ, & !(IN)
    &  xyz_VelDiv  )                   !(OUT)

    implicit none

    real(DP), intent(in)  :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZ(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyz_VelDiv(imin:imax,jmin:jmax,kmin:kmax)
    integer              :: i, j, k
    
    do  k = kmin + 1, kmax 
      do j = jmin + 1, jmax 
        do i = imin + 1, imax 
          xyz_VelDiv(i,j,k) =          &
            & + (                      &
            &     pyz_VelX(i,j,k)      &
            &   - pyz_VelX(i-1,j,k)    &
            &   ) / dx                 &
            & + (                      &
            &     xqz_VelY(i,j,k)      &
            &   - xqz_VelY(i,j-1,k)    &
            &   ) / dy                 &
            & + (                      &
            &     xyr_VelZ(i,j,k)      &
            &   - xyr_VelZ(i,j,k-1)    &
            &   ) / dz
        end do
      end do
    end do

    ! 不定を作らないために
    !
!    xyz_VelDiv(imin,:,:) = 0.0d0 
!    xyz_VelDiv(:,jmin,:) = 0.0d0 
!    xyz_VelDiv(:,:,kmin) = 0.0d0 
    xyz_VelDiv(imin,:,:) = 1.0d10 
    xyz_VelDiv(:,jmin,:) = 1.0d10 
    xyz_VelDiv(:,:,kmin) = 1.0d10 

  end subroutine VelDivC2
  
  
  
  subroutine PGrad_HE(         &
    &  xyz_VelDiv, xyz_Exner,  & !(IN)
    &  pyz_PGrad, pyz_SWF,     & !(OUT)
    &  xqz_PGrad, xqz_SWF )      !(OUT)

    
    implicit none
    
    real(DP), intent(in)  :: xyz_VelDiv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: pyz_PGrad(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xqz_PGrad(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: pyz_SWF(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xqz_SWF(imin:imax,jmin:jmax,kmin:kmax)
    integer              :: i, j, k

    !------------------------------------------------------------------
    ! X 方向

    ! 音波減衰項
    !
    do k = kmin, kmax
      do j = jmin, jmax
        do i = imin, imax-1
          
          pyz_SWF(i,j,k) =            &
            &   Alpha                     &
            &   * (                       &
            &       xyz_VelDiv(i+1,j,k)   &
            &     - xyz_VelDiv(i,j,k)     &
            &     ) / dx

        end do
      end do
    end do

    ! 圧力傾度力
    !
    do k = kmin, kmax
      do j = jmin, jmax
        do i = imin, imax-1          
          
          pyz_PGrad(i,j,k) =                &
            & - CpDry * pyz_VPTempBZ(i,j,k) &
            &   * (                         &
            &       xyz_Exner(i+1,j,k)      &
            &     - xyz_Exner(i,j,k)        &
            &     ) / dx                    
          
        end do
      end do
    end do
    
!    write(*,*) CpDry
!    write(*,*) pyz_VPTempBZ(1,1,:)

    ! 穴埋め
    !
    pyz_SWF(imax,:,:)   = 0.0d0 
    pyz_PGrad(imax,:,:) = 0.0d0 
    
    !------------------------------------------------------------------
    ! Y 方向

    ! 音波減衰項
    !
    do k = kmin, kmax
      do j = jmin, jmax-1
        do i = imin, imax
          
          xqz_SWF(i,j,k) =                &
            &   Alpha                     &
            &   * (                       &
            &       xyz_VelDiv(i,j+1,k)   &
            &     - xyz_VelDiv(i,j,k)     &
            &     ) / dy

        end do
      end do
    end do

    ! 圧力傾度力
    ! 
    do k = kmin, kmax
      do j = jmin, jmax-1
        do i = imin, imax
          
          xqz_PGrad(i,j,k) =                &
            & - CpDry * xqz_VPTempBZ(i,j,k) &
            &   * (                         &
            &       xyz_Exner(i,j+1,k)      &
            &     - xyz_Exner(i,j,k)        &
            &     ) /dy                     
          
        end do
      end do
    end do

    ! 穴埋め
    !
    xqz_SWF(:,jmax,:) = 0.0d0 
    xqz_PGrad(:,jmax,:) = 0.0d0 

  end subroutine PGrad_HE

  
  subroutine PGrad_VI(                       &
    &  Alpha, beta, dz, CpDry, xyr_VPTempBZ,  &
    &  imin, imax, jmin, jmax, kmin, kmax,   &  
    &  xyz_VelDiv, xyz_ExnerNs, xyz_ExnerAs, & !(IN)
    &  xyr_PGrad, xyr_SWF )                    !(OUT)
    
    implicit none

    real(DP), intent(in)  :: Alpha, beta, dz, CpDry
    integer,  intent(in)  :: imin, imax, jmin, jmax, kmin, kmax
    real(DP), intent(in)  :: xyr_VPTempBZ(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyz_VelDiv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyz_ExnerNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyz_ExnerAs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyr_PGrad(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyr_SWF(imin:imax,jmin:jmax,kmin:kmax)
    integer               :: i, j, k

    ! 音波減衰項
    !
    do k = kmin, kmax-1
      do j = jmin, jmax
        do i = imin, imax

          xyr_SWF(i,j,k) =                &
            & + Alpha                     &
            &   * (                       &
            &       xyz_VelDiv(i,j,k+1)   & 
            &     - xyz_VelDiv(i,j,k)     & 
            &     ) / dz

        end do
      end do
    end do

    ! 圧力傾度力
    !
    do k = kmin, kmax-1
      do j = jmin, jmax
        do i = imin, imax

          xyr_PGrad(i,j,k) =                 &
            & - CpDry * xyr_VPTempBZ(i,j,k)  &
            &   * (                          &
            &       beta                     &
            &       * (                      &
            &           xyz_ExnerAs(i,j,k+1) &
            &         - xyz_ExnerAs(i,j,k)   &
            &         )                      &
            &     + (1.0d0 - beta)           &
            &       * (                      &
            &           xyz_ExnerNs(i,j,k+1) &
            &         - xyz_ExnerNs(i,j,k)   &
            &         )                      &
            &     ) / dz

        end do
      end do
    end do

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

  end subroutine PGrad_VI



  subroutine Dynamics_Tendency_Output

    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='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


  subroutine Dynamics_Km_forcing(         &
    & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & ! (in)
    & xyz_KmBl,    xyz_KmNl,              & ! (in)
    & xyz_DKmDtNl                         & ! (inout)
    & )

    implicit none

    real(DP), intent(in)    :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax) 
    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) :: xyz_DKmDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)                :: xyz_Adv2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)                :: xyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax)

    character(15)        :: info2

    ! 移流および数値拡散
    !
    call AdvC4_nDiff_xyz(                    &
      &  xyz_KmBl, xyz_KmNl,                 & !(IN)
      &  pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN)
      &  xyz_Adv, xyz_nDiff )                  !(OUT)

    xyz_Adv2  =                                             &
      & - xyz_avr_pyz(pyz_VelXNl * pyz_c4dx_xyz(xyz_KmNl))  &
      & - xyz_avr_xqz(xqz_VelYNl * xqz_c4dy_xyz(xyz_KmNl))  &
      & - xyz_avr_xyr(xyr_VelZNl * xyr_c4dz_xyz(xyz_KmNl))    

    xyz_nDiff2 =                                                               &
      &  - NuHm * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyz_KmBl ))))) &
      &  - NuHm * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))))) &
      &  - NuVm * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))))) 


    info2 = "Km, Adv: "
    call CHK_Val( xyz_Adv(imin:imax,jmin:jmax,kmin:kmax), xyz_Adv2(imin:imax,jmin:jmax,kmin:kmax), info2 )
    info2 = "Km, nDiff: "
    call CHK_Val( xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax), xyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax), info2 )

    
    xyz_DKmDtNl = xyz_DKmDtNl + xyz_Adv + xyz_nDiff

    call HistoryAutoPut(TimeN, 'KmAdv',   xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'KmNDiff', xyz_nDiff(1:nx,1:ny,1:nz))
    
  end subroutine Dynamics_Km_forcing

  subroutine CHK_Val(xyz_Var1, xyz_Var2, info)

    real(DP), intent(in)      :: xyz_Var1(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)      :: xyz_Var2(imin:imax,jmin:jmax,kmin:kmax)
    character(15), intent(in) :: info

    integer  :: min1(3), min2(3), max1(3), max2(3)
    real(DP) :: tmp1(nx,ny,nz), tmp2(nx,ny,nz)

    real(DP), parameter :: EPS = 1.0d-12

    tmp1 = xyz_Var1(1:nx,1:ny,1:nz) - xyz_Var2(1:nx,1:ny,1:nz)
    tmp2 = (xyz_Var1(1:nx,1:ny,1:nz) - xyz_Var2(1:nx,1:ny,1:nz)) / (xyz_Var2(1:nx,1:ny,1:nz) + 1.0d-60)
    min1 = minloc(tmp1); min2 = minloc(tmp2)
    max1 = maxloc(tmp1); max2 = maxloc(tmp2)

    if (     abs(minval( tmp1 )) > EPS &
      & .OR. abs(maxval( tmp1 )) > EPS &
      & .OR. abs(minval( tmp2 )) > EPS &
      & .OR. abs(maxval( tmp2 )) > EPS ) then 

      write(*,*) "***", info
      write(*,*) minval( tmp1 ), maxval( tmp1 ), minval( tmp2 ), maxval( tmp2 )
      write(*,*) min1, tmp1(min1(1),min1(2),min1(3)), xyz_Var1(min1(1),min1(2),min1(3)), xyz_Var2(min1(1),min1(2),min1(3))
      write(*,*) max1, tmp1(max1(1),max1(2),max1(3)), xyz_Var1(max1(1),max1(2),max1(3)), xyz_Var2(max1(1),max1(2),max1(3))
      write(*,*) min2, tmp2(min2(1),min2(2),min2(3)), xyz_Var1(min2(1),min2(2),min2(3)), xyz_Var2(min2(1),min2(2),min2(3))
      write(*,*) max2, tmp2(max2(1),max2(2),max2(3)), xyz_Var1(max2(1),max2(2),max2(3)), xyz_Var2(max2(1),max2(2),max2(3))
    end if

  end subroutine CHK_Val

  subroutine CHK2_Val(xyz_Var1, xyz_Var2, info)

    real(DP), intent(in)      :: xyz_Var1(nx,ny,nz)
    real(DP), intent(in)      :: xyz_Var2(nx,ny,nz)
    character(15), intent(in) :: info

    real(DP), parameter :: EPS = 1.0d-12

    integer  :: min1(3), min2(3), max1(3), max2(3)
    real(DP) :: tmp1(nx,ny,nz), tmp2(nx,ny,nz)

    tmp1 = xyz_Var1(1:nx,1:ny,1:nz) - xyz_Var2(1:nx,1:ny,1:nz)
    tmp2 = (xyz_Var1(1:nx,1:ny,1:nz) - xyz_Var2(1:nx,1:ny,1:nz)) / (xyz_Var2(1:nx,1:ny,1:nz) + 1.0d-60)
    min1 = minloc(tmp1); min2 = minloc(tmp2)
    max1 = maxloc(tmp1); max2 = maxloc(tmp2)

    if (     abs(minval( tmp1 )) > EPS &
      & .OR. abs(maxval( tmp1 )) > EPS &
      & .OR. abs(minval( tmp2 )) > EPS &
      & .OR. abs(maxval( tmp2 )) > EPS ) then     

      write(*,*) "***", info
      write(*,*) minval( tmp1 ), maxval( tmp1 ), minval( tmp2 ), maxval( tmp2 )
      write(*,*) min1, tmp1(min1(1),min1(2),min1(3)), xyz_Var1(min1(1),min1(2),min1(3)), xyz_Var2(min1(1),min1(2),min1(3))
      write(*,*) max1, tmp1(max1(1),max1(2),max1(3)), xyz_Var1(max1(1),max1(2),max1(3)), xyz_Var2(max1(1),max1(2),max1(3))
      write(*,*) min2, tmp2(min2(1),min2(2),min2(3)), xyz_Var1(min2(1),min2(2),min2(3)), xyz_Var2(min2(1),min2(2),min2(3))
      write(*,*) max2, tmp2(max2(1),max2(2),max2(3)), xyz_Var1(max2(1),max2(2),max2(3)), xyz_Var2(max2(1),max2(2),max2(3))
    end if

  end subroutine CHK2_Val

  
end module DynamicsHEVI

