subroutine Dynamics2D_Long_forcing( pyz_VelXBl,  pyz_VelXNl, xyr_VelZBl,  xyr_VelZNl, xyz_PTempBl, xyz_PTempNl, xyzf_QMixBl, xyzf_QMixNl, xyz_KmBl,    xyz_KmNl, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, xyz_DPTempDtNl, xyzf_DQMixDtNl, xyz_DKmDtNl )

    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)    :: 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(in)    :: xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyz_KmNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(inout) :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(inout) :: xyz_DPTempDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xyzf_DQMixDtNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(inout) :: xyz_DKmDtNl(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)             :: pyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_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


    !--------------------------------------------------------------------
    ! 乱流拡散係数

    ! 移流および数値拡散
    !
    call AdvC4_nDiff_xyz( xyz_KmBl, xyz_KmNl ) !(IN)
    
    ! tendency の更新
    !
    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))
    
    !--------------------------------------------------------------------
    ! 温位

    ! 移流については, 基本場も考慮する.
    !
    xyz_PTempAll = xyz_PTempNl + xyz_PTempBZ
    call AdvC4_nDiff_xyz( xyz_PTempBl, xyz_PTempAll ) !(IN)

    ! 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)

    ! 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_xyr

    ! 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 = 0.0d0

    ! Buoyancy 
    ! 
    call BuoyancyLong_xyr

    ! 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))    

  contains
    
    subroutine AdvC4_nDiff_pyz_xyr

      implicit none

      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 = 1, ny
          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 - ( + ( 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,:,:) = 0.0d0
      pyz_Adv(imax-1:imax,:,:) = 0.0d0
      pyz_Adv(:,:,kmin:kmin+1) = 0.0d0
      pyz_Adv(:,:,kmax-1:kmax) = 0.0d0
      
      do k = kmin + 2, kmax - 2
        do j = 1, ny
          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 - 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,:,:) = 0.0d0
      xyr_Adv(imax-1:imax,:,:) = 0.0d0
      xyr_Adv(:,:,kmin:kmin+1) = 0.0d0
      xyr_Adv(:,:,kmax-1:kmax) = 0.0d0
      
      do k = kmin + 2, kmax - 2
        do j = 1, ny
          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,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,:,:) = 0.0d0
      pyz_nDiff(imax-1:imax,:,:) = 0.0d0
      pyz_nDiff(:,:,kmin:kmin+1) = 0.0d0
      pyz_nDiff(:,:,kmax-1:kmax) = 0.0d0

      do k = kmin + 2, kmax - 2
        do j = 1, ny          
          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,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,:,:) = 0.0d0
      xyr_nDiff(imax-1:imax,:,:) = 0.0d0
      xyr_nDiff(:,:,kmin:kmin+1) = 0.0d0
      xyr_nDiff(:,:,kmax-1:kmax) = 0.0d0
      
    end subroutine AdvC4_nDiff_pyz_xyr
    
    
    subroutine AdvC4_nDiff_xyz( xyz_VarBl, xyz_VarNl )
      
      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)              :: 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 = 1, ny
          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 - ( 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,:,:) = 0.0d0
      xyz_Adv(imax-1:imax,:,:) = 0.0d0
      xyz_Adv(:,:,kmin:kmin+1) = 0.0d0
      xyz_Adv(:,:,kmax-1:kmax) = 0.0d0
      
      ! 4 次の数値拡散: 2 次精度中心差分
      ! 
      do k = kmin + 2, kmax - 2
        do j = 1, ny
          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,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,:,:) = 0.0d0
      xyz_nDiff(imax-1:imax,:,:) = 0.0d0
      xyz_nDiff(:,:,kmin:kmin+1) = 0.0d0
      xyz_nDiff(:,:,kmax-1:kmax) = 0.0d0
      
    end subroutine AdvC4_nDiff_xyz
    
    
    subroutine AdvC4_nDiff_xyzf( xyzf_VarBl, xyzf_VarNl ) 

      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)              :: 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 = 1, ny
            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 - ( 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,:,:,:) = 0.0d0
      xyzf_Adv(imax-1:imax,:,:,:) = 0.0d0
      xyzf_Adv(:,:,kmin:kmin+1,:) = 0.0d0
      xyzf_Adv(:,:,kmax-1:kmax,:) = 0.0d0
      
      ! 数値拡散: 2 次精度中心差分
      ! 
      do s = 1, ncmax
        do k = kmin + 2, kmax - 2
          do j = 1, ny
            do i = imin + 2, imax - 2
              
              xyzf_nDiff(i,j,k,s) = - ( xyzf_VarBl(i+2,j,k,s) + xyzf_VarBl(i-2,j,k,s) - xyzf_VarBl(i+1,j,k,s) * 4.0d0 - xyzf_VarBl(i-1,j,k,s) * 4.0d0 + xyzf_VarBl(i  ,j,k,s) * 6.0d0 ) * NuHh / ( dx ** 4.0d0 ) - ( xyzf_VarBl(i,j,k+2,s) + xyzf_VarBl(i,j,k-2,s) - xyzf_VarBl(i,j,k+1,s) * 4.0d0 - xyzf_VarBl(i,j,k-1,s) * 4.0d0 + xyzf_VarBl(i,j,k  ,s) * 6.0d0 ) * NuVh / ( dz ** 4.0d0 )
            end do
          end do
        end do
      end do
      
      xyzf_nDiff(imin:imin+1,:,:,:) = 0.0d0
      xyzf_nDiff(imax-1:imax,:,:,:) = 0.0d0
      xyzf_nDiff(:,:,kmin:kmin+1,:) = 0.0d0
      xyzf_nDiff(:,:,kmax-1:kmax,:) = 0.0d0
      
    end subroutine AdvC4_nDiff_xyzf
  

    subroutine BuoyancyLong_xyr
      
      use composition, only: GasNum, IdxG, MolWtWet       ! 湿潤成分の分子量
      use constants,only:    MolWtDry, Grav           ! 重力加速度
      
      implicit none
      
      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
      
      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 = 1, ny
          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 = 1, ny
          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 = 1, ny
          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
      
    end subroutine BuoyancyLong_xyr
    
  end subroutine Dynamics2D_Long_forcing

  subroutine Dynamics2D_Short_integrate( pyz_VelXNs, xyr_VelZNs, xyz_ExnerNs, pyz_DVelXDtNl, xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs, pyz_VelXAs, xqz_VelYAs, xyr_VelZAs, xyz_ExnerAs )

    implicit none

    real(DP), intent(in)     :: pyz_VelXNs(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)     :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout)  :: xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax)
    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) :: 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) :: xyr_PGrad(imin:imax,jmin:jmax,kmin:kmax)

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


    !------------------------------------------------------------
    ! initialize: Divergence of velocity
    !
    call VelDivC2

    !------------------------------------------------------------
    ! VelX, VelY
    !  水平方向は陽解法で解く. 
    !
    call PGrad_HE

    ! tendency 
    !
    pyz_DVelXDtNs = pyz_PGrad + pyz_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))

    ! Time integration
    !
    pyz_VelXAs = pyz_VelXNs + DelTimeShort * (pyz_DVelXDtNl + pyz_DVelXDtNs)

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

    ! Y 方向には値がゼロ
    xqz_VelYAs = 0.0d0

    
    !------------------------------------------------------------
    ! Exner function
    !  積分値を返すことに注意. 
    !
    call Exner_HEVI

    !------------------------------------------------------------
    ! VelZ
    !
    call PGrad_VI

    ! 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)

  contains

    subroutine VelDivC2
      
      implicit none
      integer              :: i, j, k
      
      do  k = kmin + 1, kmax 
        do j = 1, ny
          do i = imin + 1, imax 
            xyz_VelDivNs(i,j,k) = + ( pyz_VelXNs(i,j,k) - pyz_VelXNs(i-1,j,k) ) / dx + ( xyr_VelZNs(i,j,k) - xyr_VelZNs(i,j,k-1) ) / dz
          end do
        end do
      end do
      
      xyz_VelDivNs(imin,:,:) = 0.0d0 
      xyz_VelDivNs(:,:,kmin) = 0.0d0 
      
    end subroutine VelDivC2


    subroutine PGrad_HE
      
      implicit none
      integer              :: i, j, k
      
      !------------------------------------------------------------------
      ! X 方向
      
      do k = kmin, kmax
        do j = 1, ny          
          do i = imin, imax - 1

            ! 音波減衰項
            !            
            pyz_SWF(i,j,k) = Alpha * ( xyz_VelDivNs(i+1,j,k) - xyz_VelDivNs(i,j,k) ) / dx            
            
            ! 圧力傾度力
            !
            pyz_PGrad(i,j,k) = - CpDry * pyz_VPTempBZ(i,j,k) * ( xyz_ExnerNs(i+1,j,k) - xyz_ExnerNs(i,j,k) ) / dx                                
          end do
        end do
      end do
      
      ! 穴埋め
      !
      pyz_SWF(imax,:,:)   = 0.0d0 
      pyz_PGrad(imax,:,:) = 0.0d0 
      
    end subroutine PGrad_HE


    subroutine PGrad_VI
    
      implicit none
      integer               :: i, j, k

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

            ! 音波減衰項
            !            
            xyr_SWF(i,j,k) = + Alpha * ( xyz_VelDivNs(i,j,k+1) - xyz_VelDivNs(i,j,k) ) / dz
            
            ! 圧力傾度力
            !
            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 Exner_HEVI
      !
      !陰解法を用いたエクスナー関数の計算. 
      !

      !暗黙の型宣言禁止
      implicit none

      !作業変数定義
      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)               :: F(1:nx,1:ny,1:nz)
      real(DP)               :: F0(1:nx,1:ny,kmin:kmax-1)  
      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'
      
      
      ! Initialize
      !
      dt = DelTimeShort

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

      !  添字の範囲は, 1:nx, 1:ny, 0:nz
      !  D(:,:,1) を求める時に D(:,:,0) の値が必要になるため. 
      !
      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) ) / dz + ( + Alpha * ( + xyz_VelDivNs(i,j,k+1) - xyz_VelDivNs(i,j,k) ) / dz + xyr_DVelZDtNl(i,j,k) ) / xyr_CpVPTempBZ(i,j,k)  
            
          end do
        end do
      end do
      
      ! 被微分関数
      !   配列 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) ) / 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 ) + 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
      
      !-----------------------------------------------------------
      !連立一次方程式の解を求める
      
      !変数の初期化
      !
      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(:,:,kmin:0) = 0.0d0
      xyz_ExnerAs(nx+1:imax,:,:) = 0.0d0
      xyz_ExnerAs(:,:,nz+1:kmax) = 0.0d0
      
      ! のり代に値を入れる
      !
      call SetMargin_xyz( xyz_ExnerAs ) ! (inout)
      
    end subroutine Exner_HEVI
    
  end subroutine Dynamics2D_Short_integrate

  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_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_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)  :: AlphaNDiff2 = 0.0d0  !2次の数値拡散の係数. CReSS マニュアルより
    real(DP)  :: NDiffRatio = 1.0d0   !速度に対する粘性を上げる場合は数字を 1 以上にする. 
    integer   :: unit            !装置番号
    integer   :: l

    NAMELIST /Dynamics_nml/ AlphaSound, AlphaNDiff, AlphaNDiff2, 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 = minval(x_dx)
    DelYMin = minval(y_dy)
    DelZMin = minval(z_dz)
    if (jmin == jmax) then 
       Alpha = AlphaSound * ( Min(DelXMin, DelZMin) ** 2.0d0 ) / DelTimeShort
    else
       Alpha = AlphaSound * ( Min(DelXMin, DelYMin, DelZMin) ** 2.0d0 ) / DelTimeShort
    end if

    if (jmin == jmax) then     
       NuHh4 = AlphaNDiff * ( DelXMin ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    else
       NuHh4 = AlphaNDiff * ( SQRT(DelXMin*DelYMin) ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    end if
    NuVh4 = AlphaNDiff * ( DelZMin ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    NuHm4 = NuHh4 * NDiffRatio
    NuVm4 = NuVh4 * NDiffRatio

    if (jmin == jmax) then     
       NuHh2 = AlphaNDiff2 * ( DelXMin ** 2.0d0 ) / (2.0d0 * DelTimeLong)
    else
       NuHh2 = AlphaNDiff2 * ( SQRT(DelXMin*DelYMin) ** 2.0d0 ) / (2.0d0 * DelTimeLong)
    end if
    NuVh2 = AlphaNDiff2 * ( DelZMin ** 2.0d0 ) / (2.0d0 * DelTimeLong)
    NuHm2 = NuHh2 * NDiffRatio
    NuVm2 = NuVh2 * NDiffRatio

    if (myrank == 0) then 
      if ( AlphaNDiff > 0.0d0 .AND. AlphaNDiff2 > 0.0d0 ) then 
        write(*,*) "CHECK! AlphaNDiff == 0.0d0 .OR. AlphaNDiff2 == 0.0d0"
        stop 
      end if
      call MessageNotify( "M", module_name, "Alpha = %f", d=(/Alpha/) )
      call MessageNotify( "M", module_name, "NuHh = %f", d=(/NuHh4/) )
      call MessageNotify( "M", module_name, "NuVh = %f", d=(/NuVh4/) )
      call MessageNotify( "M", module_name, "NuHm = %f", d=(/NuHm4/) )
      call MessageNotify( "M", module_name, "NuVm = %f", d=(/NuVm4/) )
      call MessageNotify( "M", module_name, "NuHh = %f", d=(/NuHh2/) )
      call MessageNotify( "M", module_name, "NuVh = %f", d=(/NuVh2/) )
      call MessageNotify( "M", module_name, "NuHm = %f", d=(/NuHm2/) )
      call MessageNotify( "M", module_name, "NuVm = %f", d=(/NuVm2/) )
      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 DynamicsVI_init()

    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='PTempNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of potential temperature (2 order)', 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='ExnerNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of exner function (2 order)', 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')

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

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

      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_NDiff2', dims=(/'x','y','z','t'/), longname='Diffusion term of ' //trim(SpcWetSymbol(l))//' mixing ratio (2 order)', 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='VelXNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of velocity (x) (2 order)', 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='VelXTndNs', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='Velocity Tendency (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='VelYNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of velocity (y) (2 order)', 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='VelYTndNs', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='Velocity Tendency (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='VelZNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of Velocity (z) (2 order)', 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')

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

  end subroutine Dynamics_Init

  subroutine Dynamics_Km_forcing( pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, xyz_KmBl,    xyz_KmNl, xyz_DKmDtNl )

    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, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, 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 Dynamics_Km_forcing( pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, xyz_KmBl,    xyz_KmNl, xyz_DKmDtNl )

    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_Orig(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_NDiff2(imin:imax,jmin:jmax,kmin:kmax)

    !----------------------------------
    ! 拡散係数
    !

    ! Initialize
    !
    xyz_Orig = xyz_DKmDtNl

    ! Advection term
    !
    xyz_Adv  = - 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))    

    ! Numerical diffusion term 
    !
    xyz_NDiff4 = - NuHm4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyz_KmBl ))))) - NuHm4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))))) - NuVm4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))))) 

    xyz_NDiff2 = + NuHm2 * (xyz_dx_pyz(pyz_dx_xyz( xyz_KmBl ))) + NuHm2 * (xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))) + NuVm2 * (xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))) 
    
    xyz_DKmDtNl = xyz_Orig + xyz_Adv + xyz_NDiff4 + xyz_NDiff2

    call HistoryAutoPut(TimeN, 'KmAdv',    xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'KmNDiff',  xyz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'KmNDiff2', xyz_NDiff2(1:nx,1:ny,1:nz))
    
    ! Set Margin
    !
!    call SetMargin_xyz( xyz_DKmDtNl )

  end subroutine Dynamics_Km_forcing

  subroutine Dynamics_Long_forcing( pyz_VelXBl,  pyz_VelXNl, xqz_VelYBl,  xqz_VelYNl, xyr_VelZBl,  xyr_VelZNl, xyz_PTempBl, xyz_PTempNl, xyzf_QMixBl, xyzf_QMixNl, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, xyz_DPTempDtNl, xyzf_DQMixDtNl, xyz_PTempAl, xyzf_QMixAl )

    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, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, 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, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, 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, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, pyz_Adv,    xqz_Adv,    xyr_Adv, 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, 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_forcing( pyz_VelXBl,  pyz_VelXNl, xqz_VelYBl,  xqz_VelYNl, xyr_VelZBl,  xyr_VelZNl, xyz_PTempBl, xyz_PTempNl, xyzf_QMixBl, xyzf_QMixNl, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, xyz_DPTempDtNl, xyzf_DQMixDtNl, xyz_PTempAl, xyzf_QMixAl )

    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), 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) :: xyz_DCDensDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_Orig(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_NDiff2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_Orig(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_NDiff2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_Orig(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_NDiff2(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_Orig(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_NDiff2(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyzf_Orig(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_Adv(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_NDiff4(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_NDiff2(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, 1:GasNum)
    integer              :: f

    !------------------------------
    ! tendency of cloud density    
    ! 

    ! initialize 
!    xyz_Orig = xyz_DCDensDtNl

    ! フラックス項の計算. 4 次精度中心差分を用いて計算
    !
!    xyz_Adv = &
!      & - xyz_c4dx_pyz(pyz_VelXNl * pyz_avr_xyz(xyz_CDensNl)) &
!      & - xyz_c4dy_xqz(xqz_VelYNl * xqz_avr_xyz(xyz_CDensNl)) &
!      & - xyz_c4dz_xyr(xyr_VelZNl * xyr_avr_xyz(xyz_CDensNl)) 

    ! 数値粘性項の計算
!    xyz_NDiff4 = &
!      &  - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_CDensBl))))) &
!      &  - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_CDensBl))))) &
!      &  - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_CDensBl))))) 

!    xyz_NDiff2 = &
!      &  + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz(xyz_CDensBl))) &
!      &  + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz(xyz_CDensBl))) &
!      &  + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz(xyz_CDensBl))) 
    
!    xyz_DCDensDtNl = xyz_Orig + xyz_Adv + xyz_NDiff4 + xyz_NDiff2

!    call HistoryAutoPut(TimeN, 'CDensAdv',  xyz_Adv(1:nx,1:ny,1:nz))
!    call HistoryAutoPut(TimeN, 'CDensNDiff', xyz_NDiff4(1:nx,1:ny,1:nz))
!    call HistoryAutoPut(TimeN, 'CDensNDiff2', xyz_NDiff2(1:nx,1:ny,1:nz))


    ! tendency of potential temperature
    ! 

    ! initialize 
    xyz_Orig = xyz_DPTempDtNl
    xyz_PTempAll = xyz_PTempNl + xyz_PTempBZ

    ! Advection term
    ! xyz_AdvScalar( xyz_PTempNl + xyz_PTempBZ, pyz_VelXNl, pyz_VelXNl, xyr_VelZNl) 
    !
    xyz_Adv = - 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))  

    ! numerical diffusion term
    ! xyz_Num = xyz_NumDiffScalar( xyz_PTempBl)
    !
    xyz_NDiff4 = - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))))) - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))))) - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_PTempBl ))))) 

    xyz_NDiff2 = + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))) + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))) + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz( xyz_PTempBl ))) 

    ! sum
    !
    xyz_DPTempDtNl = xyz_Orig + xyz_Adv + xyz_NDiff4 + xyz_NDiff2
    
    ! output
    !
    call HistoryAutoPut(TimeN, 'PTempAdv',   xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'PTempNDiff', xyz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'PTempNDiff2', xyz_NDiff2(1:nx,1:ny,1:nz))

    xyz_PTempAl = xyz_PTempBl + (2.0d0 * DelTimeLong) * xyz_DPTempDtNl

    ! Set Margin
    !
    call SetMargin_xyz(xyz_PTempAl)

    
    !------------------------------
    ! tendency of mixing ratio
    ! 

    ! initialize
    xyzf_Orig = xyzf_DQMixDtNl
    xyzf_QMixAll = xyzf_QMixNl + xyzf_QMixBZ

    do f = 1, ncmax
      ! Advection term
      !xyzf_Adv  = xyzf_AdvScalar(xyzf_QMixNl + xyzf_QMixBZ, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
      !
      xyzf_Adv(:,:,:,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)))    

      ! numerical diffusion term
      ! xyzf_Diff = xyzf_NumDiffScalar(xyzf_QMixBl) 
      !
      xyzf_NDiff4(:,:,:,f) = - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))))) - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))))) - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyzf_QMixBl(:,:,:,f) )))))

      xyzf_NDiff2(:,:,:,f) = + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))) + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))) + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz( xyzf_QMixBl(:,:,:,f) )))
    end do

    ! sum
    !
    xyzf_DQMixDtNl = xyzf_Orig + xyzf_Adv + xyzf_NDiff4 + xyzf_NDiff2

    ! 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_NDiff4(1:nx,1:ny,1:nz,f))
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_NDiff2', xyzf_NDiff2(1:nx,1:ny,1:nz,f))
    end do

    ! time integration
    !
    xyzf_QMixAl = xyzf_QMixBl + (2.0d0 * DelTimeLong) * xyzf_DQMixDtNl

    ! Set Margin
    ! 
    call SetMargin_xyzf(xyzf_QMixAl)

    ! 負の値を埋める
    !
    call FillNegativeQMix(xyzf_QMixAl)

    ! Set Margin
    ! 
    call SetMargin_xyzf(xyzf_QMixAl)

    !------------------------------
    ! tendency of VelX
    ! 

    ! initializa
    !
    pyz_Orig = pyz_DVelXDtNl

    ! Advection term
    !pyz_Adv  = pyz_AdvVelX(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) 
    !
    pyz_Adv  = - 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 ) )

    ! Numerical diffusion term 
    !pyz_Diff = pyz_NumDiffVelX(pyz_VelXBl)
    pyz_NDiff4 = - NuHm4 * ( pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz( pyz_VelXBl ))))) - NuHm4 * ( pyz_dy_pqz(pqz_dy_pyz(pyz_dy_pqz(pqz_dy_pyz( pyz_VelXBl ))))) - NuVm4 * ( pyz_dz_pyr(pyr_dz_pyz(pyz_dz_pyr(pyr_dz_pyz( pyz_VelXBl )))))

    pyz_NDiff2 = + NuHm2 * ( pyz_dx_xyz(xyz_dx_pyz( pyz_VelXBl ))) + NuHm2 * ( pyz_dy_pqz(pqz_dy_pyz( pyz_VelXBl ))) + NuVm2 * ( pyz_dz_pyr(pyr_dz_pyz( pyz_VelXBl )))

    ! sum
    !
    pyz_DVelXDtNl = pyz_Orig + pyz_Adv + pyz_NDiff4 + pyz_NDiff2

    call HistoryAutoPut(TimeN, 'VelXAdv',   pyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXNDiff', pyz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXNDiff2', pyz_NDiff2(1:nx,1:ny,1:nz))

    !------------------------------
    ! tendency of VelY
    !     

    ! ininitalize
    xqz_Orig = xqz_DVelYDtNl

    ! Advection term
    ! xqz_Adv  = xqz_AdvVelY(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
    !
    xqz_Adv  = - 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 ) )

    ! Numerical diffusion term
    ! xqz_Diff = xqz_NumDiffVelY(xqz_VelYBl)
    !
    xqz_NDiff4 = - NuHm4 * ( xqz_dx_pqz(pqz_dx_xqz(xqz_dx_pqz(pqz_dx_xqz( xqz_VelYBl ))))) - NuHm4 * ( xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz( xqz_VelYBl ))))) - NuVm4 * ( xqz_dz_xqr(xqr_dz_xqz(xqz_dz_xqr(xqr_dz_xqz( xqz_VelYBl )))))

    xqz_NDiff2 = + NuHm2 * ( xqz_dx_pqz(pqz_dx_xqz( xqz_VelYBl ))) + NuHm2 * ( xqz_dy_xyz(xyz_dy_xqz( xqz_VelYBl ))) + NuVm2 * ( xqz_dz_xqr(xqr_dz_xqz( xqz_VelYBl )))

    ! sum
    !
    xqz_DVelYDtNl = xqz_Orig + xqz_Adv + xqz_NDiff4 + xqz_NDiff2
    
    call HistoryAutoPut(TimeN, 'VelYAdv',   xqz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYNDiff', xqz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYNDiff2', xqz_NDiff2(1:nx,1:ny,1:nz))

    !------------------------------
    ! tendency of VelZ
    ! 

    ! Initialization
    !
    xyr_Orig = xyr_DVelZDtNl
    
    do f = 1, GasNum
      xyzf_QMixPerMolWt(:,:,:,f) = xyzf_QMixNl(:,:,:,IdxG(f)) / MolWtWet(IdxG(f))
    end do

    ! Advection term
    ! xyr_Adv  = xyr_AdvVelZ(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
    !
    xyr_Adv  = - 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 ) )

    ! Numerical diffusion term
    !xyr_Diff = xyr_NumDiffVelZ(xyr_VelZBl)
    !
    xyr_NDiff4 = - NuHm4 * ( xyr_dx_pyr(pyr_dx_xyr(xyr_dx_pyr(pyr_dx_xyr( xyr_VelZBl ))))) - NuHm4 * ( xyr_dy_xqr(xqr_dy_xyr(xyr_dy_xqr(xqr_dy_xyr( xyr_VelZBl ))))) - NuVm4 * ( xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr( xyr_VelZBl )))))

    xyr_NDiff2 = + NuHm2 * ( xyr_dx_pyr(pyr_dx_xyr( xyr_VelZBl ))) + NuHm2 * ( xyr_dy_xqr(xqr_dy_xyr( xyr_VelZBl ))) + NuVm2 * ( xyr_dz_xyz(xyz_dz_xyr( xyr_VelZBl )))

    ! Buoyancy due to temperature disturbunce
    !xyr_BuoyT = xyr_Buoy(xyz_PTempNl)
    !
    xyr_BuoyT = Grav * xyr_avr_xyz( xyz_PTempNl / xyz_PTempBZ)

    ! Buoyancy due to molecular weight
    !
    xyr_BuoyM = + 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 ) 

    ! Buoyancy due to loading
    !
    xyr_BuoyD = - Grav * xyr_avr_xyz( sum(xyzf_QMixNl(:,:,:,GasNum+1:ncmax), 4) ) / ( 1.0d0 + xyr_QMixBZ )

    ! sum
    !
    xyr_DVelZDtNl = xyr_Orig + xyr_Adv + xyr_NDiff4 + xyr_NDiff2 + 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_NDiff4(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'VelZNDiff2', xyr_NDiff2(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))    

    ! Set Margin
    !
!    call SetMargin_pyz( pyz_DVelXDtNl )
!    call SetMargin_xqz( xqz_DVelYDtNl )
!    call SetMargin_xyr( xyr_DVelZDtNl )
!    call SetMargin_xyz( xyz_DPTempDtNl )
!    call SetMargin_xyz( xyz_DCDensDtNl )
!    call SetMargin_xyzf(xyzf_DQMixDtNl )
 
  end subroutine Dynamics_Long_forcing

  subroutine Dynamics_Long_forcing( pyz_VelXBl,  pyz_VelXNl, xqz_VelYBl,  xqz_VelYNl, xyr_VelZBl,  xyr_VelZNl, xyz_PTempBl, xyz_PTempNl, xyzf_QMixBl, xyzf_QMixNl, xyz_KmBl,    xyz_KmNl, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, xyz_DPTempDtNl, xyzf_DQMixDtNl, xyz_DKmDtNl )

    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(in)    :: xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)    :: xyz_KmNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(inout) :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(inout) :: xyz_DPTempDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xyzf_DQMixDtNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(inout) :: xyz_DKmDtNl(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)             :: pyz_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


    !--------------------------------------------------------------------
    ! 乱流拡散係数

    ! 移流および数値拡散
    !
    call AdvC4_nDiff_xyz( xyz_KmBl, xyz_KmNl ) !(IN)
    
    ! tendency の更新
    !
    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))

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

    ! 移流については, 基本場も考慮する.
    !
    xyz_PTempAll = xyz_PTempNl + xyz_PTempBZ
    call AdvC4_nDiff_xyz( xyz_PTempBl, xyz_PTempAll ) !(IN)

    ! 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)

    ! 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

    ! 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

    ! 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))    

  contains
    
    subroutine AdvC4_nDiff_pyz_xqz_xyr

      implicit none

      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,:,:) = 0.0d0
      pyz_Adv(imax-1:imax,:,:) = 0.0d0
      pyz_Adv(:,jmin:jmin+1,:) = 0.0d0
      pyz_Adv(:,jmax-1:jmax,:) = 0.0d0
      pyz_Adv(:,:,kmin:kmin+1) = 0.0d0
      pyz_Adv(:,:,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,:,:) = 0.0d0
      xqz_Adv(imax-1:imax,:,:) = 0.0d0
      xqz_Adv(:,jmin:jmin+1,:) = 0.0d0
      xqz_Adv(:,jmax-1:jmax,:) = 0.0d0
      xqz_Adv(:,:,kmin:kmin+1) = 0.0d0
      xqz_Adv(:,:,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,:,:) = 0.0d0
      xyr_Adv(imax-1:imax,:,:) = 0.0d0
      xyr_Adv(:,jmin:jmin+1,:) = 0.0d0
      xyr_Adv(:,jmax-1:jmax,:) = 0.0d0
      xyr_Adv(:,:,kmin:kmin+1) = 0.0d0
      xyr_Adv(:,:,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,:,:) = 0.0d0
      pyz_nDiff(imax-1:imax,:,:) = 0.0d0
      pyz_nDiff(:,jmin:jmin+1,:) = 0.0d0
      pyz_nDiff(:,jmax-1:jmax,:) = 0.0d0
      pyz_nDiff(:,:,kmin:kmin+1) = 0.0d0
      pyz_nDiff(:,:,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,:,:) = 0.0d0
      xqz_nDiff(imax-1:imax,:,:) = 0.0d0
      xqz_nDiff(:,jmin:jmin+1,:) = 0.0d0
      xqz_nDiff(:,jmax-1:jmax,:) = 0.0d0
      xqz_nDiff(:,:,kmin:kmin+1) = 0.0d0
      xqz_nDiff(:,:,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,:,:) = 0.0d0
      xyr_nDiff(imax-1:imax,:,:) = 0.0d0
      xyr_nDiff(:,jmin:jmin+1,:) = 0.0d0
      xyr_nDiff(:,jmax-1:jmax,:) = 0.0d0
      xyr_nDiff(:,:,kmin:kmin+1) = 0.0d0
      xyr_nDiff(:,:,kmax-1:kmax) = 0.0d0
      
    end subroutine AdvC4_nDiff_pyz_xqz_xyr
    
    
    subroutine AdvC4_nDiff_xyz( xyz_VarBl, xyz_VarNl )
      
      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)              :: 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,:,:) = 0.0d0
      xyz_Adv(imax-1:imax,:,:) = 0.0d0
      xyz_Adv(:,jmin:jmin+1,:) = 0.0d0
      xyz_Adv(:,jmax-1:jmax,:) = 0.0d0
      xyz_Adv(:,:,kmin:kmin+1) = 0.0d0
      xyz_Adv(:,:,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,:,:) = 0.0d0
      xyz_nDiff(imax-1:imax,:,:) = 0.0d0
      xyz_nDiff(:,jmin:jmin+1,:) = 0.0d0
      xyz_nDiff(:,jmax-1:jmax,:) = 0.0d0
      xyz_nDiff(:,:,kmin:kmin+1) = 0.0d0
      xyz_nDiff(:,:,kmax-1:kmax) = 0.0d0
      
    end subroutine AdvC4_nDiff_xyz
    
    
    subroutine AdvC4_nDiff_xyzf( xyzf_VarBl, xyzf_VarNl ) 

      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)              :: 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,:,:,:) = 0.0d0
      xyzf_Adv(imax-1:imax,:,:,:) = 0.0d0
      xyzf_Adv(:,jmin:jmin+1,:,:) = 0.0d0
      xyzf_Adv(:,jmax-1:jmax,:,:) = 0.0d0
      xyzf_Adv(:,:,kmin:kmin+1,:) = 0.0d0
      xyzf_Adv(:,:,kmax-1:kmax,:) = 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-2,j,k,s) - xyzf_VarBl(i+1,j,k,s) * 4.0d0 - xyzf_VarBl(i-1,j,k,s) * 4.0d0 + xyzf_VarBl(i  ,j,k,s) * 6.0d0 ) * NuHh / ( dx ** 4.0d0 ) - ( xyzf_VarBl(i,j+2,k,s) + xyzf_VarBl(i,j-2,k,s) - xyzf_VarBl(i,j+1,k,s) * 4.0d0 - xyzf_VarBl(i,j-1,k,s) * 4.0d0 + xyzf_VarBl(i,j  ,k,s) * 6.0d0 ) * NuHh / ( dy ** 4.0d0 ) - ( xyzf_VarBl(i,j,k+2,s) + xyzf_VarBl(i,j,k-2,s) - xyzf_VarBl(i,j,k+1,s) * 4.0d0 - xyzf_VarBl(i,j,k-1,s) * 4.0d0 + xyzf_VarBl(i,j,k  ,s) * 6.0d0 ) * NuVh / ( dz ** 4.0d0 )
            end do
          end do
        end do
      end do
      
      xyzf_nDiff(imin:imin+1,:,:,:) = 0.0d0
      xyzf_nDiff(imax-1:imax,:,:,:) = 0.0d0
      xyzf_nDiff(:,jmin:jmin+1,:,:) = 0.0d0
      xyzf_nDiff(:,jmax-1:jmax,:,:) = 0.0d0
      xyzf_nDiff(:,:,kmin:kmin+1,:) = 0.0d0
      xyzf_nDiff(:,:,kmax-1:kmax,:) = 0.0d0
      
    end subroutine AdvC4_nDiff_xyzf
  

    subroutine BuoyancyLong_xyr
      
      use composition, only: GasNum, IdxG, MolWtWet       ! 湿潤成分の分子量
      use constants,only:    MolWtDry, Grav           ! 重力加速度
      
      implicit none
      
      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
      
      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
      
    end subroutine BuoyancyLong_xyr
        
  end subroutine Dynamics_Long_forcing

  subroutine Dynamics_Long_integrate( xyz_VarBl, xyz_DVarDtNl, xyzf_VarBl, xyzf_DVarDtNl, xyz_VarAl, xyzf_VarAl )

    implicit none

    real(DP), intent(in)   :: xyz_VarBl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(in)   :: xyz_DVarDtNl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(in)   :: xyzf_VarBl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP), intent(in)   :: xyzf_DVarDtNl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)
    real(DP), intent(out)  :: xyz_VarAl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(out)  :: xyzf_VarAl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax)

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

    ! Set Margin
    !
    call SetMargin_xyz(xyz_VarAl)

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

    call SetMargin_xyzf(xyzf_VarAl)

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

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

  end subroutine Dynamics_Long_integrate

  subroutine Dynamics_Short_forcing( pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, xyz_ExnerNs, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs, pyz_VelXAs, xqz_VelYAs, xyr_VelZAs, xyz_ExnerAs )

    real(DP), intent(in)  :: pyz_VelXNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xqz_VelYNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyz_ExnerNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax)
!    real(DP), intent(in)  :: xyz_DExnerDtNs(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) :: xyz_VelLaplaNs(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: pyz_VelXPGrad(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_VelYPGrad(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyr_VelZPGrad(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: pyz_VelXSWF(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_VelYSWF(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyr_VelZSWF(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: xyz_ExnerAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_NDiff2(imin:imax,jmin:jmax,kmin:kmax)

    ! initialize: Divergence of velocity
    !
    xyz_VelDivNs = xyz_dx_pyz( pyz_VelXNs ) + xyz_dy_xqz( xqz_VelYNs ) + xyz_dz_xyr( xyr_VelZNs )
!    call HistoryAutoPut(TimeN, 'VelDiv', xyz_VelDivNs(1:nx,1:ny,1:nz))

!    xyz_VelLaplaNs = pyz_dx_xyz( xyz_VelDivNs ) +  xqz_dy_xyz( xyz_VelDivNs ) +  xyr_dz_xyz( xyz_VelDivNs )
!    call HistoryAutoPut(TimeN, 'VelLapla', xyz_VelLaplaNs(1:nx,1:ny,1:nz))

    !--------------------------------------
    ! VelX
    !
    pyz_VelXSWF   =   Alpha * pyz_dx_xyz( xyz_VelDivNs ) 
    pyz_VelXPGrad = - pyz_avr_xyz( CpDry * xyz_VPTempBZ ) * pyz_dx_xyz( xyz_ExnerNs ) + pyz_VelXSWF
    pyz_DVelXDtNs =   pyz_VelXPGrad 

    ! Time integration
    !
    pyz_VelXAs    = pyz_VelXNs + DelTimeShort * (pyz_DVelXDtNl + pyz_DVelXDtNs)

    call HistoryAutoPut(TimeN, 'VelXPGrad', pyz_VelXPGrad(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXSWF',   pyz_VelXSWF(1:nx,1:ny,1:nz))
!    call HistoryAutoPut(TimeN, 'VelXTndNs', pyz_DVelXDtNs(1:nx,1:ny,1:nz))

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

    !--------------------------------------
    ! VelY
    !    
    xqz_VelYSWF   =   Alpha * xqz_dy_xyz( xyz_VelDivNs ) 
    xqz_VelYPGrad = - xqz_avr_xyz( CpDry * xyz_VPTempBZ ) * xqz_dy_xyz( xyz_ExnerNs ) + xqz_VelYSWF
    xqz_DVelYDtNs =   xqz_VelYPGrad
    
    ! Time integration
    !
    xqz_VelYAs = xqz_VelYNs + DelTimeShort * (xqz_DVelYDtNl + xqz_DVelYDtNs)

    call HistoryAutoPut(TimeN, 'VelYPGrad', xqz_VelYPGrad(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYSWF',   xqz_VelYSWF(1:nx,1:ny,1:nz))
!    call HistoryAutoPut(TimeN, 'VelYTndNs', xqz_DVelYDtNs(1:nx,1:ny,1:nz))

    ! Set Margin
    !
    call SetMargin_xqz( xqz_VelYAs ) ! (inout)
    
    !--------------------------------------
    ! Exner function
    !

    ! フラックス項の計算. 4 次精度中心差分を用いて計算
    !
    xyz_Adv = - xyz_avr_pyz(pyz_VelXNs * pyz_c4dx_xyz(xyz_ExnerNs)) - xyz_avr_xqz(xqz_VelYNs * xqz_c4dy_xyz(xyz_ExnerNs)) - xyz_avr_xyr(xyr_VelZNs * xyr_c4dz_xyz(xyz_ExnerNs))  !&  
!      & + CpDry / CvDry * GasRDry * xyz_ExnerNs * xyz_VelDivNs

    ! 数値粘性項の計算
    xyz_NDiff4 = - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_ExnerNs))))) - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_ExnerNs))))) - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_ExnerNs))))) 

    xyz_NDiff2 = + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz(xyz_ExnerNs))) + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz(xyz_ExnerNs))) + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz(xyz_ExnerNs))) 

    ! 短い時間ステップでのtendency 
    xyz_DExnerDtNs = xyz_DExnerDtNs + xyz_Adv + xyz_NDiff4 + xyz_NDiff2

    !!!!!!!!!!!!!!!!!!!!!!!
    !
    ! test 
    !
    !!!!!!!!!!!!!!!!!!!!!!!
    xyz_DExnerDtNs = 0.0d0

    call HistoryAutoPut(TimeN, 'ExnerAdv',   xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerNDiff', xyz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerNDiff2', xyz_NDiff2(1:nx,1:ny,1:nz))
    
    xyz_ExnerAs = xyz_Exner( pyz_VelXAs, xqz_VelYAs, xyr_VelZNs, xyz_VelDivNs, xyz_ExnerNs, xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs )

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

!    write(*,*) "+++ Exner +++", xyz_ExnerAs(imin:imax,1,kmin:kmax)

    !--------------------------------------
    ! VelZ
    !
    xyr_VelZSWF =  Alpha * xyr_dz_xyz( xyz_VelDivNs ) 
    xyr_VelZPGrad = - xyr_avr_xyz(CpDry * xyz_VPTempBZ ) * ( beta * xyr_dz_xyz( xyz_ExnerAs ) + (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs ) ) + xyr_VelZSWF
    xyr_DVelZDtNs = xyr_VelZPGrad 

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

    call HistoryAutoPut(TimeN, 'VelZPGrad', xyr_VelZPGrad(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelZSWF',   xyr_VelZSWF(1:nx,1:ny,1:nz))
!    call HistoryAutoPut(TimeN, 'VelZTndNs', xyr_DVelZDtNs(1:nx,1:ny,1:nz))

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

  end subroutine Dynamics_Short_forcing

  subroutine Dynamics_Short_forcing( pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, xyz_ExnerNs, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs, pyz_VelXAs, xqz_VelYAs, xyr_VelZAs, xyz_ExnerAs )

    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, 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, pyz_PGrad, pyz_SWF, 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, xyz_VelDivNs, xyz_ExnerNs, xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_ExnerAs )

    !------------------------------------------------------------
    ! VelZ
    !
    call PGrad_VI( Alpha, beta, dz, CpDry,  xyr_VPTempBZ, imin, imax, jmin, jmax, kmin, kmax, xyz_VelDivNs, xyz_ExnerNs, xyz_ExnerAs, 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_Short_integrate( pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, xyz_ExnerNs, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs, pyz_VelXAs, xqz_VelYAs, xyr_VelZAs, xyz_ExnerAs )

    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)
    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)


    !------------------------------------------------------------     
    ! initialize: Divergence of velocity
    !
    call VelDivC2

    !------------------------------------------------------------
    ! VelX, VelY
    !  水平方向は陽解法で解く. 
    !
    call PGrad_HE

    ! 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

    !------------------------------------------------------------
    ! VelZ
    !
    call PGrad_VI

    ! 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)

  contains

    subroutine VelDivC2
      
      implicit none
      integer              :: i, j, k
      
      do  k = kmin + 1, kmax 
        do j = jmin + 1, jmax 
          do i = imin + 1, imax 
            xyz_VelDivNs(i,j,k) = + ( pyz_VelXNs(i,j,k) - pyz_VelXNs(i-1,j,k) ) / dx + ( xqz_VelYNs(i,j,k) - xqz_VelYNs(i,j-1,k) ) / dy + ( xyr_VelZNs(i,j,k) - xyr_VelZNs(i,j,k-1) ) / dz
          end do
        end do
      end do
      
      xyz_VelDivNs(imin,:,:) = 0.0d0 
      xyz_VelDivNs(:,jmin,:) = 0.0d0  
      xyz_VelDivNs(:,:,kmin) = 0.0d0 
      
    end subroutine VelDivC2


    subroutine PGrad_HE
    
      implicit none
      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_VelDivNs(i+1,j,k) - xyz_VelDivNs(i,j,k) ) / dx            
            
            ! 圧力傾度力
            !
            pyz_PGrad(i,j,k) = - CpDry * pyz_VPTempBZ(i,j,k) * ( xyz_ExnerNs(i+1,j,k) - xyz_ExnerNs(i,j,k) ) / dx                                
          end do
        end do
      end do
      
      ! 穴埋め
      !
      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_VelDivNs(i,j+1,k) - xyz_VelDivNs(i,j,k) ) / dy
               
            ! 圧力傾度力
            !             
            xqz_PGrad(i,j,k) = - CpDry * xqz_VPTempBZ(i,j,k) * ( xyz_ExnerNs(i,j+1,k) - xyz_ExnerNs(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
    
      implicit none
      integer               :: i, j, k

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

            ! 音波減衰項
            !            
            xyr_SWF(i,j,k) = + Alpha * ( xyz_VelDivNs(i,j,k+1) - xyz_VelDivNs(i,j,k) ) / dz
            
            ! 圧力傾度力
            !
            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 Exner_HEVI
      !
      !陰解法を用いたエクスナー関数の計算. 
      !

      !暗黙の型宣言禁止
      implicit none

      !作業変数定義
      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)               :: F(1:nx,1:ny,1:nz)
      real(DP)               :: F0(1:nx,1:ny,kmin:kmax-1)  
      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'
            
      ! Initialize
      !
      dt = DelTimeShort

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

      !  添字の範囲は, 1:nx, 1:ny, 0:nz
      !  D(:,:,1) を求める時に D(:,:,0) の値が必要になるため. 
      !
      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) ) / dz + ( + Alpha * ( + xyz_VelDivNs(i,j,k+1) - xyz_VelDivNs(i,j,k) ) / dz + xyr_DVelZDtNl(i,j,k) ) / xyr_CpVPTempBZ(i,j,k)  
            
          end do
        end do
      end do
      
      ! 被微分関数
      !   配列 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) ) / 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
      
      !-----------------------------------------------------------
      !連立一次方程式の解を求める
      
      !変数の初期化
      !
      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
    
  end subroutine Dynamics_Short_integrate