!= 鉛直拡散フラックス (Mellor and Yamada, 1974)
!
!= Vertical diffusion flux (Mellor and Yamada, 1974)
!
! Authors::   Yoshiyuki O. Takahashi, Yasuhiro MORIKAWA, Yukiko YAMADA
! Version::   $Id: vdiffusion_my1974.f90,v 1.29 2011-07-08 13:18:21 yot Exp $
! Tag Name::  $Name: dcpam5-20110714 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module vdiffusion_my1974
  !
  != 鉛直拡散フラックス (Mellor and Yamada, 1974)
  !
  != Vertical diffusion flux (Mellor and Yamada, 1974)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 鉛直拡散フラックスを計算します. 
  !
  ! Vertical diffusion flux is calculated. 
  !
  !== Procedures List
  ! 
  ! VDiffusion       :: 鉛直拡散フラックスの計算
  ! VDiffusionOutPut :: フラックスの出力
  ! ------------     :: ------------
  ! VDiffusion       :: Calculate vertical diffusion fluxes
  ! VDiffusionOutPut :: Output fluxes
  !
  !== NAMELIST
  !
  ! NAMELIST#vdiffusion_my1974_nml
  !

  ! モジュール引用 ; USE statements
  !
  !== References
  !
  ! * Mellor, G. L, Yamada, T., 1974: 
  !   A hierarchy of turbulence closure models for 
  !   planetary boundary layers. 
  !   <i>J. Atmos. Sci.</i>, <b>31</b>, 1791--1806.

  ! 格子点設定
  ! 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

  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !
  use composition, only: ncmax, IndexH2OVap

  ! 種別型パラメタ
  ! 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 :: VDiffusion
  public :: VDiffusionExpTendency
  public :: VDiffusionOutPut
  public :: VDiffusionInit

  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !
  logical, save :: FlagCalcRiWithTv

  logical, save :: vdiffusion_my1974_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag

  ! 鉛直拡散フラックス算出のための変数
  ! Variables for calculation of vertical diffusion flux
  !
  logical , save:: FlagConstDiffCoef
                              ! Flag for use of constant diffusion coefficient
  real(DP), save:: ConstDiffCoefM
                              ! Diffusion coefficient for momentum that is used in the
                              ! case of FlagConstDiffCoef == .true.
  real(DP), save:: ConstDiffCoefH
                              ! Diffusion coefficient for heat that is used in the
                              ! case of FlagConstDiffCoef == .true.
  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:: VelDiffCoefMin
                              ! $ \Dvect{u} $ 拡散係数最小値. 
                              ! Minimum diffusion coefficient of $ \Dvect{u} $
  real(DP), save:: TempDiffCoefMin
                              ! $ T $ 拡散係数最小値. 
                              ! Minimum diffusion coefficient of $ T $
  real(DP), save:: QMixDiffCoefMin
                              ! $ q $ 拡散係数最小値. 
                              ! Minimum diffusion coefficient of $ q $
  real(DP), save:: VelDiffCoefMax
                              ! $ \Dvect{u} $ 拡散係数最大値. 
                              ! Maximum diffusion coefficient of $ \Dvect{u} $
  real(DP), save:: TempDiffCoefMax
                              ! $ T $ 拡散係数最大値. 
                              ! Maximum diffusion coefficient of $ T $
  real(DP), save:: QMixDiffCoefMax
                              ! $ 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: dcpam5-20110714 $' // &
    & '$Id: vdiffusion_my1974.f90,v 1.29 2011-07-08 13:18:21 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

  !--------------------------------------------------------------------------------------

  subroutine VDiffusion(                                         &
    & xyz_U,      xyz_V,      xyzf_QMix,                         & ! (in)
    & xyz_Temp,   xyr_Temp,   xyr_Press,                         & ! (in)
    & xy_SurfHeight,                                             & ! (in)
    & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner,              & ! (in)
    & xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, & ! (out)
    & xyr_VelTransCoef, xyr_TempTransCoef,                       & ! (out)
    & xyr_QMixTransCoef                                          & ! (out)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! 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
      & EpsV
                              ! $ \epsilon_v $ .
                              ! 水蒸気分子量比.
                              ! Molecular weight of water vapor

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & 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):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q $ .   質量混合比. Mass mixing ratio
    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):: xy_SurfHeight (0:imax-1,1:jmax)
                              ! $ z_s $ . 地表面高度. 
                              ! Surface height. 
    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_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward momentum flux
    real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 質量フラックス. 
                              ! Mass flux of compositions
    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_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass (composition)

    ! 作業変数
    ! 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_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity

    real(DP) :: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Virtual temperature
    real(DP) :: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              !
                              ! Virtual temperature

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

    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. vdiffusion_my1974_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


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


    !
    ! Calculate virtual temperature
    ! (If FlagCalcRiWithTv is false, virtual temperature is replaced with temperature.)
    !
    if ( FlagCalcRiWithTv ) then
      xyz_VirTemp = xyz_Temp &
        & * ( 1.0_DP + ( 1.0_DP / EpsV - 1.0_DP ) * xyzf_QMix(:,:,:,IndexH2OVap) )
      xyr_VirTemp(:,:,0) = 1.0d100
      do k = 1, kmax-1
        xyr_VirTemp(:,:,k) = ( xyz_VirTemp(:,:,k) + xyz_VirTemp(:,:,k+1) ) * 0.5_DP
      end do
      xyr_VirTemp(:,:,kmax) = 1.0d100
    else
      xyz_VirTemp = xyz_Temp
      xyr_VirTemp = xyr_Temp
    end if


    ! 拡散係数の計算
    ! Calculation of diffusion coefficient
    !
    if ( FlagConstDiffCoef ) then

      xyr_VelDiffCoef (:,:,0       ) = 0.0d0
      xyr_VelDiffCoef (:,:,1:kmax-1) = ConstDiffCoefM
      xyr_VelDiffCoef (:,:,kmax    ) = 0.0d0

      xyr_TempDiffCoef(:,:,0       ) = 0.0d0
      xyr_TempDiffCoef(:,:,1:kmax-1) = ConstDiffCoefH
      xyr_TempDiffCoef(:,:,kmax    ) = 0.0d0

      xyr_QMixDiffCoef(:,:,0       ) = 0.0d0
      xyr_QMixDiffCoef(:,:,1:kmax-1) = ConstDiffCoefH
      xyr_QMixDiffCoef(:,:,kmax    ) = 0.0d0

    else

      ! バルク $ 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 / ( xyr_VirTemp(:,:,k) / xyr_Exner(:,:,k) )  &
          &  * (   xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1)   &
          &      - xyz_VirTemp(:,:,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 VDiffCoefficient(                                  &
        & xy_SurfHeight,                                      & ! (in)
        & xyr_Height, xyr_DVelDz, xyr_BulkRiNum,              & ! (in)
        & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef & ! (out)
        & )

    end if

    ! 浅い積雲対流
    ! 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_QMixDiffCoef )


    ! 輸送係数の計算
    ! Calculate transfer coefficient
    !
    xyr_VelTransCoef (:,:,0)    = 0.
    xyr_VelTransCoef (:,:,kmax) = 0.
    xyr_TempTransCoef(:,:,0)    = 0.
    xyr_TempTransCoef(:,:,kmax) = 0.
    xyr_QMixTransCoef(:,:,0)    = 0.
    xyr_QMixTransCoef(:,:,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_QMixTransCoef(:,:,k) = &
        &                 xyr_QMixDiffCoef(:,:,k) &
        &                   * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) &
        &                   / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
    end do

    ! フラックスの計算
    ! Calculate fluxes
    !
    xyr_MomFluxX(:,:,0)    = 0.
    xyr_MomFluxX(:,:,kmax) = 0.
    xyr_MomFluxY(:,:,0)    = 0.
    xyr_MomFluxY(:,:,kmax) = 0.
    xyr_HeatFlux(:,:,0)    = 0.
    xyr_HeatFlux(:,:,kmax) = 0.
    do n = 1, ncmax
      xyrf_QMixFlux(:,:,0,n)    = 0.
      xyrf_QMixFlux(:,:,kmax,n) = 0.
    end do

    do k = 1, kmax-1
      xyr_MomFluxX(:,:,k) =   &
        & - xyr_VelTransCoef(:,:,k) * ( xyz_U(:,:,k+1) - xyz_U(:,:,k) )

      xyr_MomFluxY(:,:,k) =   &
        & - xyr_VelTransCoef(:,:,k) * ( xyz_V(:,:,k+1) - xyz_V(:,:,k) )

      xyr_HeatFlux(:,:,k) = &
        & - CpDry * xyr_TempTransCoef(:,:,k) * xyr_Exner(:,:,k) &
        &   * (   xyz_Temp(:,:,k+1) / xyz_Exner(:,:,k+1) &
        &       - xyz_Temp(:,:,k)   / xyz_Exner(:,:,k)     )
    end do

    do n = 1, ncmax
      do k = 1, kmax-1
        xyrf_QMixFlux(:,:,k,n) = &
          & - xyr_QMixTransCoef(:,:,k) &
          &   * ( xyzf_QMix(:,:,k+1,n) - xyzf_QMix(:,:,k,n)  )
      end do
    end do

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

  end subroutine VDiffusion

  !--------------------------------------------------------------------------------------

  subroutine VDiffCoefficient(                            &
    & xy_SurfHeight,                                      & ! (in)
    & xyr_Height, xyr_DVelDz, xyr_BulkRiNum,              & ! (in)
    & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef & ! (out)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! 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):: xy_SurfHeight (0:imax-1,1:jmax)
                              ! $ z_s $ . 地表面高度. 
                              ! Surface height. 
    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_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：質量
                              ! Diffusion coefficient: mass of constituents



    ! 作業変数
    ! 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 check
    !
    if ( .not. vdiffusion_my1974_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! 定数計算
    ! Calculate constants
    !
    Gamma1 = ( 1. / 3. ) - ( 2. * MYLv2ParamA1 / MYLv2ParamB1 )
    Gamma2 =   ( MYLv2ParamB2 / MYLv2ParamB1 ) &
      &      + ( 6. * MYLv2ParamA1 / MYLv2ParamB1 )
    Alpha1 = 3.  * MYLv2ParamA2 * Gamma1
    Alpha2 = 3.  * MYLv2ParamA2 * ( Gamma1 + Gamma2 )
    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
    !
    do k = 0, kmax
      xyr_MixLength(:,:,k) = &
        & FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) &
        & / (1. + FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) / MixLengthMax )
    end do

    ! 拡散係数の算出
    ! 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_QMixDiffCoef = 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_QMixDiffCoef(i,j,k) = &
            & max( min( xyr_QMixDiffCoef(i,j,k), QMixDiffCoefMax ), &
            &      QMixDiffCoefMin &
            &    )
        end do
      end do
    end do

    xyr_VelDiffCoef (:,:,0)    = 0.
    xyr_VelDiffCoef (:,:,kmax) = 0.
    xyr_TempDiffCoef(:,:,0)    = 0.
    xyr_TempDiffCoef(:,:,kmax) = 0.
    xyr_QMixDiffCoef(:,:,0)    = 0.
    xyr_QMixDiffCoef(:,:,kmax) = 0.

  end subroutine VDiffCoefficient

  !--------------------------------------------------------------------------------------

  subroutine VDiffusionExpTendency(                            &
    & xyr_Press,                                               & ! (in )
    & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (in ) optional
    & xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt            & ! (out) optional
    & )
    !
    ! 時間変化率の計算を行います.
    !
    ! Calculate tendencies.
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, &               ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル).
                              ! Air pressure (half level)
    real(DP), intent(in ), optional :: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス.
                              ! Eastward momentum flux
    real(DP), intent(in ), optional :: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス.
                              ! Northward momentum flux
    real(DP), intent(in ), optional :: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス.
                              ! Heat flux
    real(DP), intent(in ), optional :: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 比湿フラックス.
                              ! Specific humidity flux
    real(DP), intent(out), optional :: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{u}{t} $ . 東西風速変化.
                              ! Eastward wind tendency
    real(DP), intent(out), optional :: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{v}{t} $ . 南北風速変化.
                              ! Northward wind tendency
    real(DP), intent(out), optional :: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化.
                              ! Temperature tendency
    real(DP), intent(out), optional :: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \DP{q}{t} $ . 質量混合比変化.
                              ! Mass mixing ratio tendency

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

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. vdiffusion_my1974_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! Check arguments
    !
    if ( present( xyz_DUDt ) ) then
      if ( .not. present( xyr_MomFluxX ) ) then
        call MessageNotify( 'E', module_name, 'xyr_MomFluxX has to be present.' )
      end if
    end if
    if ( present( xyz_DVDt ) ) then
      if ( .not. present( xyr_MomFluxY ) ) then
        call MessageNotify( 'E', module_name, 'xyr_MomFluxY has to be present.' )
      end if
    end if
    if ( present( xyz_DTempDt ) ) then
      if ( .not. present( xyr_HeatFlux ) ) then
        call MessageNotify( 'E', module_name, 'xyr_HeatFlux has to be present.' )
      end if
    end if
    if ( present( xyzf_DQMixDt ) ) then
      if ( .not. present( xyrf_QMixFlux ) ) then
        call MessageNotify( 'E', module_name, 'xyrf_QMixFlux has to be present.' )
      end if
    end if


    if ( present( xyz_DUDt ) ) then
      do k = 1, kmax
        xyz_DUDt(:,:,k) =                                            &
          & + Grav * ( xyr_MomFluxX(:,:,k) - xyr_MomFluxX(:,:,k-1) ) &
          &        / ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) )
      end do
    end if
    if ( present( xyz_DVDt ) ) then
      do k = 1, kmax
        xyz_DVDt(:,:,k) =                                            &
          & + Grav * ( xyr_MomFluxY(:,:,k) - xyr_MomFluxY(:,:,k-1) ) &
          &        / ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) )
      end do
    end if
    if ( present( xyz_DTempDt ) ) then
      do k = 1, kmax
        xyz_DTempDt(:,:,k) =                                         &
          & + Grav / CpDry                                           &
          &        * ( xyr_HeatFlux(:,:,k) - xyr_HeatFlux(:,:,k-1) ) &
          &        / ( xyr_Press   (:,:,k) - xyr_Press   (:,:,k-1) )
      end do
    end if

    if ( present( xyzf_DQMixDt ) ) then
      do n = 1, ncmax
        do k = 1, kmax
          xyzf_DQMixDt(:,:,k,n) =                                            &
            & + Grav * ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) ) &
            &        / ( xyr_Press    (:,:,k)   - xyr_Press    (:,:,k-1)   )
        end do
      end do
    end if


  end subroutine VDiffusionExpTendency

  !--------------------------------------------------------------------------------------

  subroutine VDiffusionOutPut(                                 &
    & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (in)
    & xyz_DUDt,  xyz_DVDt,  xyz_DTempDt,  xyzf_DQMixDt,        & ! (in)
    & xyz_Exner, xyr_Exner,                                    & ! (in)
    & xyr_VelTransCoef, xyr_TempTransCoef,                     & ! (in)
    & xyr_QMixTransCoef                                        & ! (in)
    & )
    !
    ! フラックス (xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux). 
    ! について, その他の引数を用いて補正し, 出力を行う. 
    !
    ! Fluxes (xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux) are
    ! corrected by using other arguments, and the corrected values are output.
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & 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: &
      & DelTime, &            ! $ \Delta t $ [s]
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

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

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward wind flux
    real(DP), intent(in):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 質量フラックス. 
                              ! Mass flux of constituents
    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):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \DP{q}{t} $ . 混合比変化. 
                              ! Mass mixing ratio 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_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass of constituents

    ! 出力のための作業変数
    ! Work variables for output
    !
    real(DP):: xyr_MomFluxXCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward momentum flux
    real(DP):: xyr_MomFluxYCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP):: xyr_HeatFluxCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP):: xyrf_QMixFluxCor(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 質量フラックス. 
                              ! Mass flux of constituents


    ! 作業変数
    ! 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
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents
    real(DP):: LCp
                              ! $ L / C_p $ [K]. 

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. vdiffusion_my1974_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


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


    ! 風速, 温度, 比湿フラックス補正
    ! 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_MomFluxXCor( i,j,k ) = xyr_MomFluxX( i,j,k )  &
            & + ( xyz_DUDt( i,j,k ) - xyz_DUDt( i,j,k+1 ) ) &
            &       * xyr_VelTransCoef( i,j,k ) * DelTime

          xyr_MomFluxYCor( i,j,k ) = xyr_MomFluxY( i,j,k )  &
            & + ( xyz_DVDt( i,j,k ) - xyz_DVDt( i,j,k+1 ) ) &
            &       * xyr_VelTransCoef( i,j,k ) * DelTime

          xyr_HeatFluxCor( i,j,k ) = xyr_HeatFlux( 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
        end do
      end do
    end do

    do n = 1, ncmax
      do k = 1, kmax-1
        do j = 1, jmax
          do i = 0, imax-1
            xyrf_QMixFluxCor( i,j,k,n ) = xyrf_QMixFlux( i,j,k,n ) &
              & + ( xyzf_DQMixDt( i,j,k,n ) - xyzf_DQMixDt( i,j,k+1,n ) ) &
              &       * CpDry * xyr_QMixTransCoef( i,j,k ) * LCp * DelTime
          end do
        end do
      end do
    end do

    xyr_MomFluxXCor   (:,:,0)    = 0.
    xyr_MomFluxXCor   (:,:,kmax) = 0.
    xyr_MomFluxYCor   (:,:,0)    = 0.
    xyr_MomFluxYCor   (:,:,kmax) = 0.
    xyr_HeatFluxCor   (:,:,0)    = 0.
    xyr_HeatFluxCor   (:,:,kmax) = 0.
    do n = 1, ncmax
      xyrf_QMixFluxCor(:,:,0,n)    = 0.
      xyrf_QMixFluxCor(:,:,kmax,n) = 0.
    end do

    ! MEMO
    ! Output values of surface fluxes in MomFluxX, MomFluxY, HeatFlux, and QVapFlux 
    ! are not correct. (YOT, 2009/08/14)
    ! Please refer to output variables, 'TauX', 'TauY', 'Sens', and 'Evap' for those 
    ! values.  (YOT, 2011/05/28)

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'MomFluxX', xyr_MomFluxXCor  )
    call HistoryAutoPut( TimeN, 'MomFluxY', xyr_MomFluxYCor  )
    call HistoryAutoPut( TimeN, 'HeatFlux', xyr_HeatFluxCor  )
    call HistoryAutoPut( TimeN, 'QVapFlux', xyrf_QMixFluxCor )

    call HistoryAutoPut( TimeN, 'DUDtVDiff'   , xyz_DUDt                        )
    call HistoryAutoPut( TimeN, 'DVDtVDiff'   , xyz_DVDt                        )
    call HistoryAutoPut( TimeN, 'DTempDtVDiff', xyz_DTempDt                     )
    call HistoryAutoPut( TimeN, 'DQVapDtVDiff', xyzf_DQMixDt(:,:,:,IndexH2OVap) )

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

  end subroutine VDiffusionOutPut

  !--------------------------------------------------------------------------------------

  subroutine VDiffusionInit
    !
    ! 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/    &
      & FlagConstDiffCoef,              &
      & ConstDiffCoefM, ConstDiffCoefH, &
!
      & SquareVelMin, BulkRiNumMin,     &
!
      & MixLengthMax, TildeShMin, TildeSmMin, &
      & VelDiffCoefMin, TempDiffCoefMin, QMixDiffCoefMin, &
      & VelDiffCoefMax, TempDiffCoefMax, QMixDiffCoefMax, &
!
      & MYLv2ParamA1, MYLv2ParamB1, MYLv2ParamA2, MYLv2ParamB2, MYLv2ParamC1, &
      & FlagCalcRiWithTv
          !
          ! デフォルト値については初期化手続 "vdiffusion_my1974#VDiffInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "vdiffusion_my1974#VDiffInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( vdiffusion_my1974_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    FlagConstDiffCoef = .false.
    ConstDiffCoefM    = 0.0d0
    ConstDiffCoefH    = 0.0d0

    SquareVelMin    =     0.1
    BulkRiNumMin    = - 100.

    MixLengthMax    = 300.
    TildeShMin      =   0.
    TildeSmMin      =   0.
    VelDiffCoefMin  =   0.1
    TempDiffCoefMin =   0.1
    QMixDiffCoefMin =   0.1
    VelDiffCoefMax  = 10000.
    TempDiffCoefMax = 10000.
    QMixDiffCoefMax = 10000.

    ! Parameters proposed by Mellor and Yamada (1982).
    !
    MYLv2ParamA1 =  0.92
    MYLv2ParamB1 = 16.6
    MYLv2ParamA2 =  0.74
    MYLv2ParamB2 = 10.1
    MYLv2ParamC1 =  0.08

    FlagCalcRiWithTv = .false.


    ! 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 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( 'MomFluxX', &
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &
      & 'eastward momentum flux', 'N m-2' )
    call HistoryAutoAddVariable( 'MomFluxY', &
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &
      & 'northward momentum flux', 'N m-2' )
    call HistoryAutoAddVariable( 'HeatFlux', &
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &
      & 'heat flux', 'W m-2' )
    call HistoryAutoAddVariable( 'QVapFlux', &
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &
      & 'moisture flux', 'W m-2' )

    call HistoryAutoAddVariable( 'DUDtVDiff', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'tendency of zonal wind by vertical diffusion', 'm s-2' )
    call HistoryAutoAddVariable( 'DVDtVDiff', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'tendency of meridional wind by vertical diffusion', 'm s-2' )
    call HistoryAutoAddVariable( 'DTempDtVDiff', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'tendency of temperature by vertical diffusion', 'K s-1' )
    call HistoryAutoAddVariable( 'DQVapDtVDiff', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'tendency of specific humidity by vertical diffusion', 's-1' )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'For vertical diffusion flux:' )
    call MessageNotify( 'M', module_name, '  FlagConstDiffCoef = %b', l = (/ FlagConstDiffCoef /) )
    call MessageNotify( 'M', module_name, '  ConstDiffCoefM    = %f', d = (/ ConstDiffCoefM /) )
    call MessageNotify( 'M', module_name, '  ConstDiffCoefH    = %f', d = (/ ConstDiffCoefH /) )
    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, '  QMixDiffCoefMin   = %f', d = (/ QMixDiffCoefMin /) )
    call MessageNotify( 'M', module_name, '  VelDiffCoefMax    = %f', d = (/ VelDiffCoefMax  /) )
    call MessageNotify( 'M', module_name, '  TempDiffCoefMax   = %f', d = (/ TempDiffCoefMax /) )
    call MessageNotify( 'M', module_name, '  QMixDiffCoefMax   = %f', d = (/ QMixDiffCoefMax /) )
    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, '  FlagCalcRiWithTvf = %b', l = (/ FlagCalcRiWithTv /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    vdiffusion_my1974_inited = .true.

  end subroutine VDiffusionInit

  !--------------------------------------------------------------------------------------

end module vdiffusion_my1974
