subroutine VerticalDiffusion( xyz_U,      xyz_V,      xyz_QVap, xyz_Temp,   xyr_Temp,   xyr_Press, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyr_UFlux,  xyr_VFlux,  xyr_TempFlux, xyr_QVapFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QVapTransCoef )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated. 
    !

    ! モジュール引用 ; USE statements
    !

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: FKarm, Grav, GasRDry, CpDry, LatentHeat
                              ! $ L $ [J kg-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimeN, TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u $ .   東西風速. Eastward wind
    real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v $ .   南北風速. Northward wind
    real(DP), intent(in):: xyz_QVap (0:imax-1, 1:jmax, 1:kmax)
                              ! $ q $ .   比湿. Specific humidity
    real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .   温度. Temperature
    real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T} $ . 温度 (半整数レベル). 
                              ! Temperature (half level)
    real(DP), intent(in):: xyr_Press  (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
                              ! 高度 (整数レベル). 
                              ! Height (full level)
    real(DP), intent(in):: xyr_Height (0:imax-1, 1:jmax, 0:kmax)
                              ! 高度 (半整数レベル). 
                              ! Height (half level)
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(out):: xyr_UFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西風速フラックス. 
                              ! Eastward wind flux
    real(DP), intent(out):: xyr_VFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北風速フラックス. 
                              ! Northward wind flux
    real(DP), intent(out):: xyr_TempFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 温度フラックス. 
                              ! Temperature flux
    real(DP), intent(out):: xyr_QVapFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux
    real(DP), intent(out):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP), intent(out):: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(out):: xyr_QVapTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：比湿. 
                              ! Transfer coefficient: specific humidity

    ! 作業変数
    ! Work variables
    !
    real(DP):: xyr_DVelDz (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \DD{|\Dvect{v}|}{z} $
    real(DP):: xyr_BulkRiNum (0:imax-1, 1:jmax, 0:kmax)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $
    real(DP):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP):: xyr_QVapDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity

    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement
    !

    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )

    ! 初期化
    ! Initialization
    !
    if ( .not. vdiffusion_my1974_inited ) call VtclDiffInit

    ! バルク $ R_i $ 数算出
    ! Calculate bulk $ R_i $
    !
    xyr_DVelDz(:,:,0)       = 0.
    xyr_DVelDz(:,:,kmax)    = 0.
    xyr_BulkRiNum(:,:,0)    = 0.
    xyr_BulkRiNum(:,:,kmax) = 0.

    do k = 1, kmax-1
      xyr_DVelDz(:,:,k) = sqrt( max( SquareVelMin , ( xyz_U(:,:,k+1) - xyz_U(:,:,k) )**2 + ( xyz_V(:,:,k+1) - xyz_V(:,:,k) )**2 ) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
      
      xyr_BulkRiNum(:,:,k) = Grav / BasePotTemp * (   xyz_Temp(:,:,k+1) / xyz_Exner(:,:,k+1) - xyz_Temp(:,:,k)   / xyz_Exner(:,:,k)   ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) / xyr_DVelDz(:,:,k)**2
      
      xyr_BulkRiNum(:,:,k) = max( xyr_BulkRiNum(:,:,k) , BulkRiNumMin )
      
    end do
    
    ! 拡散係数の計算
    ! Calculate diffusion coefficients
    !
    call VtclDiffCoefficient( xyr_Height, xyr_DVelDz, xyr_BulkRiNum, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QVapDiffCoef ) ! (out)
    
    ! 浅い積雲対流
    ! Shallow cumulus convection
    !
    ! (AGCM5 から導入予定)
    
    
    ! 拡散係数の出力
    ! Output diffusion coefficients
    !
    ! (上記の「浅い積雲対流」導入後に作成)

    ! 拡散係数出力
    ! Diffusion coeffficients output
    !
    call HistoryAutoPut( TimeN, 'VelDiffCoef',  xyr_VelDiffCoef )
    call HistoryAutoPut( TimeN, 'TempDiffCoef', xyr_TempDiffCoef )
    call HistoryAutoPut( TimeN, 'QVapDiffCoef', xyr_QVapDiffCoef )

    
    ! 輸送係数の計算
    ! Calculate transfer coefficient
    !
    xyr_VelTransCoef (:,:,0)    = 0.
    xyr_VelTransCoef (:,:,kmax) = 0.
    xyr_TempTransCoef(:,:,0)    = 0.
    xyr_TempTransCoef(:,:,kmax) = 0.
    xyr_QVapTransCoef(:,:,0)    = 0.
    xyr_QVapTransCoef(:,:,kmax) = 0.
    
    do k = 1, kmax-1
      xyr_VelTransCoef(:,:,k) = xyr_VelDiffCoef(:,:,k) * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )

      xyr_TempTransCoef(:,:,k) = xyr_TempDiffCoef(:,:,k) * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )

      xyr_QVapTransCoef(:,:,k) = xyr_QVapDiffCoef(:,:,k) * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
    end do
    
    ! フラックスの計算
    ! Calculate fluxes
    !
    xyr_UFlux(:,:,0)       = 0.
    xyr_UFlux(:,:,kmax)    = 0.
    xyr_VFlux(:,:,0)       = 0.
    xyr_VFlux(:,:,kmax)    = 0.
    xyr_TempFlux(:,:,0)    = 0.
    xyr_TempFlux(:,:,kmax) = 0.
    xyr_QVapFlux(:,:,0)    = 0.
    xyr_QVapFlux(:,:,kmax) = 0.

    do k = 1, kmax-1
      xyr_UFlux(:,:,k) = - xyr_VelTransCoef(:,:,k) * ( xyz_U(:,:,k+1) - xyz_U(:,:,k) )
      
      xyr_VFlux(:,:,k) = - xyr_VelTransCoef(:,:,k) * ( xyz_V(:,:,k+1) - xyz_V(:,:,k) )
      
      xyr_TempFlux(:,:,k) = - CpDry * xyr_TempTransCoef(:,:,k) * xyr_Exner(:,:,k) * (   xyz_Temp(:,:,k+1) / xyz_Exner(:,:,k+1) - xyz_Temp(:,:,k)   / xyz_Exner(:,:,k)     )
      
      xyr_QVapFlux(:,:,k) = - LatentHeat * xyr_QVapTransCoef(:,:,k) * ( xyz_QVap(:,:,k+1) - xyz_QVap(:,:,k)  )
      
    end do
    
    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine VerticalDiffusion

  subroutine VerticalDiffusionOutPut( xyr_UFlux, xyr_VFlux, xyr_TempFlux, xyr_QVapFlux, xyz_DUDt,  xyz_DVDt,  xyz_DTempDt,  xyz_DQVapDt, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QVapTransCoef )
    !
    ! フラックス (xyr_UFlux, xyr_VFlux, xyr_TempFlux, xyr_QVapFlux). 
    ! について, その他の引数を用いて補正し, 出力を行う. 
    !
    ! Fluxes (xyr_UFlux, xyr_VFlux, xyr_TempFlux, xyr_QVapFlux) are
    ! corrected by using other arguments, and the corrected values are output.
    !

    ! モジュール引用 ; USE statements
    !

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: GasRDry, CpDry, LatentHeat
                              ! $ L $ [J kg-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xyr_UFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西風速フラックス. 
                              ! Eastward wind flux
    real(DP), intent(in):: xyr_VFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北風速フラックス. 
                              ! Northward wind flux
    real(DP), intent(in):: xyr_TempFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 温度フラックス. 
                              ! Temperature flux
    real(DP), intent(in):: xyr_QVapFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux
    real(DP), intent(in):: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{u}{t} $ . 東西風速変化. 
                              ! Eastward wind tendency
    real(DP), intent(in):: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{v}{t} $ . 南北風速変化. 
                              ! Northward wind tendency
    real(DP), intent(in):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化. 
                              ! Temperature tendency
    real(DP), intent(in):: xyz_DQVapDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{q}{t} $ . 比湿変化. 
                              ! Temperature tendency
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)
    real(DP), intent(in):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP), intent(in):: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xyr_QVapTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：比湿. 
                              ! Transfer coefficient: specific humidity

    ! 出力のための作業変数
    ! Work variables for output
    !
    real(DP):: xyr_UFluxCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西風速フラックス. 
                              ! Eastward wind flux
    real(DP):: xyr_VFluxCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北風速フラックス. 
                              ! Northward wind flux
    real(DP):: xyr_TempFluxCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 温度フラックス. 
                              ! Temperature flux
    real(DP):: xyr_QVapFluxCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux


    ! 作業変数
    ! Work variables
    !
    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    real(DP):: LCp
                              ! $ L / C_p $ [K]. 

    ! 実行文 ; Executable statement
    !

    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )

    ! 初期化
    ! Initialization
    !
    if ( .not. vdiffusion_my1974_inited ) call VtclDiffInit

    ! 風速, 温度, 比湿フラックス補正
    ! Correct fluxes of wind, temperature, specific humidity
    !
    LCp = LatentHeat / CpDry

    do k = 1, kmax-1
      do j = 1, jmax
        do i = 0, imax-1

          xyr_UFluxCor( i,j,k ) = xyr_UFlux( i,j,k ) + ( xyz_DUDt( i,j,k ) - xyz_DUDt( i,j,k+1 ) ) * xyr_VelTransCoef( i,j,k ) * DelTime

          xyr_VFluxCor( i,j,k ) = xyr_VFlux( i,j,k ) + ( xyz_DVDt( i,j,k ) - xyz_DVDt( i,j,k+1 ) ) * xyr_VelTransCoef( i,j,k ) * DelTime

          xyr_TempFluxCor( i,j,k ) = xyr_TempFlux( i,j,k ) + (   xyz_DTempDt( i,j,k   ) / xyz_Exner( i,j,k   ) - xyz_DTempDt( i,j,k+1 ) / xyz_Exner( i,j,k+1 ) ) * CpDry * xyr_TempTransCoef( i,j,k ) * xyr_Exner( i,j,k ) * DelTime

          xyr_QVapFluxCor( i,j,k ) = xyr_QVapFlux( i,j,k ) + ( xyz_DQVapDt( i,j,k ) - xyz_DQVapDt( i,j,k+1 ) ) * CpDry * xyr_QVapTransCoef( i,j,k ) * LCp * DelTime

        end do
      end do
    end do

    xyr_UFluxCor   (:,:,0)    = 0.
    xyr_UFluxCor   (:,:,kmax) = 0.
    xyr_VFluxCor   (:,:,0)    = 0.
    xyr_VFluxCor   (:,:,kmax) = 0.
    xyr_TempFluxCor(:,:,0)    = 0.
    xyr_TempFluxCor(:,:,kmax) = 0.
    xyr_QVapFluxCor(:,:,0)    = 0.
    xyr_QVapFluxCor(:,:,kmax) = 0.

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'UFlux', xyr_UFluxCor )
    call HistoryAutoPut( TimeN, 'VFlux', xyr_VFluxCor )
    call HistoryAutoPut( TimeN, 'TempFlux', xyr_TempFluxCor )
    call HistoryAutoPut( TimeN, 'QVapFlux', xyr_QVapFluxCor )

    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine VerticalDiffusionOutPut

  subroutine InitCheck
    !
    ! 依存モジュールの初期化チェック
    !
    ! Check initialization of dependency modules

    ! モジュール引用 ; USE statements
    !

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_util_inited

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: gridset_inited

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: constants_inited

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: axesset_inited

    ! 時刻管理
    ! Time control
    !
    use timeset, only: timeset_inited

    ! 実行文 ; Executable statement
    !

    if ( .not. namelist_util_inited ) call MessageNotify( 'E', module_name, '"namelist_util" module is not initialized.' )

    if ( .not. gridset_inited ) call MessageNotify( 'E', module_name, '"gridset" module is not initialized.' )

    if ( .not. constants_inited ) call MessageNotify( 'E', module_name, '"constants" module is not initialized.' )

    if ( .not. axesset_inited ) call MessageNotify( 'E', module_name, '"axesset" module is not initialized.' )

    if ( .not. timeset_inited ) call MessageNotify( 'E', module_name, '"timeset" module is not initialized.' )

  end subroutine InitCheck

  subroutine VtclDiffCoefficient( xyr_Height, xyr_DVelDz, xyr_BulkRiNum, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QVapDiffCoef )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated. 
    !

    ! モジュール引用 ; USE statements
    !

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: FKarm
                              ! $ k $ .
                              ! カルマン定数. 
                              ! Karman constant

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xyr_Height (0:imax-1, 1:jmax, 0:kmax)
                              ! 高度 (半整数レベル). 
                              ! Height (half level)
    real(DP), intent(in):: xyr_DVelDz (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \DD{|\Dvect{v}|}{z} $
    real(DP), intent(in):: xyr_BulkRiNum (0:imax-1, 1:jmax, 0:kmax)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $
    real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(out):: xyr_QVapDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity



    ! 作業変数
    ! Work variables
    !
    real(DP):: xyr_FluxRiNum (0:imax-1, 1:jmax, 0:kmax)
                              ! フラックス $ R_i $ 数. 
                              ! Flux $ R_i $ number
    real(DP):: xyr_TildeSh (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \tilde{S_h} $ (温度, 比湿). 
                              ! $ \tilde{S_h} $ (temperature, specific humidity)
    real(DP):: xyr_TildeSm (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \tilde{S_m} $ (運動量). 
                              ! $ \tilde{S_m} $ (momentum)
    real(DP):: xyr_MixLength (0:imax-1, 1:jmax, 0:kmax)
                              ! 混合距離. 
                              ! Mixing length

    real(DP):: Alpha1, Alpha2
    real(DP):: Beta1, Beta2, Beta3, Beta4
    real(DP):: Gamma1, Gamma2
    real(DP):: CrtlFluxRiNum

    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. vdiffusion_my1974_inited ) call VtclDiffInit

    ! 定数計算
    ! Calculate constants
    !
    Gamma1 = ( 1. / 3. ) - ( 2. * MYLv2ParamA1 / MYLv2ParamB1 )
    Gamma2 =   ( MYLv2ParamB2 / MYLv2ParamB1 ) + ( 6. * MYLv2ParamA1 / MYLv2ParamB1 )
    Alpha1 = 3.  * MYLv2ParamA2 * Gamma1
    Alpha2 = 3.  * MYLv2ParamA2 * ( Gamma1 + Gamma1 )
    Beta1  = MYLv2ParamA1 * MYLv2ParamB1 * ( Gamma1 - MYLv2ParamC1 )
    Beta2  = MYLv2ParamA1 * (   MYLv2ParamB1 * ( Gamma1 - MYLv2ParamC1 ) + 6. * MYLv2ParamA1 + 3. * MYLv2ParamA2 )
    Beta3  = MYLv2ParamA2 * MYLv2ParamB1 * Gamma1
    Beta4  = MYLv2ParamA2 * (   MYLv2ParamB1 * ( Gamma1 + Gamma2 ) - 3. * MYLv2ParamA1 )
    CrtlFluxRiNum = Gamma1 / ( Gamma1 + Gamma2 )

    ! フラックス $ R_i $ 数の算出
    ! Calculate flux $ R_i $ number
    !
    xyr_FluxRiNum = (   Beta1 + Beta4 * xyr_BulkRiNum - sqrt(   ( Beta1 + Beta4 * xyr_BulkRiNum )**2 - 4. * Beta2 * Beta3 * xyr_BulkRiNum ) ) / ( 2. * Beta2 )
 
    ! $ \tilde{S_h} $ と $ \tilde{S_m} $ の算出
    ! Calculate $ \tilde{S_h} $ and $ \tilde{S_m} $
    !
    xyr_TildeSh(:,:,kmax) = 0.
    xyr_TildeSm(:,:,kmax) = 0.
    
    do k = 0, kmax-1
      do i = 0, imax-1
        do j = 1, jmax
          
          if ( xyr_FluxRiNum(i,j,k) < CrtlFluxRiNum ) then 
            
            xyr_TildeSh(i,j,k) = (   Alpha1 - Alpha2 * xyr_FluxRiNum(i,j,k) ) / ( 1. - 1. * xyr_FluxRiNum(i,j,k) )

            xyr_TildeSm(i,j,k) = (   Beta1 - Beta2 * xyr_FluxRiNum(i,j,k) ) / ( Beta3 - Beta4 * xyr_FluxRiNum(i,j,k) ) * xyr_TildeSh(i,j,k)

            xyr_TildeSh(i,j,k) = max( xyr_TildeSh(i,j,k), TildeShMin )
            xyr_TildeSm(i,j,k) = max( xyr_TildeSm(i,j,k), TildeSmMin )
            
          else
            
            xyr_TildeSh(i,j,k) = TildeShMin
            xyr_TildeSm(i,j,k) = TildeSmMin
            
          end if
          
        end do
      end do
    end do
    
    
    ! 混合距離の算出
    ! Calculate mixing length
    !
    xyr_MixLength = FKarm * xyr_Height / (1. + FKarm * xyr_Height / MixLengthMax )
    
    ! 拡散係数の算出
    ! Calculate diffusion constants
    !
    xyr_VelDiffCoef = xyr_MixLength**2 * xyr_DVelDz * sqrt ( MYLv2ParamB1 * ( 1. - xyr_FluxRiNum ) * xyr_TildeSm ) * xyr_TildeSm
    
    xyr_TempDiffCoef = xyr_MixLength ** 2 * xyr_DVelDz * sqrt ( MYLv2ParamB1 * ( 1. - xyr_FluxRiNum ) * xyr_TildeSm ) * xyr_TildeSh
    
    xyr_QVapDiffCoef = xyr_TempDiffCoef
    
    do k = 0, kmax-1
      do i = 0, imax-1
        do j = 1, jmax
          xyr_VelDiffCoef(i,j,k) = max( min( xyr_VelDiffCoef(i,j,k), VelDiffCoefMax ), VelDiffCoefMin )
          xyr_TempDiffCoef(i,j,k) = max( min( xyr_TempDiffCoef(i,j,k), TempDiffCoefMax ), TempDiffCoefMin )
          xyr_QVapDiffCoef(i,j,k) = max( min( xyr_QVapDiffCoef(i,j,k), QVapDiffCoefMax ), QVapDiffCoefMin )
        end do
      end do
    end do
    
    xyr_VelDiffCoef (:,:,0)    = 0.
    xyr_VelDiffCoef (:,:,kmax) = 0.
    xyr_TempDiffCoef(:,:,0)    = 0.
    xyr_TempDiffCoef(:,:,kmax) = 0.
    xyr_QVapDiffCoef(:,:,0)    = 0.
    xyr_QVapDiffCoef(:,:,kmax) = 0.

  end subroutine VtclDiffCoefficient

  subroutine VtclDiffInit
    !
    ! vdiffusion_my1974 モジュールの初期化を行います. 
    ! NAMELIST#vdiffusion_my1974_nml の読み込みはこの手続きで行われます. 
    !
    ! "vdiffusion_my1974" module is initialized. 
    ! "NAMELIST#vdiffusion_my1974_nml" is loaded in this procedure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: StoA

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 宣言文 ; Declaration statements
    !
    implicit none

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /vdiffusion_my1974_nml/ BasePotTemp, SquareVelMin, BulkRiNumMin, MixLengthMax, TildeShMin, TildeSmMin, VelDiffCoefMin, TempDiffCoefMin, QVapDiffCoefMin, VelDiffCoefMax, TempDiffCoefMax, QVapDiffCoefMax, MYLv2ParamA1, MYLv2ParamB1, MYLv2ParamA2, MYLv2ParamB2, MYLv2ParamC1
          !
          ! デフォルト値については初期化手続 "vdiffusion_my1974#VtclDiffInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "vdiffusion_my1974#VtclDiffInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( vdiffusion_my1974_inited ) return
    call InitCheck

    ! デフォルト値の設定
    ! Default values settings
    !
    BasePotTemp     =   300.
    SquareVelMin    =     0.1
    BulkRiNumMin    = - 100.

    MixLengthMax    = 300.
    TildeShMin      =   0.
    TildeSmMin      =   0.
    VelDiffCoefMin  =   0.1
    TempDiffCoefMin =   0.1
    QVapDiffCoefMin =   0.1
    VelDiffCoefMax  = 10000.
    TempDiffCoefMax = 10000.
    QVapDiffCoefMax = 10000.

    MYLv2ParamA1 =  0.92
    MYLv2ParamB1 = 16.6
    MYLv2ParamA2 =  0.74
    MYLv2ParamB2 = 10.1
    MYLv2ParamC1 =  0.08


    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = vdiffusion_my1974_nml, iostat = iostat_nml )           ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
      if ( iostat_nml == 0 ) write( STDOUT, nml = vdiffusion_my1974_nml )
    end if


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'VelDiffCoef', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'diffusion coef. momentum', 'm2 s-1' )
    call HistoryAutoAddVariable( 'TempDiffCoef', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'diffusion coef. heat    ', 'm2 s-1' )
    call HistoryAutoAddVariable( 'QVapDiffCoef', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'diffusion coef. moisture', 'm2 s-1' )

    call HistoryAutoAddVariable( 'UFlux', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'eastward momentum flux', 'N m-2' )
    call HistoryAutoAddVariable( 'VFlux', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'northward momentum flux', 'N m-2' )
    call HistoryAutoAddVariable( 'TempFlux', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'heat flux', 'W m-2' )
    call HistoryAutoAddVariable( 'QVapFlux', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'moisture flux', 'W m-2' )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'For vertical diffusion flux:' )
    call MessageNotify( 'M', module_name, '  BasePotTemp  = %f', d = (/ BasePotTemp /) )
    call MessageNotify( 'M', module_name, '  SquareVelMin = %f', d = (/ SquareVelMin /) )
    call MessageNotify( 'M', module_name, '  BulkRiNumMin = %f', d = (/ BulkRiNumMin /) )
    call MessageNotify( 'M', module_name, 'For diffusion coefficients:' )
    call MessageNotify( 'M', module_name, '  MixLengthMax     = %f', d = (/ MixLengthMax     /) )
    call MessageNotify( 'M', module_name, '  TildeShMin       = %f', d = (/ TildeShMin       /) )
    call MessageNotify( 'M', module_name, '  TildeSmMin       = %f', d = (/ TildeSmMin       /) )
    call MessageNotify( 'M', module_name, '  VelDiffCoefMin   = %f', d = (/ VelDiffCoefMin  /) )
    call MessageNotify( 'M', module_name, '  TempDiffCoefMin  = %f', d = (/ TempDiffCoefMin /) )
    call MessageNotify( 'M', module_name, '  QVapDiffCoefMin  = %f', d = (/ QVapDiffCoefMin /) )
    call MessageNotify( 'M', module_name, '  VelDiffCoefMax   = %f', d = (/ VelDiffCoefMax  /) )
    call MessageNotify( 'M', module_name, '  TempDiffCoefMax  = %f', d = (/ TempDiffCoefMax /) )
    call MessageNotify( 'M', module_name, '  QVapDiffCoefMax  = %f', d = (/ QVapDiffCoefMax /) )
    call MessageNotify( 'M', module_name, '  MYLv2ParamA1     = %f', d = (/ MYLv2ParamA1     /) )
    call MessageNotify( 'M', module_name, '  MYLv2ParamB1     = %f', d = (/ MYLv2ParamB1     /) )
    call MessageNotify( 'M', module_name, '  MYLv2ParamA2     = %f', d = (/ MYLv2ParamA2     /) )
    call MessageNotify( 'M', module_name, '  MYLv2ParamB2     = %f', d = (/ MYLv2ParamB2     /) )
    call MessageNotify( 'M', module_name, '  MYLv2ParamC1     = %f', d = (/ MYLv2ParamC1     /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    vdiffusion_my1974_inited = .true.
  end subroutine VtclDiffInit