!= 鉛直拡散フラックス
!
!= Vertical diffusion flux
!
! Authors::   Yasuhiro MORIKAWA, Yukiko YAMADA
! Version::   $Id: $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module vdiffusion_my1974
  !
  != 鉛直拡散フラックス
  !
  != Vertical diffusion flux
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 鉛直拡散フラックスを計算します. 
  !
  ! Vertical diffusion flux is calculated. 
  !
  !== Procedures List
  ! 
  ! VerticalDiffusion :: 鉛直拡散フラックスの計算
  ! ------------      :: ------------
  ! VerticalDiffusion :: Calculate vertical diffusion flux
  !
  !== NAMELIST
  !
  ! NAMELIST#vdiffusion_my1974_nml
  !

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

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, & ! 経度格子点数. 
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数. 
                             ! Number of grid points in latitude
    &                kmax    ! 鉛直層数. 
                             ! Number of vertical level

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING     ! 文字列.       Strings. 

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public:: VerticalDiffusion

  ! 公開変数
  ! Public variables
  !
  logical, save, public:: vdiffusion_my1974_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag

  ! 非公開変数
  ! Private variables
  !


  ! 鉛直拡散フラックス算出のための変数
  ! Variables for calculation of vertical diffusion flux
  !
  real(DP), save:: RefPress
                              ! 参照気圧. 
                              ! Reference air pressure
  real(DP), save:: BasePotTemp
                              ! 基本温位. 
                              ! Base potential temperature
  real(DP), save:: SquareVelMin
                              ! 風二乗差最小値. 
                              ! Minimum value of square of velocity
  real(DP), save:: BulkRiNumMin
                              ! バルク $ R_i $ 数最小値. 
                              ! Minimum value of bulk $ R_i $

  ! 拡散係数算出のための変数
  ! Veriables for calculation of diffusion coefficients
  !
  real(DP), save:: MixLengthMax
                              ! 最大混合距離. 
                              ! Maximum mixing length
  real(DP), save:: TildeShMin
                              ! $ \tilde{S_h} $ 最小値. 
                              ! Minimum $ \tilde{S_h} $
  real(DP), save:: TildeSmMin
                              ! $ \tilde{S_m} $ 最小値. 
                              ! Minimum $ \tilde{S_m} $
  real(DP), save:: VelDiffCoeffMin
                              ! $ \Dvect{u} $ 拡散係数最小値. 
                              ! Minimum diffusion coefficient of $ \Dvect{u} $
  real(DP), save:: TempDiffCoeffMin
                              ! $ T $ 拡散係数最小値. 
                              ! Minimum diffusion coefficient of $ T $
  real(DP), save:: QvapDiffCoeffMin
                              ! $ q $ 拡散係数最小値. 
                              ! Minimum diffusion coefficient of $ q $
  real(DP), save:: VelDiffCoeffMax
                              ! $ \Dvect{u} $ 拡散係数最大値. 
                              ! Maximum diffusion coefficient of $ \Dvect{u} $
  real(DP), save:: TempDiffCoeffMax
                              ! $ T $ 拡散係数最大値. 
                              ! Maximum diffusion coefficient of $ T $
  real(DP), save:: QvapDiffCoeffMax
                              ! $ q $ 拡散係数最大値. 
                              ! Maximum diffusion coefficient of $ q $

  ! Mellor Yamada Level 2 定数
  ! Constants for Mellor Yamada Level 2
  !
  real(DP), save:: MYLv2ParamA1
  real(DP), save:: MYLv2ParamB1
  real(DP), save:: MYLv2ParamA2
  real(DP), save:: MYLv2ParamB2
  real(DP), save:: MYLv2ParamC1


  character(*), parameter:: module_name = 'vdiffusion_my1974'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id:  $'
                              ! モジュールのバージョン
                              ! Module version

  ! INTERFACE 文 ; INTERFACE statements
  !
  interface VerticalDiffusion
    module procedure VerticalDiffusion
  end interface

contains

  subroutine VerticalDiffusion( &
    & xyz_U,     xyz_V,     xyz_QVap, &                   ! (in)
    & xyz_Temp,  xyr_Temp, &                              ! (in)
    & xyz_Press, xyr_Press, xyz_GeoPot,   xyr_GeoPot, &   ! (in)
    & xyr_UFlux, xyr_VFlux, xyr_TempFlux, xyr_QVapFlux, & ! (inout)
    & xyza_UVMtx, xyra_TempMtx, xyza_QVapMtx &            ! (inout)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated. 
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & FKarm, &
                              ! $ k $ .
                              ! カルマン定数. 
                              ! Karman constant
      & Grav, & 
                              ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration
      & GasRDry, &
                              ! $ R $ [J kg-1 K-1]. 
                              ! 乾燥大気の気体定数. 
                              ! Gas constant of air
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure
      & LatentHeat
                              ! $ L $ [J kg-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation

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

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

    ! 宣言文 ; 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):: xyz_Press  (0:imax-1, 1:jmax, 1:kmax)
                              ! $ p $ . 気圧 (整数レベル). 
                              ! Air pressure (full 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_GeoPot (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \phi $ . ジオポテンシャル (整数レベル). 
                              ! Geo-potential (full level)
    real(DP), intent(in):: xyr_GeoPot (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{\phi} $ . ジオポテンシャル (半整数レベル). 
                              ! Geo-potential (half level)

    real(DP), intent(inout):: xyr_UFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西風速フラックス. 
                              ! Zonal wind flux
    real(DP), intent(inout):: xyr_VFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北風速フラックス. 
                              ! Meridional wind flux
    real(DP), intent(inout):: xyr_TempFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 温度フラックス. 
                              ! Temperature flux
    real(DP), intent(inout):: xyr_QVapFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux
    real(DP), intent(inout):: xyza_UVMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 速度陰解行列. 
                              ! Implicit matrix about velocity 
    real(DP), intent(inout):: xyra_TempMtx (0:imax-1, 1:jmax, 0:kmax, -1:1)
                              ! 温度陰解行列. 
                              ! Implicit matrix about temperature
    real(DP), intent(inout):: xyza_QVapMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 比湿陰解行列. 
                              ! Implicit matrix about 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_VelTransCoeff (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP):: xyr_TempTransCoeff (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP):: xyr_QvapTransCoeff (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：比湿. 
                              ! Transfer coefficient: specific humidity
    real(DP):: xyr_VelDiffCoeff (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP):: xyr_TempDiffCoeff (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP):: xyr_QvapDiffCoeff (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity
    real(DP):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    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 VtcllDiffInit

    ! Exner 関数算出
    ! Calculate Exner functions
    !
    xyz_Exner = ( xyz_Press / RefPress ) ** ( GasRDry / CpDry )
    xyr_Exner = ( xyr_Press / RefPress ) ** ( GasRDry / CpDry )

    ! バルク $ R_i $ 数算出
    ! Calculate bulk $ R_i $
    !
    xyr_DVelDz    = 0.0_DP
    xyr_BulkRiNum = 0.0_DP
    
    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_GeoPot(:,:,k+1) - xyz_GeoPot(:,:,k) )
      
      xyr_BulkRiNum(:,:,k) = &
        &  Grav / BasePotTemp &
        &  * (   xyz_Temp(:,:,k+1) / xyz_Exner(:,:,k+1)   &
        &      - xyz_Temp(:,:,k)   / xyz_Exner(:,:,k)   ) &
        &  / ( xyz_GeoPot(:,:,k+1) - xyz_GeoPot(:,:,k) ) &
        &  / xyr_DVelDz(:,:,k)**2
      
      xyr_BulkRiNum(:,:,k) = max( xyr_BulkRiNum(:,:,k) , BulkRiNumMin )
      
    end do
    
    ! 拡散係数の計算
    ! Calculate diffusion coefficients
    !
    call VtclDiffCoefficient( &
      & xyr_GeoPot, xyr_DVelDz, xyr_BulkRiNum, &                 ! (in)
      & xyr_VelDiffCoeff, xyr_TempDiffCoeff, xyr_QvapDiffCoeff ) ! (out)
    
    ! 浅い積雲対流
    ! Shallow cumulus convection
    !
    ! (AGCM5 から導入予定)
    
    
    ! 拡散係数の出力
    ! Output diffusion coefficients
    !
    ! (上記の「浅い積雲対流」導入後に作成)

    
    ! 輸送係数の計算
    ! Calculate transfer coefficient
    !
    xyr_VelTransCoeff  = 0.0_DP
    xyr_TempTransCoeff = 0.0_DP
    xyr_QvapTransCoeff = 0.0_DP
    
    do k = 1, kmax-1
      xyr_VelTransCoeff(:,:,k) = &
        &                 xyr_VelDiffCoeff(:,:,k) &
        &                   * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) &
        &                   / ( xyz_GeoPot(:,:,k+1) - xyz_GeoPot(:,:,k) )

      xyr_TempTransCoeff(:,:,k) = &
        &                 xyr_TempDiffCoeff(:,:,k) &
        &                   * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) &
        &                   / ( xyz_GeoPot(:,:,k+1) - xyz_GeoPot(:,:,k) )

      xyr_QvapTransCoeff(:,:,k) = &
        &                 xyr_QvapDiffCoeff(:,:,k) &
        &                   * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) &
        &                   / ( xyz_GeoPot(:,:,k+1) - xyz_GeoPot(:,:,k) )
    end do
    
    ! フラックスの計算
    ! Calculate fluxes
    !
    do k = 1, kmax-1
      xyr_UFlux(:,:,k) =   &
        &   xyr_UFlux(:,:,k) &
        & + xyr_VelTransCoeff(:,:,k) * ( xyz_U(:,:,k) - xyz_U(:,:,k+1) )
      
      xyr_VFlux(:,:,k) =   &
        &   xyr_VFlux(:,:,k) &
        & + xyr_VelTransCoeff(:,:,k) * ( xyz_V(:,:,k) - xyz_V(:,:,k+1) )
      
      xyr_TempFlux(:,:,k) = &
        &   xyr_TempFlux(:,:,k) &
        & + CpDry * xyr_TempTransCoeff(:,:,k) * xyr_Exner(:,:,k) &
        &   * (   xyz_Temp(:,:,k)   / xyz_Exner(:,:,k)      &
        &       - xyz_Temp(:,:,k+1) / xyz_Exner(:,:,k+1) )
      
      xyr_QvapFlux(:,:,k) = &
        &   xyr_QvapFlux(:,:,k) &
        & + LatentHeat * xyr_QvapTransCoeff(:,:,k) &
        &   * ( xyz_Qvap(:,:,k) - xyz_Qvap(:,:,k+1) )
      
    end do
    

    ! 陰解行列の計算
    ! Calculate implicit matrices
    !
    do k = 2, kmax
      xyza_UVMtx(:,:,k, 0)   = xyza_UVMtx(:,:,k,0) + xyr_VelTransCoeff(:,:,k-1)

      xyza_UVMtx(:,:,k,-1)   =                     - xyr_VelTransCoeff(:,:,k-1)
      
      xyra_TempMtx(:,:,k, 0) =   xyra_TempMtx(:,:,k,0) &
        &                      + CpDry * xyr_TempTransCoeff(:,:,k-1) &
        &                          * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k)

      xyra_TempMtx(:,:,k,-1) = - CpDry * xyr_TempTransCoeff(:,:,k-1) &
        &                          * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k-1)
      
      xyza_QvapMtx(:,:,k, 0) =   xyza_QvapMtx(:,:,k,0) &
        &                      + CpDry * xyr_QvapTransCoeff(:,:,k-1)

      xyza_QvapMtx(:,:,k,-1) = - CpDry * xyr_QvapTransCoeff(:,:,k-1)
    end do
    
    do k = 1, kmax-1
      xyza_UVMtx(:,:,k,0)    =   xyza_UVMtx(:,:,k,0) &
        &                      + xyr_VelTransCoeff(:,:,k)

      xyza_UVMtx(:,:,k,1)    = - xyr_VelTransCoeff(:,:,k)
      
      xyra_TempMtx(:,:,k,0)  =   xyra_TempMtx(:,:,k,0) &
        &                      + CpDry * xyr_TempTransCoeff(:,:,k) &
        &                          * xyr_Exner(:,:,k) / xyz_Exner(:,:,k)

      xyra_TempMtx(:,:,k,1)  = - CpDry * xyr_TempTransCoeff(:,:,k) &
        &                          * xyr_Exner(:,:,k) / xyz_Exner(:,:,k+1)
      
      xyza_QvapMtx(:,:,k,0)  =   xyza_QvapMtx(:,:,k,0) &
        &                      + CpDry * xyr_QvapTransCoeff(:,:,k)

      xyza_QvapMtx(:,:,k,1)  = - CpDry * xyr_QvapTransCoeff(:,:,k)
    end do

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryFileOutput( 'UDiffFlux',    xyr_UFlux )
    call HistoryFileOutput( 'VDiffFlux',    xyr_VFlux )
    call HistoryFileOutput( 'TempDiffFlux', xyr_TempFlux )
    call HistoryFileOutput( 'QVapDiffFlux', xyr_QVapFlux )

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

  end subroutine VerticalDiffusion



  subroutine VtclDiffCoefficient( &
    & xyr_GeoPot, xyr_DVelDz, xyr_BulkRiNum, &                 ! (in)
    & xyr_VelDiffCoeff, xyr_TempDiffCoeff, xyr_QvapDiffCoeff & ! (out)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated. 
    !

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

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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xyr_GeoPot (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{\phi} $ . ジオポテンシャル (半整数レベル). 
                              ! Geo-potential (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_VelDiffCoeff (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: xyr_TempDiffCoeff (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(out):: xyr_QvapDiffCoeff (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 VtcllDiffInit

    ! 定数計算
    ! Calculate constants
    !
    Gamma1 = ( 1.0_DP / 3.0_DP ) - ( 2.0_DP * MYLv2ParamA1 / MYLv2ParamB1 )
    Gamma2 =   ( MYLv2ParamB2 / MYLv2ParamB1 ) &
      &      + ( 6.0_DP * MYLv2ParamA1 / MYLv2ParamB1 )
    Alpha1 = 3.0_DP  * MYLv2ParamA2 * Gamma1
    Alpha2 = 3.0_DP  * MYLv2ParamA2 * ( Gamma1 + Gamma1 )
    Beta1  = MYLv2ParamA1 * MYLv2ParamB1 * ( Gamma1 - MYLv2ParamC1 )
    Beta2  = MYLv2ParamA1 &
      &      * (   MYLv2ParamB1 * ( Gamma1 - MYLv2ParamC1 ) &
      &          + 6.0_DP * MYLv2ParamA1 &
      &          + 3.0_DP * MYLv2ParamA2 )
    Beta3  = MYLv2ParamA2 * MYLv2ParamB1 * Gamma1
    Beta4  = MYLv2ParamA2 * (   MYLv2ParamB1 * ( Gamma1 + Gamma2 ) &
      &                       - 3.0_DP * 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.0_DP * Beta2 * Beta3 * xyr_BulkRiNum &
      &         ) &
      & ) / ( 2.0_DP * Beta2 )
    
    ! $ \tilde{S_h} $ と $ \tilde{S_m} $ の算出
    ! Calculate $ \tilde{S_h} $ and $ \tilde{S_m} $
    !
    xyr_TildeSh = 0.0_DP
    xyr_TildeSm = 0.0_DP
    
    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.0_DP - 1.0_DP * 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_GeoPot / (1.0_DP + FKarm * xyr_GeoPot / MixLengthMax )
    
    ! 拡散係数の算出
    ! Calculate diffusion constants
    !
    xyr_VelDiffCoeff = &
      &   xyr_MixLength**2 &
      &   * xyr_DVelDz &
      &   * sqrt ( MYLv2ParamB1 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_TildeSm ) &
      &   * xyr_TildeSm
    
    xyr_TempDiffCoeff = &
      &   xyr_MixLength ** 2 &
      &   * xyr_DVelDz &
      &   * sqrt ( MYLv2ParamB1 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_TildeSm ) &
      &   * xyr_TildeSh
    
    xyr_QvapDiffCoeff = xyr_TempDiffCoeff
    
    do k = 0, kmax-1
      do i = 0, imax-1
        do j = 1, jmax
          xyr_VelDiffCoeff(i,j,k) = &
            & max( min( xyr_VelDiffCoeff(i,j,k), VelDiffCoeffMax ), &
            &      VelDiffCoeffMin &
            &    )
          xyr_TempDiffCoeff(i,j,k) = &
            & max( min( xyr_TempDiffCoeff(i,j,k), TempDiffCoeffMax ), &
            &      TempDiffCoeffMin &
            &    )
          xyr_QvapDiffCoeff(i,j,k) = &
            & max( min( xyr_QvapDiffCoeff(i,j,k), QvapDiffCoeffMax ), &
            &      QvapDiffCoeffMin &
            &    )
        end do
      end do
    end do
    
    xyr_VelDiffCoeff(:,:,0)     = 0.0_DP
    xyr_TempDiffCoeff(:,:,0)    = 0.0_DP
    xyr_QvapDiffCoeff(:,:,0)    = 0.0_DP
    xyr_VelDiffCoeff(:,:,kmax)  = 0.0_DP
    xyr_TempDiffCoeff(:,:,kmax) = 0.0_DP
    xyr_QvapDiffCoeff(:,:,kmax) = 0.0_DP

  end subroutine VtclDiffCoefficient



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

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

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

    ! 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 history_file_io, only: HistoryFileRegister

    ! 宣言文 ; 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/ &
      & RefPress, BasePotTemp, SquareVelMin, BulkRiNumMin, &
!
      & MixLengthMax, TildeShMin, TildeSmMin, &
      & VelDiffCoeffMin, TempDiffCoeffMin, QvapDiffCoeffMin, &
      & VelDiffCoeffMax, TempDiffCoeffMax, QvapDiffCoeffMax, &
!
      & MYLv2ParamA1, MYLv2ParamB1, MYLv2ParamA2, MYLv2ParamB2, MYLv2ParamC1
          !
          ! デフォルト値については初期化手続 "vdiffusion_my1974#VtcllDiffInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "vdiffusion_my1974#VtcllDiffInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( vdiffusion_my1974_inited ) return
    call InitCheck

    ! デフォルト値の設定
    ! Default values settings
    !
    RefPress        =   1.0e+6_DP
    BasePotTemp     =   300.0_DP
    SquareVelMin    =     0.1_DP
    BulkRiNumMin    = - 100.0_DP

    MixLengthMax     = 300.0_DP
    TildeShMin       =   0.0_DP
    TildeSmMin       =   0.0_DP
    VelDiffCoeffMin  =   0.1_DP
    TempDiffCoeffMin =   0.1_DP
    QvapDiffCoeffMin =   0.1_DP
    VelDiffCoeffMax  = 10000.0_DP
    TempDiffCoeffMax = 10000.0_DP
    QvapDiffCoeffMax = 10000.0_DP

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


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

      rewind( unit_nml )
      read( unit_nml, &                   ! (in)
        & nml = vdiffusion_my1974_nml, &  ! (out)
        & 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 HistoryFileRegister( 'UDiffFlux',    'eastward wind flux by vertical diffusion', &
      & 'N s-2', StoA('lon', 'lat', 'sigm', 'time') )
    call HistoryFileRegister( 'VDiffFlux',    'northward wind flux by vertical diffusion', &
      & 'N s-2', StoA('lon', 'lat', 'sigm', 'time') )
    call HistoryFileRegister( 'TempDiffFlux', 'temperature flux by vertical diffusion', &
      & 'W s-2', StoA('lon', 'lat', 'sigm', 'time') )
    call HistoryFileRegister( 'QVapDiffFlux', 'specific humidity flux by vertical diffusion', &
      & 'W s-2', StoA('lon', 'lat', 'sigm', 'time') )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'For vertical diffusion flux:' )
    call MessageNotify( 'M', module_name, '  RefPress     = %f', d = (/ RefPress /) )
    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, '  VelDiffCoeffMin  = %f', d = (/ VelDiffCoeffMin  /) )
    call MessageNotify( 'M', module_name, '  TempDiffCoeffMin = %f', d = (/ TempDiffCoeffMin /) )
    call MessageNotify( 'M', module_name, '  QvapDiffCoeffMin = %f', d = (/ QvapDiffCoeffMin /) )
    call MessageNotify( 'M', module_name, '  VelDiffCoeffMax  = %f', d = (/ VelDiffCoeffMax  /) )
    call MessageNotify( 'M', module_name, '  TempDiffCoeffMax = %f', d = (/ TempDiffCoeffMax /) )
    call MessageNotify( 'M', module_name, '  QvapDiffCoeffMax = %f', d = (/ QvapDiffCoeffMax /) )
    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 VtcllDiffInit



  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

end module vdiffusion_my1974