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

    !暗黙の型宣言禁止
    implicit none

    ! CReSS マニュアルでは, 2 次精度中心差分の場合, Alpha < 1/8 くらいが適当と述べている.
    NuH  = Alpha * ( DelX ** 2.0d0 ) / DelTimeLong
    NuV  = Alpha * ( DelZ ** 2.0d0 ) / DelTimeLong

    !確認
    if (cpurank == 0) then
      call MessageNotify( "M", "NumDiffusion_init", "Alpha = %f", d=(/Alpha/) )
      call MessageNotify( "M", "NumDiffusion_init", "NuH = %f", d=(/NuH/) )
      call MessageNotify( "M", "NumDiffusion_init", "NuV = %f", d=(/NuV/) )
      if (Alpha > 0.125) then
        call MessageNotify( "E", "NumDiffusion_init", "Alpha = %f > 1/8", d=(/Alpha/) )
        stop
      end if
    end if
  end subroutine NumDiffusion_init

  subroutine NumDiffusion_init2(cfgfile)
    !
    ! NumDiffusion モジュールの初期化ルーチン. Alpha はデフォルト値を使う場合. 
    !
    
    !暗黙の型宣言禁止
    implicit none

    !
    !NAMELIST から放射強制の設定を取得
    !
    character(*), intent(in) :: cfgfile

    ! NAMELIST から情報を取得
    NAMELIST /numdiffusion/ Alpha_Velocity

    open (10, FILE=cfgfile)
    read(10, NML=numdiffusion)
    close(10)    

    ! 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/) )
    call MessageNotify( "M", "NumDiffusion_init", "Alpha_Velocity = %f", d=(/Alpha_Velocity/) )

  end subroutine NumDiffusion_init2

  subroutine NumDiffusion_init3(cfgfile)
    !
    ! NumDiffusion モジュールの初期化ルーチン
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !入力変数
    character(*), intent(in) :: cfgfile
    real(8)                  :: Coeff

    ! Namelist から Alpha の値を決め直す.
    NAMELIST /numdiffusion/ Coeff
    open (10, FILE=cfgfile)
    read(10, NML=numdiffusion)
    close(10)
    
    Alpha = Coeff
    if (cpurank == 0) then 
       call MessageNotify( "M", "NumDiffusion_init", "Coeff = %f", d=(/Coeff/) )
    end if
    
    ! CReSS マニュアルでは, 2 次精度中心差分の場合, Alpha < 1/8 くらいが適当と述べている. 
    NuH  = Alpha * ( DelX ** 2.0d0 ) / DelTimeLong
    NuV  = Alpha * ( DelZ ** 2.0d0 ) / DelTimeLong

    ! Alpha_Velocity は 1.0 にしておく. 整合性のため. 
    Alpha_Velocity = 1.0d0
    
    !確認
    if (cpurank == 0) then 
      call MessageNotify( "M", "NumDiffusion_init", "Alpha = %f", d=(/Alpha/) )
      call MessageNotify( "M", "NumDiffusion_init", "NuH = %f", d=(/NuH/) )
      call MessageNotify( "M", "NumDiffusion_init", "NuV = %f", d=(/NuV/) )
      if (Alpha > 0.125) then 
        call MessageNotify( "E", "NumDiffusion_init", "Alpha = %f > 1/8", d=(/Alpha/) )
        stop
      end if
    end if
  end subroutine NumDiffusion_init3

  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)
                                                    !数値拡散
    
    !速度の拡散係数だけ変更する場合にAlpha_Velocityの値を1以外の値にする
!   pz_NumDiffVelX = 0.0d0  !初期化
    pz_NumDiffVelX = Alpha_Velocity * NuH * ( pz_dx_xz( xz_dx_pz( pz_VarX ) ) ) + Alpha_Velocity * NuV * ( pz_dz_pr( pr_dz_pz( pz_VarX ) ) )
    
    call StoreMomDiff( pz_NumDiffVelX )

  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)
                                                    !数値拡散

    !速度の拡散係数だけ変更する場合にAlpha_Velocityの値を1以外の値にする
!    xr_NumDiffVelZ = 0.0d0 
    xr_NumDiffVelZ = Alpha_Velocity * NuH * ( xr_dx_pr( pr_dx_xr( xr_VarZ ) ) ) + Alpha_Velocity * 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