subroutine NumDiffusion_init()
    !
    ! NumDiffusion モジュールの初期化ルーチン
    !

    !暗黙の型宣言禁止
    implicit none
    
    ! CReSS マニュアルでは, 2 次精度中心差分の場合, Alpha < 1/8 くらいが適当と述べている. 
    ! deepconv では NuH として 500 以上の値が入っていたそうだ,
    NuH  = Alpha * ( DelX ** 2.0d0 ) / DelTimeLong
    NuV  = Alpha * ( DelZ ** 2.0d0 ) / DelTimeLong
    
    !確認
    call MessageNotify( "M", "NumDiffusion_init", "NuH = %f", d=(/NuH/) )
    call MessageNotify( "M", "NumDiffusion_init", "NuV = %f", d=(/NuV/) )

  end subroutine NumDiffusion_init

  function pz_NumDiffVelX(pz_VarX)
    !
    ! z 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in) :: pz_VarX(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !物理量
    real(8)             :: pz_NumDiffVelX(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !数値拡散
    
!    pz_NumDiffVelX = 0.0d0  !初期化
    pz_NumDiffVelX = NuH * ( pz_dx_xz( xz_dx_pz( pz_VarX ) ) ) + NuV * ( pz_dz_pr( pr_dz_pz( pz_VarX ) ) )
    
  end function pz_NumDiffVelX

  function xr_NumDiffVelZ(xr_VarZ)
    !
    ! x 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in) :: xr_VarZ(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !物理量
    real(8)             :: xr_NumDiffVelZ(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !数値拡散

!    xr_NumDiffVelZ = 0.0d0 
    xr_NumDiffVelZ = NuH * ( xr_dx_pr( pr_dx_xr( xr_VarZ ) ) ) + NuV * ( xr_dz_xz( xz_dz_xr( xr_VarZ ) ) )
    
  end function xr_NumDiffVelZ

  function xz_NumDiffKm(xz_Scalar)
    !
    ! x, z 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !スカラー量
    real(8)             :: xz_NumDiffKm(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !水平方向の数値拡散
    
    xz_NumDiffKm = NuH * (xz_dx_pz(pz_dx_xz( xz_Scalar ))) + NuV * (xz_dz_xr(xr_dz_xz( xz_Scalar ))) 
    
  end function xz_NumDiffKm

  function xz_NumDiffScalar(xz_Scalar)
    !
    ! x, z 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !スカラー量
    real(8)             :: xz_NumDiffScalar(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !水平方向の数値拡散
    
!    xz_NumDiffScalar = 0.0d0  !初期化
    xz_NumDiffScalar = NuH * (xz_dx_pz(pz_dx_xz( xz_Scalar ))) + NuV * (xz_dz_xr(xr_dz_xz( xz_Scalar ))) 
    
    call StorePotTempDiff( xz_NumDiffScalar )

  end function xz_NumDiffScalar

  function xza_NumDiffScalar(xza_Scalar)
    !
    ! x, z 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xza_Scalar(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                    !スカラー量
    real(8)             :: xza_NumDiffScalar(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                    !水平方向の数値拡散
    integer             :: s
    

    do s = 1, SpcNum
      xza_NumDiffScalar(:,:,s) = NuH * (xz_dx_pz(pz_dx_xz( xza_Scalar(:,:,s) ))) + NuV * (xz_dz_xr(xr_dz_xz( xza_Scalar(:,:,s) ))) 
    end do

    call StoreMixRtDiff( xza_NumDiffScalar )

  end function xza_NumDiffScalar