!= 鉛直拡散フラックス (Mellor and Yamada, 1974, 1982)
!
!= Vertical diffusion flux (Mellor and Yamada, 1974, 1982)
!
! Authors::   Yasuhiro MORIKAWA, Yukiko YAMADA, Yoshiyuki O. Takahashi, Shin-ichi Takehiro
! Version::   $Id: vdiffusion_my.f90,v 1.5 2025/09/17 22:00:00 takepiro Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2025. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module vdiffusion_my
  !
  != 鉛直拡散フラックス (Mellor and Yamada, 1974, 1982)
  !
  != Vertical diffusion flux (Mellor and Yamada, 1974, 1982)
  !
  ! <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_my_nml
  !

  ! モジュール引用 ; USE statements
  !
  !== References
  !
  ! * Mellor, G. L., and T. Yamada, 1974: 
  !   A hierarchy of turbulence closure models for 
  !   planetary boundary layers. 
  !   <i>J. Atmos. Sci.</i>, <b>31</b>, 1791--1806.
  ! * Mellor, G. L., and T. Yamada, 1982: 
  !   Development of a turbulence closure model for geophysical fluid problems,
  !   <i>Rev. Geophys. Space Phys.</i>, <b>20</b>, 851--875.
  ! * Galperin, B., L. H. Kantha, S. Hassid, and T. Rosati, 1988:
  !   A quasi-equilibrium turbulent energy model for geophysical flows,
  !   <i>J. Atmos. Sci.</i>, <b>45</b>, 55--62.

  ! 格子点設定
  ! 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 :: VDiffusionMY25
  public :: VDiffusionMY251DWrapper3D
  public :: VDiffusionMY25GBT94
  public :: VDiffusion
  public :: VDiffusionExpTendency
  public :: VDiffusionOutPut
  public :: VDiffusionInit

  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !
  logical, save :: vdiffusion_my_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:: ShMin
                              ! $ S_h $ 最小値. 
                              ! Minimum $ S_h $
  real(DP), save:: SmMin
                              ! $ S_m $ 最小値. 
                              ! Minimum $ 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:: VelDiffCoefMax
                              ! $ \Dvect{u} $ 拡散係数最大値. 
                              ! Maximum diffusion coefficient of $ \Dvect{u} $
  real(DP), save:: TempDiffCoefMax
                              ! $ T $ 拡散係数最大値. 
                              ! Maximum diffusion coefficient of $ T $

  ! Mellor Yamada Level 2 定数
  ! Constants for Mellor Yamada Level 2
  !
  real(DP), save:: MYConstA1
  real(DP), save:: MYConstB1
  real(DP), save:: MYConstA2
  real(DP), save:: MYConstB2
  real(DP), save:: MYConstC1


  character(*), parameter:: module_name = 'vdiffusion_my'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: vdiffusion_my.f90,v 1.5 2025/09/17 22:00:00 takepiro Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

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

  subroutine VDiffusionMY25(                                     &
    & xyz_U, xyz_V, xyzf_QMix,                                   & ! (in)
    & xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press,   & ! (in)
    & xy_SurfHeight,                                             & ! (in)
    & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner,              & ! (in)
    & xyz_TurKinEne,                                             & ! (in)
    & xy_SurfMomFluxX, xy_SurfMomFluxY,                          & ! (in)
    & xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, & ! (out)
    & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef,       & ! (in)
    & xyz_DTurKinEneDt                                           & ! (out)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated by use of MY2.5 model.
    !

    ! モジュール引用 ; 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

    ! 時刻管理
    ! 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

    ! MPI 関連ルーチン
    ! MPI related routines
    !
    use mpi_wrapper, only : MPIWrapperChkTrue

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

    ! 宣言文 ; 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):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T_v $ .   仮温度. Virtual temperature
    real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual 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(in):: xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Turbulent kinetic energy (m2 s-2)
    real(DP), intent(in):: xy_SurfMomFluxX (0:imax-1, 1:jmax)
                              ! 
                              ! Eastward momentum flux at surface
    real(DP), intent(in):: xy_SurfMomFluxY (0:imax-1, 1:jmax)
                              ! 
                              ! Northward momentum flux at surface

    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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Diffusion coefficient: temperature
    real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity

    real(DP), intent(out):: xyz_DTurKinEneDt (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Tendency of turbulent kinetic energy

    ! 作業変数
    ! Work variables
    !

    real(DP) :: xyz_MixLength(0:imax-1, 1:jmax, 1:kmax)
                              ! 混合距離. 
                              ! Mixing length
    real(DP) :: xyz_DVelDzSq(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Vertical shear squared (s-2)
    real(DP) :: xyz_StatStab(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Static stability (s-2)
    real(DP) :: GhMin
                              !
                              ! Minimum of G_h
    real(DP) :: GhMax
                              !
                              ! Maximum of G_h
    real(DP) :: xyz_Gm(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! G_m
    real(DP) :: xyz_Gh(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! G_h
!!$    real(DP) :: xyz_SmHat(0:imax-1, 1:jmax, 1:kmax)
!!$                              !
!!$                              ! \hat{S}_M
!!$    real(DP) :: xyz_ShHat(0:imax-1, 1:jmax, 1:kmax)
!!$                              !
!!$                              ! \hat{S}_h
!!$    real(DP) :: xyz_DSmHatDTKE(0:imax-1, 1:jmax, 1:kmax)
!!$                              !
!!$                              ! derivative of \hat{S}_M
!!$    real(DP) :: xyz_DShHatDTKE(0:imax-1, 1:jmax, 1:kmax)
!!$                              !
!!$                              ! derivative of \hat{S}_h
    real(DP) :: xyz_Sm(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! S_M
    real(DP) :: xyz_Sh(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! S_h

    real(DP), parameter :: Stke = 0.2_DP
                              !
                              ! S_{TKE} = 0.2

    real(DP) :: xyz_VelDiffCoef (0:imax-1, 1:jmax, 1:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP) :: xyz_TempDiffCoef (0:imax-1, 1:jmax, 1:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature

    real(DP) :: xyr_TurKinEneDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Diffusion coefficient: turbulent kinetic energy
    real(DP) :: xyz_TurKinEneDiffCoef (0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Diffusion coefficient: turbulent kinetic energy
    real(DP) :: xyr_TurKinEneTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Transfer coefficient: turbulent kinetic energy

    real(DP) :: xyr_TurKinEneFlux(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Turbulent energy flux

    real(DP) :: xyz_CShe1(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CShe2(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CBuo1(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CBuo2(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CDis1(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CDis2(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_TurKinEneProShear(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_TurKinEneProBuoya(0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xy_FricVelSq    (0:imax-1, 1:jmax)
    real(DP) :: xy_TurKinEneAtLB(0:imax-1, 1:jmax)

    real(DP) :: xyza_TurKinEneMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 
                              ! Implicit matrix for turbulent kinetic energy
    real(DP) :: xyz_TurKinEneVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Implicit vector for turbulent kinetic energy

    real(DP) :: xyza_TurKinEneLUMtx  (0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP) :: xyz_DelTurKinEneLUVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Tendency of turbulent kinetic energy

    real(DP) :: xyz_TurKinEneDiss(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Dissipation rate of turbulent kinetic energy (m2 s-3)

    real(DP) :: xyz_TurKinEneNonZero(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Turbulent kinetic energy with offset (m2 s-2)

    real(DP), parameter :: TurKinEneOffset  = ( 1.0e-3_DP )**2 / 2.0_DP

    logical :: FlagReCalc
                              !
                              ! Flag for recalculation
    logical :: a_FlagReCalcLocal (1)
    logical :: a_FlagReCalcGlobal(1)
    integer :: iloop
    integer :: nloop

    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

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. vdiffusion_my_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 turbulent kinetic energy with offset
    !
    xyz_TurKinEneNonZero = xyz_TurKinEne + TurKinEneOffset

    !
    ! Calculation of vertical shear squared
    do k = 1, kmax
      if ( k == 1 ) then
        xyz_DVelDzSq(:,:,k) =                              &
          &   (   ( xyz_U(:,:,k+1) - xyz_U(:,:,k  ) )**2   &
          &     + ( xyz_V(:,:,k+1) - xyz_V(:,:,k  ) )**2 ) &
          & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k  ) )**2
!!$        xyz_DVelDzSq(:,:,k) =                              &
!!$          &   (   ( xyz_U(:,:,k+1) - 0.0_DP )**2   &
!!$          &     + ( xyz_V(:,:,k+1) - 0.0_DP )**2 ) &
!!$          & / ( xyz_Height(:,:,k+1) - xy_SurfHeight )**2
      else if ( k == kmax ) then
        xyz_DVelDzSq(:,:,k) =                               &
          &    (   ( xyz_U(:,:,k  ) - xyz_U(:,:,k-1) )**2   &
          &      + ( xyz_V(:,:,k  ) - xyz_V(:,:,k-1) )**2 ) &
          &  / ( xyz_Height(:,:,k  ) - xyz_Height(:,:,k-1) )**2
      else
        xyz_DVelDzSq(:,:,k) =                              &
          &   (   ( xyz_U(:,:,k+1) - xyz_U(:,:,k-1) )**2   &
          &     + ( xyz_V(:,:,k+1) - xyz_V(:,:,k-1) )**2 ) &
          & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) )**2
      end if
    end do
    ! Calculation of static stability
    do k = 1, kmax
      if ( k == 1 ) then
        xyz_StatStab(:,:,k) = &
          &  Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) )  &
          &  * (   xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1)   &
          &      - xyz_VirTemp(:,:,k  ) / xyz_Exner(:,:,k  ) ) &
          &  / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k  ) )
      else if ( k == kmax ) then
        xyz_StatStab(:,:,k) = &
          &  Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) )  &
          &  * (   xyz_VirTemp(:,:,k  ) / xyz_Exner(:,:,k  )   &
          &      - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) &
          &  / ( xyz_Height(:,:,k  ) - xyz_Height(:,:,k-1) )
      else
        xyz_StatStab(:,:,k) = &
          &  Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) )  &
          &  * (   xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1)   &
          &      - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) &
          &  / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) )
      end if
    end do

    ! 混合距離の算出
    ! Calculate mixing length
    !
    do k = 1, kmax
      xyz_MixLength(:,:,k) =                            &
        & FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) &
        & / (1.0_DP + FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) / MixLengthMax )
    end do
    !   Limit mixing length (Galperin et al., 1988) and avoid zero
    xyz_MixLength = &
      &   min( xyz_MixLength, &
      &        0.53_DP &
      &          * sqrt( 2.0_DP * xyz_TurKinEneNonZero / max( xyz_StatStab, 1.0e-10_DP ) ) ) &
      & + 1.0e-10_DP

    xyz_Gh = - xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) &
      & * xyz_StatStab
    ! Actually, xyz_Gm is not used below.
    xyz_Gm =   xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) &
      & * xyz_DVelDzSq


    ! Limit Gh (Galperin et al., 1988)
    GhMin = - 0.53_DP**2
    GhMax = 1.0_DP                                         &
      & / ( MYConstA2                                      &
      &     * (  12.0_DP * MYConstA1 + MYConstB1 + 3.0_DP * MYConstB2 ) )
    xyz_Gh = max( GhMin, min( xyz_Gh, GhMax ) )


    xyz_Sh =                                                        &
      &     MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) &
      & / (   1.0_DP                                                &
      &     - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * xyz_Gh )

    xyz_Sm =                                                        &
      & (                                                           &
      &     MYConstA1 * (   1.0_DP - 3.0_DP * MYConstC1             &
      &                   - 6.0_DP * MYConstA1 / MYConstB1 )        &
      &   + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) &
      &     * xyz_Gh * xyz_Sh                                       &
      & )                                                           &
      & / (   1.0_DP                                                &
      &     - 9.0_DP * MYConstA1 * MYConstA2 * xyz_Gh )


!!$    xyz_DShHatDTKE =                                                       &
!!$      & - 2.0_DP * MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) &
!!$      & / (   2.0_DP * xyz_TurKinEneNonZero                                &
!!$      &     - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * xyz_GhPrime )**2
!!$    xyz_DSmHatDTKE =                                                &
!!$      & (                                                           &
!!$      &     9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) &
!!$      &     * xyz_GhPrime * xyz_DShHatDTKE                          &
!!$      &     * (   2.0_DP * xyz_TurKinEneNonZero                     &
!!$      &         - 9.0_DP * MYConstA1 * MYConstA2 * xyz_GhPrime )    &
!!$      &   - 2.0_DP                                                        &
!!$      &     * (                                                           &
!!$      &           MYConstA1 * (   1.0_DP - 3.0_DP * MYConstC1             &
!!$      &                         - 6.0_DP * MYConstA1 / MYConstB1 )        &
!!$      &         + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) &
!!$      &           * xyz_GhPrime * xyz_ShHat                               &
!!$      &       )                                                           &
!!$      & )                                                                 &
!!$      & / (   2.0_DP * xyz_TurKinEneNonZero                               &
!!$      &     - 9.0_DP * MYConstA1 * MYConstA2 * xyz_GhPrime )**2


    ! 拡散係数の計算
    ! Calculation of diffusion coefficient
    !
    xyz_VelDiffCoef  = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sm
    xyz_TempDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sh
    !
    do k = 0, kmax
      if ( ( k == 0 ) .or. ( k == kmax ) ) then
        xyr_VelDiffCoef (:,:,k) = 0.0_DP
        xyr_TempDiffCoef(:,:,k) = 0.0_DP
      else
        xyr_VelDiffCoef (:,:,k) = ( xyz_VelDiffCoef (:,:,k) + xyz_VelDiffCoef (:,:,k+1) ) / 2.0_DP
        xyr_TempDiffCoef(:,:,k) = ( xyz_TempDiffCoef(:,:,k) + xyz_TempDiffCoef(:,:,k+1) ) / 2.0_DP
      end if
    end do
    !
    do k = 1, kmax-1
      do j = 1, jmax
        do i = 0, imax-1
          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                                  &
            &    )
        end do
      end do
    end do
    !
    xyr_QMixDiffCoef      = xyr_TempDiffCoef


    ! 輸送係数とフラックスの計算
    ! Calculate transfer coefficient and flux
    !
    call VDiffusionCalcFlux(                                       &
      & xyz_U, xyz_V, xyzf_QMix,                                   & ! (in)
      & xyz_Temp, xyr_VirTemp, xyr_Press,                          & ! (in)
      & xyz_Height, xyz_Exner, xyr_Exner,                          & ! (in)
      & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef,       & ! (in)
      & xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux  & ! (out)
      & )


    ! Calculate tendency of turbulent kinetic energy

    !   Set diffusion coefficient for turbulent kinetic energy
    xyz_TurKinEneDiffCoef = &
      & xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * Stke
    !
    do k = 0, kmax
      if ( k == 0 ) then
        xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,1)
      else if ( k == kmax ) then
        xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,kmax)
      else
        xyr_TurKinEneDiffCoef(:,:,k) = &
          & ( xyz_TurKinEneDiffCoef(:,:,k) + xyz_TurKinEneDiffCoef(:,:,k+1) ) &
          & / 2.0_DP
      end if
    end do

    !   Calculate turbulent kinetic energy at lower boundary
    !
    xy_FricVelSq = sqrt( xy_SurfMomFluxX**2 + xy_SurfMomFluxY**2 ) &
      & / ( xyr_Press(:,:,0) / ( GasRDry * xyr_VirTemp(:,:,0) ) )
    xy_TurKinEneAtLB = MYConstB1**(2.0_DP/3.0_DP) / 2.0_DP * xy_FricVelSq
    xy_TurKinEneAtLB = xy_TurKinEneAtLB + TurKinEneOffset

    !   Calculate transfer coefficient and flux of turbulent kinetic energy
    !
    !    When transfer coefficient at lower boundary is calculated, 
    !    diffusion coefficient at mid-point of 1st layer is used. 
    !    In addition, transfer coefficient at upper boundary is assumed 
    !    to be zero.
    k = 0
    xyr_TurKinEneTransCoef(:,:,k) =                               &
      &   xyz_TurKinEneDiffCoef(:,:,1)                            &
      &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
      &     / ( xyz_Height(:,:,1) - xy_SurfHeight )
    do k = 1, kmax-1
      xyr_TurKinEneTransCoef(:,:,k) =                               &
        &   xyr_TurKinEneDiffCoef(:,:,k)                            &
        &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
        &     / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
    end do
    k = kmax
    xyr_TurKinEneTransCoef(:,:,k) = 0.0_DP
    !
    do k = 1, kmax-1
      xyr_TurKinEneFlux(:,:,k) =               &
        & - xyr_TurKinEneTransCoef(:,:,k)      &
        &   * ( xyz_TurKinEneNonZero(:,:,k+1) - xyz_TurKinEneNonZero(:,:,k) )
    end do
    k = 0
    xyr_TurKinEneFlux(:,:,k) =               &
      & - xyr_TurKinEneTransCoef(:,:,k)      &
      &   * ( xyz_TurKinEneNonZero(:,:,k+1) - xy_TurKinEneAtLB )
    k = kmax
    xyr_TurKinEneFlux(:,:,k) = 0.0_DP


!!$    xyz_CShe1 =          2.0d0**1.5 * xyz_MixLength * xyz_SmHat &
!!$      & * xyz_DVelDzSq                                          &
!!$!      & * xyz_TurKinEne**1.5
!!$      & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$    xyz_CShe2 = 1.5_DP * 2.0d0**1.5 * xyz_MixLength * xyz_SmHat &
!!$      & * xyz_DVelDzSq                                          &
!!$!      & * xyz_TurKinEne**0.5
!!$      & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5
!!$    xyz_CBuo1 = -          2.0d0**1.5 * xyz_MixLength * xyz_ShHat &
!!$      & * xyz_StatStab                                            &
!!$!      & * xyz_TurKinEne**1.5
!!$      & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$    xyz_CBuo2 = - 1.5_DP * 2.0d0**1.5 * xyz_MixLength * xyz_ShHat &
!!$      & * xyz_StatStab                                            &
!!$!      & * xyz_TurKinEne**0.5
!!$      & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5
!!$    xyz_CDis1 =          2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$!      & * xyz_TurKinEne**1.5
!!$      & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$    xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$!      & * xyz_TurKinEne**0.5
!!$      & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5


!!$!    xyz_CShe1 =          2.0d0**1.5 * xyz_MixLength * xyz_SmHat &
!!$!      & * xyz_DVelDzSq                                          &
!!$!!      & * xyz_TurKinEne**1.5
!!$!      & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$    xyz_CShe1 =          2.0d0**1.5 * xyz_MixLength * xyz_Sm &
!!$      & * xyz_DVelDzSq                                       &
!!$!      & * xyz_TurKinEne**1.5
!!$      & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5
!!$    xyz_CShe2 = 0.0_DP
!!$!    xyz_CBuo1 = -          2.0d0**1.5 * xyz_MixLength * xyz_ShHat &
!!$!      & * xyz_StatStab                                            &
!!$!!      & * xyz_TurKinEne**1.5
!!$!      & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$    xyz_CBuo1 = -          2.0d0**1.5 * xyz_MixLength * xyz_Sh &
!!$      & * xyz_StatStab                                         &
!!$!      & * xyz_TurKinEne**1.5
!!$      & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5
!!$    xyz_CBuo2 = 0.0_DP
!!$    xyz_CDis1 =          2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$      & * xyz_TurKinEne**1.5
!!$!      & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5
!!$    xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$      & * xyz_TurKinEne**0.5
!!$!      & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5


!!$    xyz_CShe1 = 2.0d0**1.5 * xyz_MixLength * xyz_DVelDzSq       &
!!$      & * xyz_SmHat                                             &
!!$      & * xyz_TurKinEneNonZero**1.5
!!$    xyz_CShe2 = 2.0d0**1.5 * xyz_MixLength * xyz_DVelDzSq       &
!!$      & * (   xyz_DSmHatDTKE                                    &
!!$      &         * xyz_TurKinEneNonZero**1.5                     &
!!$      &     + 1.5_DP * xyz_SmHat                                &
!!$      &         * xyz_TurKinEneNonZero**0.5                     &
!!$      &   )
!!$    xyz_CBuo1 = - 2.0d0**1.5 * xyz_MixLength * xyz_StatStab     &
!!$      & * xyz_ShHat                                             &
!!$      & * xyz_TurKinEneNonZero**1.5
!!$    xyz_CBuo2 = - 2.0d0**1.5 * xyz_MixLength * xyz_StatStab     &
!!$      & * (   xyz_DShHatDTKE                                    &
!!$      &         * xyz_TurKinEneNonZero**1.5                     &
!!$      &     + 1.5_DP * xyz_ShHat                                &
!!$      &         * xyz_TurKinEneNonZero**0.5                     &
!!$      &   )
!!$    xyz_CDis1 =          2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$      & * xyz_TurKinEneNonZero**1.5
!!$    xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) &
!!$      & * xyz_TurKinEneNonZero**0.5


    xyz_CShe1 =            sqrt( 2.0_DP ) * xyz_MixLength * xyz_DVelDzSq  &
      & * xyz_Sm                                                          &
      & * sqrt( xyz_TurKinEneNonZero )
!!$    xyz_CShe2 =   0.5_DP * sqrt( 2.0_DP ) * xyz_MixLength * xyz_DVelDzSq  &
!!$      & * xyz_Sm                                                          &
!!$      & / sqrt( xyz_TurKinEneNonZero )
    xyz_CShe2 = 0.0_DP
    xyz_CBuo1 = -          sqrt( 2.0_DP ) * xyz_MixLength * xyz_StatStab  &
      & * xyz_Sh                                                          &
      & * sqrt( xyz_TurKinEneNonZero )
!!$    xyz_CBuo2 = - 0.5_DP * sqrt( 2.0_DP ) * xyz_MixLength * xyz_StatStab  &
!!$      & * xyz_Sh                                                          &
!!$      & / sqrt( xyz_TurKinEneNonZero )
    xyz_CBuo2 = 0.0_DP
    xyz_CDis1 =          2.0_DP**1.5_DP / ( MYConstB1 * xyz_MixLength ) &
      & * xyz_TurKinEneNonZero**1.5_DP
    xyz_CDis2 = 1.5_DP * 2.0_DP**1.5_DP / ( MYConstB1 * xyz_MixLength ) &
      & * xyz_TurKinEneNonZero**0.5_DP


    nloop = kmax
    loop_solve : do iloop = 1, nloop

      !
      ! Construct implicit matrix from transfer coefficient of vertical 
      ! diffusion scheme (turbulent kinetic energy)
      !
      k = 1
      xyza_TurKinEneMtx(:,:,k,-1) = 0.0_DP
      xyza_TurKinEneMtx(:,:,k, 0) =                                             &
        & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
        & + xyr_TurKinEneTransCoef(:,:,k-1)                                     &
        & + xyr_TurKinEneTransCoef(:,:,k  )                                     &
        & + ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav                    &
        &     * ( xyz_CShe2(:,:,k) + xyz_CBuo2(:,:,k) - xyz_CDis2(:,:,k) )
      xyza_TurKinEneMtx(:,:,k, 1) =                                             &
        & - xyr_TurKinEneTransCoef(:,:,k  )
      !
      do k = 2, kmax-1
        xyza_TurKinEneMtx(:,:,k,-1) =&
          & - xyr_TurKinEneTransCoef(:,:,k-1)
        xyza_TurKinEneMtx(:,:,k, 0) =                                             &
          & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
          & + xyr_TurKinEneTransCoef(:,:,k-1)                                     &
          & + xyr_TurKinEneTransCoef(:,:,k  )                                     &
          & + ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav                    &
          &     * ( xyz_CShe2(:,:,k) + xyz_CBuo2(:,:,k) - xyz_CDis2(:,:,k) )
        xyza_TurKinEneMtx(:,:,k, 1) =                                             &
          & - xyr_TurKinEneTransCoef(:,:,k  )
      end do
      !
      k = kmax
      xyza_TurKinEneMtx(:,:,k,-1) =&
        & - xyr_TurKinEneTransCoef(:,:,k-1)
      xyza_TurKinEneMtx(:,:,k, 0) =                                             &
        & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
        & + xyr_TurKinEneTransCoef(:,:,k-1)                                     &
        & + xyr_TurKinEneTransCoef(:,:,k  )                                     &
        & + ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav                    &
        &     * ( xyz_CShe2(:,:,k) + xyz_CBuo2(:,:,k) - xyz_CDis2(:,:,k) )
      xyza_TurKinEneMtx(:,:,k, 1) = 0.0_DP

      do k = 1, kmax
        xyz_TurKinEneVec(:,:,k) =                                               &
          & - ( xyr_TurKinEneFlux(:,:,k) - xyr_TurKinEneFlux(:,:,k-1) )         &
          & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav                  &
          &     * ( xyz_CShe1(:,:,k) + xyz_CBuo1(:,:,k) - xyz_CDis1(:,:,k) )
      end do


      !
      ! Solve simultaneous linear equations by use of LU decomposition technique
      !
      xyza_TurKinEneLUMtx = xyza_TurKinEneMtx
      !
      call PhyImplLUDecomp3(   &
        & xyza_TurKinEneLUMtx, &   ! (inout)
        & imax * jmax, kmax    &   ! (in)
        & )

      xyz_DelTurKinEneLUVec = xyz_TurKinEneVec
      !
      call PhyImplLUSolve3(      &
        & xyz_DelTurKinEneLUVec, & ! (inout)
        & xyza_TurKinEneLUMtx,   & ! (in)
        & 1, imax * jmax , kmax  & ! (in)
        & )

      xyz_DTurKinEneDt = xyz_DelTurKinEneLUVec / ( 2.0_DP * DelTime )


      ! Calculation of dissipation rate of turbulent kinetic energy
      !
      ! Calculate production rate of turbulent kinetic energy
      ! by shear and buoyancy
      xyz_TurKinEneProShear = xyz_CShe1 + xyz_CShe2 * xyz_DTurKinEneDt * 2.0_DP * DelTime
      xyz_TurKinEneProBuoya = xyz_CBuo1 + xyz_CBuo2 * xyz_DTurKinEneDt * 2.0_DP * DelTime
      xyz_TurKinEneDiss     = xyz_CDis1 + xyz_CDis2 * xyz_DTurKinEneDt * 2.0_DP * DelTime

      ! Check of turbulent kinetic energy dissipation rate
      ! If it is negative, tendency is recalculated without dissipation.
      !
      FlagReCalc = .false.
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xyz_TurKinEneDiss(i,j,k) < 0.0_DP ) then
              xyz_CDis1(i,j,k) = 0.0_DP
              xyz_CDis2(i,j,k) = 0.0_DP
              FlagReCalc = .true.
            end if
          end do
        end do
      end do

      ! Check convergence
      a_FlagReCalcLocal = FlagReCalc
      call MPIWrapperChkTrue(   &
        & 1, a_FlagReCalcLocal, & ! (in)
        & a_FlagReCalcGlobal    & ! (out)
        & )
      if ( .not. a_FlagReCalcGlobal(1) ) exit loop_solve

    end do loop_solve


!!$    write( 6, * ) TimeN, iloop
!!$    write( 6, * ) xyz_TurKinEne(0,jmax/2+1,1:4)
!!$    write( 6, * ) xyz_TempDiffCoef(0,jmax/2+1,1:4)
!!$    write( 6, * ) xyz_TurKinEneProShear(0,jmax/2+1,1:4)
!!$    write( 6, * ) xyz_TurKinEneProBuoya(0,jmax/2+1,1:4)
!!$    write( 6, * ) xyz_TurKinEneDiss(0,jmax/2+1,1:4)
!!$    write( 6, * ) xyz_DTurKinEneDt(0,jmax/2+1,1:4)



    ! 拡散係数の出力
    ! Output diffusion coefficients
    !

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

    call HistoryAutoPut( TimeN, 'TKEPShear', xyz_TurKinEneProShear )
    call HistoryAutoPut( TimeN, 'TKEPBuoy' , xyz_TurKinEneProBuoya )
    call HistoryAutoPut( TimeN, 'TKEDiss'  , xyz_TurKinEneDiss     )

    call HistoryAutoPut( TimeN, 'MixLength' , xyz_MixLength )



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

  end subroutine VDiffusionMY25

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

  subroutine VDiffusionMY251DWrapper3D(                          &
    & xyz_U, xyz_V, xyzf_QMix,                                   & ! (in)
    & xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press,   & ! (in)
    & xy_SurfHeight,                                             & ! (in)
    & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner,              & ! (in)
    & xyz_TurKinEne,                                             & ! (in)
    & xy_SurfMomFluxX, xy_SurfMomFluxY,                          & ! (in)
    & xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, & ! (out)
    & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef,       & ! (in)
    & xyz_DTurKinEneDt                                           & ! (out)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated by use of MY2.5 model.
    !

    ! モジュール引用 ; 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

    ! 時刻管理
    ! 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

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

    ! 宣言文 ; 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):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T_v $ .   仮温度. Virtual temperature
    real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual 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(in):: xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Turbulent kinetic energy (m2 s-2)
    real(DP), intent(in):: xy_SurfMomFluxX (0:imax-1, 1:jmax)
                              ! 
                              ! Eastward momentum flux at surface
    real(DP), intent(in):: xy_SurfMomFluxY (0:imax-1, 1:jmax)
                              ! 
                              ! Northward momentum flux at surface

    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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Diffusion coefficient: temperature
    real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity
    real(DP), intent(out):: xyz_DTurKinEneDt (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Tendency of turbulent kinetic energy

    ! 作業変数
    ! Work variables
    !
    real(DP) :: z_U (1:kmax)
    real(DP) :: z_V (1:kmax)
    real(DP) :: zf_QMix(1:kmax, 1:ncmax)
    real(DP) :: z_Temp (1:kmax)
    real(DP) :: r_Temp (0:kmax)
    real(DP) :: z_VirTemp (1:kmax)
    real(DP) :: r_VirTemp (0:kmax)
    real(DP) :: r_Press  (0:kmax)
    real(DP) :: SurfHeight
    real(DP) :: z_Height (1:kmax)
    real(DP) :: r_Height (0:kmax)
    real(DP) :: z_Exner (1:kmax)
    real(DP) :: r_Exner (0:kmax)
    real(DP) :: z_TurKinEne(1:kmax)
    real(DP) :: SurfMomFluxX
    real(DP) :: SurfMomFluxY
    real(DP) :: r_MomFluxX (0:kmax)
    real(DP) :: r_MomFluxY (0:kmax)
    real(DP) :: r_HeatFlux (0:kmax)
    real(DP) :: rf_QMixFlux(0:kmax, 1:ncmax)
    real(DP) :: r_VelDiffCoef (0:kmax)
    real(DP) :: r_TempDiffCoef(0:kmax)
    real(DP) :: r_QMixDiffCoef(0:kmax)
    real(DP) :: z_DTurKinEneDt (1:kmax)
                              !
                              ! Tendency of turbulent kinetic energy


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


    ! 実行文 ; Executable statement
    !

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


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


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

        z_U          = xyz_U          (i,j,:)
        z_V          = xyz_V          (i,j,:)
        zf_QMix      = xyzf_QMix      (i,j,:,:)
        z_Temp       = xyz_Temp       (i,j,:)
        r_Temp       = xyr_Temp       (i,j,:)
        z_VirTemp    = xyz_VirTemp    (i,j,:)
        r_VirTemp    = xyr_VirTemp    (i,j,:)
        r_Press      = xyr_Press      (i,j,:)
        SurfHeight   = xy_SurfHeight  (i,j)
        z_Height     = xyz_Height     (i,j,:)
        r_Height     = xyr_Height     (i,j,:)
        z_Exner      = xyz_Exner      (i,j,:)
        r_Exner      = xyr_Exner      (i,j,:)
        z_TurKinEne  = xyz_TurKinEne  (i,j,:)
        SurfMomFluxX = xy_SurfMomFluxX(i,j)
        SurfMomFluxY = xy_SurfMomFluxY(i,j)


        call VDiffusionMY251D(                               &
          & z_U, z_V, zf_QMix,                               & ! (in)
          & z_Temp, r_Temp, z_VirTemp, r_VirTemp, r_Press,   & ! (in)
          & SurfHeight,                                      & ! (in)
          & z_Height, r_Height, z_Exner, r_Exner,            & ! (in)
          & z_TurKinEne,                                     & ! (in)
          & SurfMomFluxX, SurfMomFluxY,                      & ! (in)
          & r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux, & ! (out)
          & r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef,   & ! (out)
          & z_DTurKinEneDt                                   & ! (out)
          & )

        xyr_MomFluxX     (i,j,:)   = r_MomFluxX
        xyr_MomFluxY     (i,j,:)   = r_MomFluxY
        xyr_HeatFlux     (i,j,:)   = r_HeatFlux
        xyrf_QMixFlux    (i,j,:,:) = rf_QMixFlux
        xyr_VelDiffCoef (i,j,:) = r_VelDiffCoef
        xyr_TempDiffCoef(i,j,:) = r_TempDiffCoef
        xyr_QMixDiffCoef(i,j,:) = r_QMixDiffCoef
        xyz_DTurKinEneDt (i,j,:)   = z_DTurKinEneDt

      end do
    end do



!!$    ! 拡散係数の出力
!!$    ! Output diffusion coefficients
!!$    !
!!$
!!$    ! 拡散係数出力
!!$    ! Diffusion coeffficients output
!!$    !
!!$    call HistoryAutoPut( TimeN, 'VelDiffCoef',  xyr_VelDiffCoef  )
!!$    call HistoryAutoPut( TimeN, 'TempDiffCoef', xyr_TempDiffCoef )
!!$    call HistoryAutoPut( TimeN, 'QVapDiffCoef', xyr_QMixDiffCoef )
!!$
!!$    call HistoryAutoPut( TimeN, 'TKEPShear', xyz_TurKinEneProShear )
!!$    call HistoryAutoPut( TimeN, 'TKEPBuoy' , xyz_TurKinEneProBuoya )
!!$    call HistoryAutoPut( TimeN, 'TKEDiss'  , xyz_TurKinEneDiss     )
!!$
!!$    call HistoryAutoPut( TimeN, 'MixLength' , xyz_MixLength )



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

  end subroutine VDiffusionMY251DWrapper3D

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

  subroutine VDiffusionMY251D(                         &
    & z_U, z_V, zf_QMix,                               & ! (in)
    & z_Temp, r_Temp, z_VirTemp, r_VirTemp, r_Press,   & ! (in)
    & SurfHeight,                                      & ! (in)
    & z_Height, r_Height, z_Exner, r_Exner,            & ! (in)
    & z_TurKinEne,                                     & ! (in)
    & SurfMomFluxX, SurfMomFluxY,                      & ! (in)
    & r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux, & ! (out)
    & r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef,   & ! (out)
    & z_DTurKinEneDt                                   & ! (out)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated by use of MY2.5 model.
    !

    ! モジュール引用 ; 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

    ! 時刻管理
    ! 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

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: z_U (1:kmax)
                              ! $ u $ .   東西風速. Eastward wind
    real(DP), intent(in):: z_V (1:kmax)
                              ! $ v $ .   南北風速. Northward wind
    real(DP), intent(in):: zf_QMix(1:kmax, 1:ncmax)
                              ! $ q $ .   質量混合比. Mass mixing ratio
    real(DP), intent(in):: z_Temp (1:kmax)
                              ! $ T $ .   温度. Temperature
    real(DP), intent(in):: r_Temp (0:kmax)
                              ! $ \hat{T} $ . 温度 (半整数レベル). 
                              ! Temperature (half level)
    real(DP), intent(in):: z_VirTemp (1:kmax)
                              ! $ T_v $ .   仮温度. Virtual temperature
    real(DP), intent(in):: r_VirTemp (0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual temperature (half level)
    real(DP), intent(in):: r_Press  (0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: SurfHeight
                              ! $ z_s $ . 地表面高度. 
                              ! Surface height. 
    real(DP), intent(in):: z_Height (1:kmax)
                              ! 高度 (整数レベル). 
                              ! Height (full level)
    real(DP), intent(in):: r_Height (0:kmax)
                              ! 高度 (半整数レベル). 
                              ! Height (half level)
    real(DP), intent(in):: z_Exner (1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: r_Exner (0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in):: z_TurKinEne(1:kmax)
                              ! 
                              ! Turbulent kinetic energy (m2 s-2)
    real(DP), intent(in):: SurfMomFluxX
                              ! 
                              ! Eastward momentum flux at surface
    real(DP), intent(in):: SurfMomFluxY
                              ! 
                              ! Northward momentum flux at surface

    real(DP), intent(out):: r_MomFluxX (0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward momentum flux
    real(DP), intent(out):: r_MomFluxY (0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP), intent(out):: r_HeatFlux (0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(out):: rf_QMixFlux(0:kmax, 1:ncmax)
                              ! 質量フラックス. 
                              ! Mass flux of compositions
    real(DP), intent(out):: r_VelDiffCoef (0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: r_TempDiffCoef(0:kmax)
                              ! 拡散係数：温度. 
                              ! Diffusion coefficient: temperature
    real(DP), intent(out):: r_QMixDiffCoef(0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity
    real(DP), intent(out):: z_DTurKinEneDt (1:kmax)
                              !
                              ! Tendency of turbulent kinetic energy

    ! 作業変数
    ! Work variables
    !

    real(DP) :: z_MixLength(1:kmax)
                              ! 混合距離. 
                              ! Mixing length
    real(DP) :: z_DVelDzSq(1:kmax)
                              !
                              ! Vertical shear squared (s-2)
    real(DP) :: z_StatStab(1:kmax)
                              !
                              ! Static stability (s-2)
    real(DP) :: GhMin
                              !
                              ! Minimum of G_h
    real(DP) :: GhMax
                              !
                              ! Maximum of G_h
    real(DP) :: z_Gm(1:kmax)
                              !
                              ! G_m
    real(DP) :: z_Gh(1:kmax)
                              !
                              ! G_h
    real(DP) :: z_Sm(1:kmax)
                              !
                              ! S_M
    real(DP) :: z_Sh(1:kmax)
                              !
                              ! S_h

    real(DP), parameter :: Stke = 0.2_DP
                              !
                              ! S_{TKE} = 0.2

    real(DP) :: z_VelDiffCoef (1:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP) :: z_TempDiffCoef (1:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature

    real(DP) :: r_TurKinEneDiffCoef (0:kmax)
                              ! 
                              ! Diffusion coefficient: turbulent kinetic energy
    real(DP) :: z_TurKinEneDiffCoef (1:kmax)
                              ! 
                              ! Diffusion coefficient: turbulent kinetic energy
    real(DP) :: r_TurKinEneTransCoef(0:kmax)
                              ! 
                              ! Transfer coefficient: turbulent kinetic energy

    real(DP) :: r_TurKinEneFlux(0:kmax)
                              ! 
                              ! Turbulent energy flux

    real(DP) :: z_CShe1(1:kmax)
    real(DP) :: z_CShe2(1:kmax)
    real(DP) :: z_CBuo1(1:kmax)
    real(DP) :: z_CBuo2(1:kmax)
    real(DP) :: z_CDis1(1:kmax)
    real(DP) :: z_CDis2(1:kmax)
    real(DP) :: z_TurKinEneProShear(1:kmax)
    real(DP) :: z_TurKinEneProBuoya(1:kmax)

    real(DP) :: FricVelSq
    real(DP) :: TurKinEneAtLB

    real(DP) :: za_TurKinEneMtx(1:kmax, -1:1)
                              ! 
                              ! Implicit matrix for turbulent kinetic energy
    real(DP) :: z_TurKinEneVec(1:kmax)
                              ! 
                              ! Implicit vector for turbulent kinetic energy

    real(DP) :: aaza_TurKinEneLUMtx  (1:1, 1:1, 1:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP) :: aaz_DelTurKinEneLUVec(1:1, 1:1, 1:kmax)
                              ! 
                              ! Tendency of turbulent kinetic energy

    real(DP) :: z_TurKinEneDiss(1:kmax)
                              !
                              ! Dissipation rate of turbulent kinetic energy (m2 s-3)

    real(DP) :: z_TurKinEneNonZero(1:kmax)
                              ! 
                              ! Turbulent kinetic energy with offset (m2 s-2)

    real(DP), parameter :: TurKinEneOffset  = ( 1.0e-3_DP )**2 / 2.0_DP

    logical :: FlagReCalc
                              !
                              ! Flag for recalculation
    logical :: a_FlagReCalcLocal (1)
    logical :: a_FlagReCalcGlobal(1)
    integer :: iloop
    integer :: nloop

    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_my_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 turbulent kinetic energy with offset
    !
    z_TurKinEneNonZero = z_TurKinEne + TurKinEneOffset

    !
    ! Calculation of vertical shear squared
    do k = 1, kmax
      if ( k == 1 ) then
        z_DVelDzSq(k) =                              &
          &   (   ( z_U(k+1) - z_U(k  ) )**2   &
          &     + ( z_V(k+1) - z_V(k  ) )**2 ) &
          & / ( z_Height(k+1) - z_Height(k  ) )**2
      else if ( k == kmax ) then
        z_DVelDzSq(k) =                               &
          &    (   ( z_U(k  ) - z_U(k-1) )**2   &
          &      + ( z_V(k  ) - z_V(k-1) )**2 ) &
          &  / ( z_Height(k  ) - z_Height(k-1) )**2
      else
        z_DVelDzSq(k) =                              &
          &   (   ( z_U(k+1) - z_U(k-1) )**2   &
          &     + ( z_V(k+1) - z_V(k-1) )**2 ) &
          & / ( z_Height(k+1) - z_Height(k-1) )**2
      end if
    end do
    ! Calculation of static stability
    do k = 1, kmax
      if ( k == 1 ) then
        z_StatStab(k) = &
          &  Grav / ( z_VirTemp(k) / z_Exner(k) )  &
          &  * (   z_VirTemp(k+1) / z_Exner(k+1)   &
          &      - z_VirTemp(k  ) / z_Exner(k  ) ) &
          &  / ( z_Height(k+1) - z_Height(k  ) )
      else if ( k == kmax ) then
        z_StatStab(k) = &
          &  Grav / ( z_VirTemp(k) / z_Exner(k) )  &
          &  * (   z_VirTemp(k  ) / z_Exner(k  )   &
          &      - z_VirTemp(k-1) / z_Exner(k-1) ) &
          &  / ( z_Height(k  ) - z_Height(k-1) )
      else
        z_StatStab(k) = &
          &  Grav / ( z_VirTemp(k) / z_Exner(k) )  &
          &  * (   z_VirTemp(k+1) / z_Exner(k+1)   &
          &      - z_VirTemp(k-1) / z_Exner(k-1) ) &
          &  / ( z_Height(k+1) - z_Height(k-1) )
      end if
    end do

    ! 混合距離の算出
    ! Calculate mixing length
    !
    do k = 1, kmax
      z_MixLength(k) =                            &
        & FKarm * ( z_Height(k) - SurfHeight ) &
        & / (1.0_DP + FKarm * ( z_Height(k) - SurfHeight ) / MixLengthMax )
    end do
    !   Limit mixing length (Galperin et al., 1988) and avoid zero
    z_MixLength = &
      &   min( z_MixLength, &
      &        0.53_DP &
      &          * sqrt( 2.0_DP * z_TurKinEneNonZero / max( z_StatStab, 1.0e-10_DP ) ) ) &
      & + 1.0e-10_DP

    z_Gh = - z_MixLength**2 / ( 2.0_DP * z_TurKinEneNonZero ) &
      & * z_StatStab
    ! Actually, xyz_Gm is not used below.
    z_Gm =   z_MixLength**2 / ( 2.0_DP * z_TurKinEneNonZero ) &
      & * z_DVelDzSq


    ! Limit Gh (Galperin et al., 1988)
    GhMin = - 0.53_DP**2
    GhMax = 1.0_DP                                         &
      & / ( MYConstA2                                      &
      &     * (  12.0_DP * MYConstA1 + MYConstB1 + 3.0_DP * MYConstB2 ) )
    z_Gh = max( GhMin, min( z_Gh, GhMax ) )


    z_Sh =                                                          &
      &     MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) &
      & / (   1.0_DP                                                &
      &     - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * z_Gh )

    z_Sm =                                                          &
      & (                                                           &
      &     MYConstA1 * (   1.0_DP - 3.0_DP * MYConstC1             &
      &                   - 6.0_DP * MYConstA1 / MYConstB1 )        &
      &   + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) &
      &     * z_Gh * z_Sh                                           &
      & )                                                           &
      & / (   1.0_DP                                                &
      &     - 9.0_DP * MYConstA1 * MYConstA2 * z_Gh )



    ! 拡散係数の計算
    ! Calculation of diffusion coefficient
    !
    z_VelDiffCoef  = z_MixLength * sqrt( 2.0_DP * z_TurKinEneNonZero ) * z_Sm
    z_TempDiffCoef = z_MixLength * sqrt( 2.0_DP * z_TurKinEneNonZero ) * z_Sh
    !
    do k = 0, kmax
      if ( ( k == 0 ) .or. ( k == kmax ) ) then
        r_VelDiffCoef (k) = 0.0_DP
        r_TempDiffCoef(k) = 0.0_DP
      else
        r_VelDiffCoef (k) = ( z_VelDiffCoef (k) + z_VelDiffCoef (k+1) ) / 2.0_DP
        r_TempDiffCoef(k) = ( z_TempDiffCoef(k) + z_TempDiffCoef(k+1) ) / 2.0_DP
      end if
    end do
    !
    do k = 1, kmax-1
      r_VelDiffCoef(k) =                                &
        & max( min( r_VelDiffCoef(k), VelDiffCoefMax ), &
        &      VelDiffCoefMin                                 &
        &    )
      r_TempDiffCoef(k) =                                 &
        & max( min( r_TempDiffCoef(k), TempDiffCoefMax ), &
        &      TempDiffCoefMin                                  &
        &    )
    end do
    !
    r_QMixDiffCoef      = r_TempDiffCoef


    ! 輸送係数とフラックスの計算
    ! Calculate transfer coefficient and flux
    !
    call VDiffusionCalcFlux1D(                           &
      & z_U, z_V, zf_QMix,                               & ! (in)
      & z_Temp, r_VirTemp, r_Press,                      & ! (in)
      & z_Height, z_Exner, r_Exner,                      & ! (in)
      & r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef,   & ! (in)
      & r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux  & ! (out)
      & )


    ! Calculate tendency of turbulent kinetic energy

    !   Set diffusion coefficient for turbulent kinetic energy
    z_TurKinEneDiffCoef = &
      & z_MixLength * sqrt( 2.0_DP * z_TurKinEneNonZero ) * Stke
    !
    do k = 0, kmax
      if ( k == 0 ) then
        r_TurKinEneDiffCoef(k) = z_TurKinEneDiffCoef(1)
      else if ( k == kmax ) then
        r_TurKinEneDiffCoef(k) = z_TurKinEneDiffCoef(kmax)
      else
        r_TurKinEneDiffCoef(k) = &
          & ( z_TurKinEneDiffCoef(k) + z_TurKinEneDiffCoef(k+1) ) &
          & / 2.0_DP
      end if
    end do

    !   Calculate turbulent kinetic energy at lower boundary
    !
    FricVelSq = sqrt( SurfMomFluxX**2 + SurfMomFluxY**2 ) &
      & / ( r_Press(0) / ( GasRDry * r_VirTemp(0) ) )
    TurKinEneAtLB = MYConstB1**(2.0_DP/3.0_DP) / 2.0_DP * FricVelSq
    TurKinEneAtLB = TurKinEneAtLB + TurKinEneOffset

    !   Calculate transfer coefficient and flux of turbulent kinetic energy
    !
    !    When transfer coefficient at lower boundary is calculated, 
    !    diffusion coefficient at mid-point of 1st layer is used. 
    !    In addition, transfer coefficient at upper boundary is assumed 
    !    to be zero.
    k = 0
    r_TurKinEneTransCoef(k) =                               &
      &   z_TurKinEneDiffCoef(1)                            &
      &     * r_Press(k) / ( GasRDry * r_VirTemp(k) ) &
      &     / ( z_Height(1) - SurfHeight )
    do k = 1, kmax-1
      r_TurKinEneTransCoef(k) =                               &
        &   r_TurKinEneDiffCoef(k)                            &
        &     * r_Press(k) / ( GasRDry * r_VirTemp(k) ) &
        &     / ( z_Height(k+1) - z_Height(k) )
    end do
    k = kmax
    r_TurKinEneTransCoef(k) = 0.0_DP
    !
    do k = 1, kmax-1
      r_TurKinEneFlux(k) =               &
        & - r_TurKinEneTransCoef(k)      &
        &   * ( z_TurKinEneNonZero(k+1) - z_TurKinEneNonZero(k) )
    end do
    k = 0
    r_TurKinEneFlux(k) =               &
      & - r_TurKinEneTransCoef(k)      &
      &   * ( z_TurKinEneNonZero(k+1) - TurKinEneAtLB )
    k = kmax
    r_TurKinEneFlux(k) = 0.0_DP


    z_CShe1 =            sqrt( 2.0_DP ) * z_MixLength * z_DVelDzSq  &
      & * z_Sm                                                          &
      & * sqrt( z_TurKinEneNonZero )
    z_CShe2 = 0.0_DP
    z_CBuo1 = -          sqrt( 2.0_DP ) * z_MixLength * z_StatStab  &
      & * z_Sh                                                          &
      & * sqrt( z_TurKinEneNonZero )
    z_CBuo2 = 0.0_DP
    z_CDis1 =          2.0_DP**1.5_DP / ( MYConstB1 * z_MixLength ) &
      & * z_TurKinEneNonZero**1.5_DP
    z_CDis2 = 1.5_DP * 2.0_DP**1.5_DP / ( MYConstB1 * z_MixLength ) &
      & * z_TurKinEneNonZero**0.5_DP


    nloop = kmax
    loop_solve : do iloop = 1, nloop

      !
      ! Construct implicit matrix from transfer coefficient of vertical 
      ! diffusion scheme (turbulent kinetic energy)
      !
      k = 1
      za_TurKinEneMtx(k,-1) = 0.0_DP
      za_TurKinEneMtx(k, 0) =                                             &
        & - ( r_Press(k) - r_Press(k-1) ) / Grav / ( 2.0_DP * DelTime ) &
        & + r_TurKinEneTransCoef(k-1)                                     &
        & + r_TurKinEneTransCoef(k  )                                     &
        & + ( r_Press(k) - r_Press(k-1) ) / Grav                    &
        &     * ( z_CShe2(k) + z_CBuo2(k) - z_CDis2(k) )
      za_TurKinEneMtx(k, 1) =                                             &
        & - r_TurKinEneTransCoef(k  )
      !
      do k = 2, kmax-1
        za_TurKinEneMtx(k,-1) =&
          & - r_TurKinEneTransCoef(k-1)
        za_TurKinEneMtx(k, 0) =                                             &
          & - ( r_Press(k) - r_Press(k-1) ) / Grav / ( 2.0_DP * DelTime ) &
          & + r_TurKinEneTransCoef(k-1)                                     &
          & + r_TurKinEneTransCoef(k  )                                     &
          & + ( r_Press(k) - r_Press(k-1) ) / Grav                    &
          &     * ( z_CShe2(k) + z_CBuo2(k) - z_CDis2(k) )
        za_TurKinEneMtx(k, 1) =                                             &
          & - r_TurKinEneTransCoef(k  )
      end do
      !
      k = kmax
      za_TurKinEneMtx(k,-1) =&
        & - r_TurKinEneTransCoef(k-1)
      za_TurKinEneMtx(k, 0) =                                             &
        & - ( r_Press(k) - r_Press(k-1) ) / Grav / ( 2.0_DP * DelTime ) &
        & + r_TurKinEneTransCoef(k-1)                                     &
        & + r_TurKinEneTransCoef(k  )                                     &
        & + ( r_Press(k) - r_Press(k-1) ) / Grav                    &
        &     * ( z_CShe2(k) + z_CBuo2(k) - z_CDis2(k) )
      za_TurKinEneMtx(k, 1) = 0.0_DP

      do k = 1, kmax
        z_TurKinEneVec(k) =                                               &
          & - ( r_TurKinEneFlux(k) - r_TurKinEneFlux(k-1) )         &
          & - ( r_Press(k) - r_Press(k-1) ) / Grav                  &
          &     * ( z_CShe1(k) + z_CBuo1(k) - z_CDis1(k) )
      end do


      !
      ! Solve simultaneous linear equations by use of LU decomposition technique
      !
      aaza_TurKinEneLUMtx(1,1,:,:) = za_TurKinEneMtx
      !
      call PhyImplLUDecomp3(   &
        & aaza_TurKinEneLUMtx, &   ! (inout)
        & 1 * 1, kmax          &   ! (in)
        & )

      aaz_DelTurKinEneLUVec(1,1,:) = z_TurKinEneVec
      !
      call PhyImplLUSolve3(      &
        & aaz_DelTurKinEneLUVec, & ! (inout)
        & aaza_TurKinEneLUMtx,   & ! (in)
        & 1, 1 * 1 , kmax        & ! (in)
        & )

      z_DTurKinEneDt = aaz_DelTurKinEneLUVec(1,1,:) / ( 2.0_DP * DelTime )


      ! Calculation of dissipation rate of turbulent kinetic energy
      !
      ! Calculate production rate of turbulent kinetic energy
      ! by shear and buoyancy
      z_TurKinEneProShear = z_CShe1 + z_CShe2 * z_DTurKinEneDt * 2.0_DP * DelTime
      z_TurKinEneProBuoya = z_CBuo1 + z_CBuo2 * z_DTurKinEneDt * 2.0_DP * DelTime
      z_TurKinEneDiss     = z_CDis1 + z_CDis2 * z_DTurKinEneDt * 2.0_DP * DelTime

      ! Check of turbulent kinetic energy dissipation rate
      ! If it is negative, tendency is recalculated without dissipation.
      !
      FlagReCalc = .false.
      do k = 1, kmax
        if ( z_TurKinEneDiss(k) < 0.0_DP ) then
          z_CDis1(k) = 0.0_DP
          z_CDis2(k) = 0.0_DP
          FlagReCalc = .true.
        end if
      end do

      ! Check convergence
      a_FlagReCalcLocal = FlagReCalc
!!$      call MPIWrapperChkTrue(   &
!!$        & 1, a_FlagReCalcLocal, & ! (in)
!!$        & a_FlagReCalcGlobal    & ! (out)
!!$        & )
      a_FlagReCalcGlobal = a_FlagReCalcLocal
      if ( .not. a_FlagReCalcGlobal(1) ) exit loop_solve

    end do loop_solve



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

  end subroutine VDiffusionMY251D

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

  subroutine VDiffusionMY25GBT94(                                &
    & xyz_U, xyz_V, xyzf_QMix,                                   & ! (in)
    & xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press,   & ! (in)
    & xy_SurfHeight,                                             & ! (in)
    & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner,              & ! (in)
    & xyz_TurKinEne,                                             & ! (in)
    & xy_SurfMomFluxX, xy_SurfMomFluxY,                          & ! (in)
    & xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, & ! (out)
    & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef,       & ! (in)
    & xyz_DTurKinEneDt                                           & ! (out)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated by use of MY2.5 model.
    !

    ! モジュール引用 ; 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

    ! 時刻管理
    ! 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

    ! MPI 関連ルーチン
    ! MPI related routines
    !
    use mpi_wrapper, only : MPIWrapperChkTrue

    ! 陰解法による時間積分のためのルーチン
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3

    ! 宣言文 ; 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):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T_v $ .   仮温度. Virtual temperature
    real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual 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(in):: xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Turbulent kinetic energy (m2 s-2)
    real(DP), intent(in):: xy_SurfMomFluxX (0:imax-1, 1:jmax)
                              ! 
                              ! Eastward momentum flux at surface
    real(DP), intent(in):: xy_SurfMomFluxY (0:imax-1, 1:jmax)
                              ! 
                              ! Northward momentum flux at surface

    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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Diffusion coefficient: temperature
    real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity

    real(DP), intent(out):: xyz_DTurKinEneDt (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Tendency of turbulent kinetic energy

    ! 作業変数
    ! Work variables
    !

    real(DP) :: xyz_MixLength(0:imax-1, 1:jmax, 1:kmax)
                              ! 混合距離. 
                              ! Mixing length
    real(DP) :: xyz_DVelDzSq(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Vertical shear squared (s-2)
    real(DP) :: xyz_StatStab(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Static stability (s-2)
!!$    real(DP) :: GhMin
!!$                              !
!!$                              ! Minimum of G_h
!!$    real(DP) :: GhMax
!!$                              !
!!$                              ! Maximum of G_h
    real(DP) :: xyz_Gm(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! G_m
    real(DP) :: xyz_Gh(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! G_h
!!$    real(DP) :: xyz_SmHat(0:imax-1, 1:jmax, 1:kmax)
!!$                              !
!!$                              ! \hat{S}_M
!!$    real(DP) :: xyz_ShHat(0:imax-1, 1:jmax, 1:kmax)
!!$                              !
!!$                              ! \hat{S}_h
!!$    real(DP) :: xyz_DSmHatDTKE(0:imax-1, 1:jmax, 1:kmax)
!!$                              !
!!$                              ! derivative of \hat{S}_M
!!$    real(DP) :: xyz_DShHatDTKE(0:imax-1, 1:jmax, 1:kmax)
!!$                              !
!!$                              ! derivative of \hat{S}_h
    real(DP) :: xyz_Sm(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! S_M
    real(DP) :: xyz_Sh(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! S_h

    real(DP), parameter :: Stke = 0.2_DP
                              !
                              ! S_{TKE} = 0.2

    real(DP) :: xyz_VelDiffCoef (0:imax-1, 1:jmax, 1:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP) :: xyz_TempDiffCoef (0:imax-1, 1:jmax, 1:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature

    real(DP) :: xyr_TurKinEneDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Diffusion coefficient: turbulent kinetic energy
    real(DP) :: xyz_TurKinEneDiffCoef (0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Diffusion coefficient: turbulent kinetic energy
    real(DP) :: xyr_TurKinEneTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Transfer coefficient: turbulent kinetic energy

    real(DP) :: xyr_TurKinEneFlux(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Turbulent energy flux

    real(DP) :: xy_FricVelSq    (0:imax-1, 1:jmax)
    real(DP) :: xy_TurKinEneAtLB(0:imax-1, 1:jmax)

    real(DP) :: xyza_TurKinEneMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! 
                              ! Implicit matrix for turbulent kinetic energy
    real(DP) :: xyz_TurKinEneVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Implicit vector for turbulent kinetic energy

    real(DP) :: xyza_TurKinEneLUMtx  (0:imax-1, 1:jmax, 1:kmax, -1:1)
                              ! LU 行列.
                              ! LU matrix
    real(DP) :: xyz_DelTurKinEneLUVec(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Tendency of turbulent kinetic energy

    real(DP) :: xyz_TurKinEneDiss(0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Dissipation rate of turbulent kinetic energy (m2 s-3)

    real(DP) :: xyz_TurKinEneNonZero(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Turbulent kinetic energy with offset (m2 s-2)

    real(DP), parameter :: TurKinEneOffset  = ( 1.0e-3_DP )**2 / 2.0_DP


    logical :: xyz_FlagUseRiNum       (0:imax-1, 1:jmax, 1:kmax)
    logical :: xyz_FlagTKEAsymptToZero(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_RiNum(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: Beta1
    real(DP) :: Beta2
    real(DP) :: Beta3
    real(DP) :: Beta4
    real(DP) :: MixLength
    real(DP) :: xyz_A1    (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_A2    (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_R1    (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_R2    (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_SqrtA1(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_TKEInit     (0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Turbulent kinetic energy for current time step
                              ! with offset (m2 s-2)
    real(DP) :: xyz_TKEx2       (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_SqrtTKEx2   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_TKETentative(0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: Alpha
    real(DP) :: BetaSq
    real(DP) :: StatStab
    real(DP) :: RiNum
    real(DP) :: DVelDzSq

    real(DP), parameter :: Epsilon   = 1.0e-10_DP
    real(DP), parameter :: CrtlRiNum = 0.195_DP
    real(DP), parameter :: CrtlShear = 0.001_DP / 1000.0_DP

!!$    logical :: FlagReCalc
!!$                              !
!!$                              ! Flag for recalculation
!!$    logical :: a_FlagReCalcLocal (1)
!!$    logical :: a_FlagReCalcGlobal(1)
!!$    integer :: iloop
!!$    integer :: nloop

    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

    ! 実行文 ; Executable statement
    !

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


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


    !
    ! Calculation of vertical shear squared
    do k = 1, kmax
      if ( k == 1 ) then
        xyz_DVelDzSq(:,:,k) =                              &
          &   (   ( xyz_U(:,:,k+1) - xyz_U(:,:,k  ) )**2   &
          &     + ( xyz_V(:,:,k+1) - xyz_V(:,:,k  ) )**2 ) &
          & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k  ) )**2
!!$        xyz_DVelDzSq(:,:,k) =                              &
!!$          &   (   ( xyz_U(:,:,k+1) - 0.0_DP )**2   &
!!$          &     + ( xyz_V(:,:,k+1) - 0.0_DP )**2 ) &
!!$          & / ( xyz_Height(:,:,k+1) - xy_SurfHeight )**2
      else if ( k == kmax ) then
        xyz_DVelDzSq(:,:,k) =                               &
          &    (   ( xyz_U(:,:,k  ) - xyz_U(:,:,k-1) )**2   &
          &      + ( xyz_V(:,:,k  ) - xyz_V(:,:,k-1) )**2 ) &
          &  / ( xyz_Height(:,:,k  ) - xyz_Height(:,:,k-1) )**2
      else
        xyz_DVelDzSq(:,:,k) =                              &
          &   (   ( xyz_U(:,:,k+1) - xyz_U(:,:,k-1) )**2   &
          &     + ( xyz_V(:,:,k+1) - xyz_V(:,:,k-1) )**2 ) &
          & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) )**2
      end if
    end do
    ! Calculation of static stability
    do k = 1, kmax
      if ( k == 1 ) then
        xyz_StatStab(:,:,k) = &
          &  Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) )  &
          &  * (   xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1)   &
          &      - xyz_VirTemp(:,:,k  ) / xyz_Exner(:,:,k  ) ) &
          &  / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k  ) )
      else if ( k == kmax ) then
        xyz_StatStab(:,:,k) = &
          &  Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) )  &
          &  * (   xyz_VirTemp(:,:,k  ) / xyz_Exner(:,:,k  )   &
          &      - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) &
          &  / ( xyz_Height(:,:,k  ) - xyz_Height(:,:,k-1) )
      else
        xyz_StatStab(:,:,k) = &
          &  Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) )  &
          &  * (   xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1)   &
          &      - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) &
          &  / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) )
      end if
    end do


    ! 混合距離の算出
    ! Calculate mixing length
    !
    do k = 1, kmax
      xyz_MixLength(:,:,k) =                            &
        & FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) &
        & / (1.0_DP + FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) / MixLengthMax )
    end do
!!$    !   Limit mixing length (Galperin et al., 1988) and avoid zero
!!$    xyz_MixLength = &
!!$      &   min( xyz_MixLength, &
!!$      &        0.53_DP &
!!$      &          * sqrt( 2.0_DP * xyz_TurKinEneNonZero / max( xyz_StatStab, 1.0d-10 ) ) ) &
!!$      & + 1.0d-10

    !********************************************************************
    ! Gerrity et al. (1994)

    ! Set flag for using Richardson number
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_DVelDzSq(i,j,k) < CrtlShear**2 ) then
            xyz_FlagUseRiNum(i,j,k) = .false.
          else
            xyz_FlagUseRiNum(i,j,k) = .true.
          end if
        end do
      end do
    end do
    ! Calculation of Richardson number
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FlagUseRiNum(i,j,k) ) then
            xyz_RiNum(i,j,k) = xyz_StatStab(i,j,k) / xyz_DVelDzSq(i,j,k)
          else
            xyz_RiNum(i,j,k) = - 1.0e100_DP
          end if
        end do
      end do
    end do

    ! Set flag for selecting an asymptotic equation
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FlagUseRiNum(i,j,k) ) then
            if ( CrtlRiNum <= xyz_RiNum(i,j,k) ) then
              ! Ric <= Ri
              xyz_FlagTKEAsymptToZero = .true.
            else
              ! Ri  <  Ric
              xyz_FlagTKEAsymptToZero = .false.
            end if
          else
            if ( xyz_StatStab(i,j,k) > 0.0_DP ) then
              xyz_FlagTKEAsymptToZero = .true.
            else
              xyz_FlagTKEAsymptToZero = .false.
            end if
          end if
        end do
      end do
    end do

    ! Calculation of roos for equations for steady state
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FlagUseRiNum(i,j,k) ) then
            ! Calculation with Richardson number
            MixLength = xyz_MixLength(i,j,k)
            DVelDzSq  = xyz_DVelDzSq (i,j,k)
            RiNum     = xyz_RiNum    (i,j,k)
            !
            Beta1 =                              &
              & - MixLength**2 * DVelDzSq        &
              &   * (  6.53_DP - 49.0_DP  * RiNum )
            Beta2 =                              &
              & - MixLength**4 * DVelDzSq**2     &
              &   * ( 51.2_DP  - 262.7_DP * RiNum ) * RiNum
            Beta3 =                              &
              &   MixLength**2 * DVelDzSq        &
              &   * (  5.08_DP + 36.7_DP  * RiNum )
            Beta4 =                              &
              &   MixLength**4 * DVelDzSq**2     &
              &   * ( 88.8_DP  + 187.4_DP * RiNum ) * RiNum
            !
            xyz_A1(i,j,k) = - Beta1 / 2.0_DP &
              &             + sqrt( Beta1**2 - 4.0_DP * Beta2 ) / 2.0_DP
            xyz_A2(i,j,k) = - Beta1 / 2.0_DP &
              &             - sqrt( Beta1**2 - 4.0_DP * Beta2 ) / 2.0_DP
            xyz_R1(i,j,k) = - Beta3 / 2.0_DP &
              &             + sqrt( Beta3**2 - 4.0_DP * Beta4 ) / 2.0_DP
            xyz_R2(i,j,k) = - Beta3 / 2.0_DP &
              &             - sqrt( Beta3**2 - 4.0_DP * Beta4 ) / 2.0_DP
          else
            ! Calculation without Richardson number
            MixLength = xyz_MixLength(i,j,k)
            StatStab  = xyz_StatStab (i,j,k)
            !
            xyz_A1(i,j,k) = 24.49_DP * MixLength**2 *    ( - StatStab ) &
              &           + 18.36_DP * MixLength**2 * abs( - StatStab )
            xyz_A2(i,j,k) = 24.49_DP * MixLength**2 *    ( - StatStab ) &
              &           - 18.36_DP * MixLength**2 * abs( - StatStab )
            xyz_R1(i,j,k) = 18.35_DP * MixLength**2 *    ( - StatStab ) &
              &           + 12.22_DP * MixLength**2 * abs( - StatStab )
            xyz_R2(i,j,k) = 18.35_DP * MixLength**2 *    ( - StatStab ) &
              &           - 12.22_DP * MixLength**2 * abs( - StatStab )
          end if
        end do
      end do
    end do

    ! Set turbulent kinetic energy at current time step
    !
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FlagUseRiNum(i,j,k) ) then
            if ( xyz_RiNum(i,j,k) < 0.0_DP ) then
              ! Ri < 0 (or CT > 0)
              xyz_TKEInit(i,j,k) = &
                & max( xyz_TurKinEne(i,j,k), xyz_R1(i,j,k) / 2.0_DP )
            else if ( xyz_RiNum(i,j,k) < CrtlRiNum ) then
              ! 0 <= Ri < Ric
              xyz_TKEInit(i,j,k) = xyz_TurKinEne(i,j,k)
            else
              ! Ric <= Ri
              xyz_TKEInit(i,j,k) = xyz_TurKinEne(i,j,k)
            end if
          else
            if ( xyz_StatStab(i,j,k) < 0.0_DP ) then
              ! CT >  0
              xyz_TKEInit(i,j,k) = &
                & max( xyz_TurKinEne(i,j,k), xyz_R1(i,j,k) / 2.0_DP )
            else
              ! CT <= 0
              xyz_TKEInit(i,j,k) = xyz_TurKinEne(i,j,k)
            end if
          end if
          xyz_TKEInit(i,j,k) = xyz_TKEInit(i,j,k) + TurKinEneOffset
        end do
      end do
    end do
    !
    xyz_TKEx2     = 2.0_DP * xyz_TKEInit
    xyz_SqrtTKEx2 = sqrt( xyz_TKEx2 )
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FlagTKEAsymptToZero(i,j,k) ) then
            xyz_SqrtA1(i,j,k) = 1.0e100_DP
          else
            xyz_SqrtA1(i,j,k) = sqrt( xyz_A1(i,j,k) )
          end if
        end do
      end do
    end do
    !
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FlagTKEAsymptToZero(i,j,k) ) then
            ! Use equations asymptotic to zero TKE, Eq. (10)
            ! Eq. (11)
            BetaSq =                               &
              &   ( xyz_TKEx2(i,j,k) - xyz_A1(i,j,k)**2 )     &
              & * ( xyz_TKEx2(i,j,k) - xyz_A2(i,j,k)**2 )     &
              & / ( MYConstB1 * xyz_MixLength(i,j,k)          &
              &     * ( xyz_TKEx2(i,j,k) - xyz_R1(i,j,k) )    &
              &     * ( xyz_TKEx2(i,j,k) - xyz_R2(i,j,k) ) )
            ! Eq. (10)
            xyz_TKETentative(i,j,k) =                   &
              & xyz_TKEInit(i,j,k) / 2.0_DP             &
              & / ( 1.0_DP                              &
              &     + xyz_SqrtTKEx2(i,j,k) * BetaSq     &
              &       * ( 2.0_DP * DelTime ) )**2
          else
            ! Use equations asymptotic to certain TKE, Eq. (8)
            ! Eq. (7b)
            BetaSq =                                          &
              &   xyz_TKEx2(i,j,k)                            &
              & * ( xyz_TKEx2(i,j,k) - xyz_A2(i,j,k)**2 )     &
              & / ( MYConstB1 * xyz_MixLength(i,j,k)          &
              &     * ( xyz_TKEx2(i,j,k) - xyz_R1(i,j,k) )    &
              &     * ( xyz_TKEx2(i,j,k) - xyz_R2(i,j,k) ) )
            ! Eq. (9)
            Alpha =                                            &
              &   ( xyz_SqrtTKEx2(i,j,k) - xyz_SqrtA1(i,j,k) ) &
              & / ( xyz_SqrtTKEx2(i,j,k) + xyz_SqrtA1(i,j,k) ) &
              & * exp( - 2.0_DP * xyz_SqrtA1(i,j,k) * BetaSq   &
              &          * ( 2.0_DP * DelTime ) )
            Alpha = max( min( Alpha, 1.0_DP - Epsilon ), -1.0_DP )
            ! Eq. (8)
            xyz_TKETentative(i,j,k) =                 &
              & xyz_A1(i,j,k) / 2.0_DP                &
              &   * ( ( 1.0_DP + Alpha ) / ( 1.0_DP - Alpha ) )**2
          end if
        end do
      end do
    end do


    !********************************************************************

    ! Set turbulent kinetic energy for diffusion calculation
    !
!!$    xyz_TurKinEneNonZero = xyz_TurKinEne + TurKinEneOffset
    xyz_TurKinEneNonZero = xyz_TKETentative + TurKinEneOffset


    xyz_Gh = - xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) &
      & * xyz_StatStab
    ! Actually, xyz_Gm is not used below.
    xyz_Gm =   xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) &
      & * xyz_DVelDzSq


    ! Limit Gh (Galperin et al., 1988)
!!$    GhMin = - 0.53d0**2
!!$    GhMax = 1.0_DP                                         &
!!$      & / ( MYConstA2                                      &
!!$      &     * (  12.0_DP * MYConstA1 + MYConstB1 + 3.0_DP * MYConstB2 ) )
!!$    xyz_Gh = max( GhMin, min( xyz_Gh, GhMax ) )


    xyz_Sh =                                                        &
      &     MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) &
      & / (   1.0_DP                                                &
      &     - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * xyz_Gh )

    xyz_Sm =                                                        &
      & (                                                           &
      &     MYConstA1 * (   1.0_DP - 3.0_DP * MYConstC1             &
      &                   - 6.0_DP * MYConstA1 / MYConstB1 )        &
      &   + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) &
      &     * xyz_Gh * xyz_Sh                                       &
      & )                                                           &
      & / (   1.0_DP                                                &
      &     - 9.0_DP * MYConstA1 * MYConstA2 * xyz_Gh )


    ! 拡散係数の計算
    ! Calculation of diffusion coefficient
    !
    xyz_VelDiffCoef  = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sm
    xyz_TempDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sh
    !
    do k = 0, kmax
      if ( ( k == 0 ) .or. ( k == kmax ) ) then
        xyr_VelDiffCoef (:,:,k) = 0.0_DP
        xyr_TempDiffCoef(:,:,k) = 0.0_DP
      else
        xyr_VelDiffCoef (:,:,k) = ( xyz_VelDiffCoef (:,:,k) + xyz_VelDiffCoef (:,:,k+1) ) / 2.0_DP
        xyr_TempDiffCoef(:,:,k) = ( xyz_TempDiffCoef(:,:,k) + xyz_TempDiffCoef(:,:,k+1) ) / 2.0_DP
      end if
    end do
    !
    do k = 1, kmax-1
      do j = 1, jmax
        do i = 0, imax-1
          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                                  &
            &    )
        end do
      end do
    end do
    !
    xyr_QMixDiffCoef      = xyr_TempDiffCoef


    ! 輸送係数とフラックスの計算
    ! Calculate transfer coefficient and flux
    !
    call VDiffusionCalcFlux(                                       &
      & xyz_U, xyz_V, xyzf_QMix,                                   & ! (in)
      & xyz_Temp, xyr_VirTemp, xyr_Press,                          & ! (in)
      & xyz_Height, xyz_Exner, xyr_Exner,                          & ! (in)
      & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef,       & ! (in)
      & xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux  & ! (out)
      & )


    ! Calculate tendency of turbulent kinetic energy

    !   Set diffusion coefficient for turbulent kinetic energy
    xyz_TurKinEneDiffCoef = &
      & xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * Stke
    !
    do k = 0, kmax
      if ( k == 0 ) then
        xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,1)
      else if ( k == kmax ) then
        xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,kmax)
      else
        xyr_TurKinEneDiffCoef(:,:,k) = &
          & ( xyz_TurKinEneDiffCoef(:,:,k) + xyz_TurKinEneDiffCoef(:,:,k+1) ) &
          & / 2.0_DP
      end if
    end do

    !   Calculate turbulent kinetic energy at lower boundary
    !
    xy_FricVelSq = sqrt( xy_SurfMomFluxX**2 + xy_SurfMomFluxY**2 ) &
      & / ( xyr_Press(:,:,0) / ( GasRDry * xyr_VirTemp(:,:,0) ) )
    xy_TurKinEneAtLB = MYConstB1**(2.0_DP/3.0_DP) / 2.0_DP * xy_FricVelSq
    xy_TurKinEneAtLB = xy_TurKinEneAtLB + TurKinEneOffset

    !   Calculate transfer coefficient and flux of turbulent kinetic energy
    !
    !    When transfer coefficient at lower boundary is calculated, 
    !    diffusion coefficient at mid-point of 1st layer is used. 
    !    In addition, transfer coefficient at upper boundary is assumed 
    !    to be zero.
    k = 0
    xyr_TurKinEneTransCoef(:,:,k) =                               &
      &   xyz_TurKinEneDiffCoef(:,:,1)                            &
      &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
      &     / ( xyz_Height(:,:,1) - xy_SurfHeight )
    do k = 1, kmax-1
      xyr_TurKinEneTransCoef(:,:,k) =                               &
        &   xyr_TurKinEneDiffCoef(:,:,k)                            &
        &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
        &     / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
    end do
    k = kmax
    xyr_TurKinEneTransCoef(:,:,k) = 0.0_DP
    !
    do k = 1, kmax-1
      xyr_TurKinEneFlux(:,:,k) =               &
        & - xyr_TurKinEneTransCoef(:,:,k)      &
        &   * ( xyz_TurKinEneNonZero(:,:,k+1) - xyz_TurKinEneNonZero(:,:,k) )
    end do
    k = 0
    xyr_TurKinEneFlux(:,:,k) =               &
      & - xyr_TurKinEneTransCoef(:,:,k)      &
      &   * ( xyz_TurKinEneNonZero(:,:,k+1) - xy_TurKinEneAtLB )
    k = kmax
    xyr_TurKinEneFlux(:,:,k) = 0.0_DP


    !
    ! Construct implicit matrix from transfer coefficient of vertical 
    ! diffusion scheme (turbulent kinetic energy)
    !
    k = 1
    xyza_TurKinEneMtx(:,:,k,-1) = 0.0_DP
    xyza_TurKinEneMtx(:,:,k, 0) =                                             &
      & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + xyr_TurKinEneTransCoef(:,:,k-1)                                     &
      & + xyr_TurKinEneTransCoef(:,:,k  )
    xyza_TurKinEneMtx(:,:,k, 1) =                                             &
      & - xyr_TurKinEneTransCoef(:,:,k  )
    !
    do k = 2, kmax-1
      xyza_TurKinEneMtx(:,:,k,-1) =&
        & - xyr_TurKinEneTransCoef(:,:,k-1)
      xyza_TurKinEneMtx(:,:,k, 0) =                                             &
        & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
        & + xyr_TurKinEneTransCoef(:,:,k-1)                                     &
        & + xyr_TurKinEneTransCoef(:,:,k  )
      xyza_TurKinEneMtx(:,:,k, 1) =                                             &
        & - xyr_TurKinEneTransCoef(:,:,k  )
    end do
    !
    k = kmax
    xyza_TurKinEneMtx(:,:,k,-1) =&
      & - xyr_TurKinEneTransCoef(:,:,k-1)
    xyza_TurKinEneMtx(:,:,k, 0) =                                             &
      & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) &
      & + xyr_TurKinEneTransCoef(:,:,k-1)                                     &
      & + xyr_TurKinEneTransCoef(:,:,k  )
    xyza_TurKinEneMtx(:,:,k, 1) = 0.0_DP

    do k = 1, kmax
      xyz_TurKinEneVec(:,:,k) =                                               &
        & - ( xyr_TurKinEneFlux(:,:,k) - xyr_TurKinEneFlux(:,:,k-1) )
    end do


    !
    ! Solve simultaneous linear equations by use of LU decomposition technique
    !
    xyza_TurKinEneLUMtx = xyza_TurKinEneMtx
    !
    call PhyImplLUDecomp3(   &
      & xyza_TurKinEneLUMtx, &   ! (inout)
      & imax * jmax, kmax    &   ! (in)
      & )

    xyz_DelTurKinEneLUVec = xyz_TurKinEneVec
    !
    call PhyImplLUSolve3(      &
      & xyz_DelTurKinEneLUVec, & ! (inout)
      & xyza_TurKinEneLUMtx,   & ! (in)
      & 1, imax * jmax , kmax  & ! (in)
      & )

!!$    xyz_DTurKinEneDt = xyz_DelTurKinEneLUVec / ( 2.0_DP * DelTime )

    xyz_TKETentative = xyz_TKETentative + xyz_DelTurKinEneLUVec
    xyz_DTurKinEneDt = &
      & ( xyz_TKETentative - xyz_TurKinEne ) / ( 2.0_DP * DelTime )


!!$    write( 6, * ) TimeN, iloop
!!$    write( 6, * ) xyz_TurKinEne(0,jmax/2+1,1:4)
!!$    write( 6, * ) xyz_TempDiffCoef(0,jmax/2+1,1:4)
!!$    write( 6, * ) xyz_DTurKinEneDt(0,jmax/2+1,1:4)



    ! 拡散係数の出力
    ! Output diffusion coefficients
    !

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

    call HistoryAutoPut( TimeN, 'MixLength' , xyz_MixLength )



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

  end subroutine VDiffusionMY25GBT94

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

  subroutine VDiffusion(                                         &
    & xyz_U, xyz_V, xyzf_QMix,                                   & ! (in)
    & xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, 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_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef        & ! (out)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated by use of MY2 model.
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & FKarm, &
                              ! $ k $ .
                              ! カルマン定数. 
                              ! Karman constant
      & Grav
                              ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration

    ! 時刻管理
    ! 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):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T_v $ .   仮温度. Virtual temperature
    real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual 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_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: 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 $

    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_my_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


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


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

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

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

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

    else

      ! バルク $ R_i $ 数算出
      ! Calculate bulk $ R_i $
      !
      xyr_DVelDz(:,:,0)       = 0.0_DP
      xyr_DVelDz(:,:,kmax)    = 0.0_DP
      xyr_BulkRiNum(:,:,0)    = 0.0_DP
      xyr_BulkRiNum(:,:,kmax) = 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_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 and flux
    !
    call VDiffusionCalcFlux(                                       &
      & xyz_U, xyz_V, xyzf_QMix,                                   & ! (in)
      & xyz_Temp, xyr_VirTemp, xyr_Press,                          & ! (in)
      & xyz_Height, xyz_Exner, xyr_Exner,                          & ! (in)
      & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef,       & ! (in)
      & xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux  & ! (out)
      & )

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

  end subroutine VDiffusion

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

  subroutine VDiffusionCalcFlux(                                 &
    & xyz_U, xyz_V, xyzf_QMix,                                   & ! (in)
    & xyz_Temp, xyr_VirTemp, xyr_Press,                          & ! (in)
    & xyz_Height, xyz_Exner, xyr_Exner,                          & ! (in)
    & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef,       & ! (in)
    & xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux  & ! (out)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated. 
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & FKarm, &
                              ! $ k $ .
                              ! カルマン定数. 
                              ! Karman constant
      & 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

    ! 時刻管理
    ! 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_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual 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):: 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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(in):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity

    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

    ! 作業変数
    ! Work variables
    !
    real(DP) :: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP) :: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP) :: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass (composition)

    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_my_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! 輸送係数の計算
    ! Calculate transfer coefficient
    !
    xyr_VelTransCoef (:,:,0)    = 0.0_DP
    xyr_VelTransCoef (:,:,kmax) = 0.0_DP
    xyr_TempTransCoef(:,:,0)    = 0.0_DP
    xyr_TempTransCoef(:,:,kmax) = 0.0_DP
    xyr_QMixTransCoef(:,:,0)    = 0.0_DP
    xyr_QMixTransCoef(:,:,kmax) = 0.0_DP

    do k = 1, kmax-1
      xyr_VelTransCoef(:,:,k) =                                     &
        &   xyr_VelDiffCoef(:,:,k)                                  &
        &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
        &     / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )

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

      xyr_QMixTransCoef(:,:,k) =                                    &
        &   xyr_QMixDiffCoef(:,:,k)                                 &
        &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
        &     / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
    end do

    ! フラックスの計算
    ! Calculate fluxes
    !
    xyr_MomFluxX(:,:,0)    = 0.0_DP
    xyr_MomFluxX(:,:,kmax) = 0.0_DP
    xyr_MomFluxY(:,:,0)    = 0.0_DP
    xyr_MomFluxY(:,:,kmax) = 0.0_DP
    xyr_HeatFlux(:,:,0)    = 0.0_DP
    xyr_HeatFlux(:,:,kmax) = 0.0_DP
    do n = 1, ncmax
      xyrf_QMixFlux(:,:,0,n)    = 0.0_DP
      xyrf_QMixFlux(:,:,kmax,n) = 0.0_DP
    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


  end subroutine VDiffusionCalcFlux

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

  subroutine VDiffusionCalcFlux1D(                     &
    & z_U, z_V, zf_QMix,                               & ! (in)
    & z_Temp, r_VirTemp, r_Press,                      & ! (in)
    & z_Height, z_Exner, r_Exner,                      & ! (in)
    & r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef,   & ! (in)
    & r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux  & ! (out)
    & )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated. 
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & FKarm, &
                              ! $ k $ .
                              ! カルマン定数. 
                              ! Karman constant
      & 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

    ! 時刻管理
    ! 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):: z_U (1:kmax)
                              ! $ u $ .   東西風速. Eastward wind
    real(DP), intent(in):: z_V (1:kmax)
                              ! $ v $ .   南北風速. Northward wind
    real(DP), intent(in):: zf_QMix(1:kmax, 1:ncmax)
                              ! $ q $ .   質量混合比. Mass mixing ratio
    real(DP), intent(in):: z_Temp (1:kmax)
                              ! $ T $ .   温度. Temperature
    real(DP), intent(in):: r_VirTemp (0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual temperature (half level)
    real(DP), intent(in):: r_Press  (0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: z_Height (1:kmax)
                              ! 高度 (整数レベル). 
                              ! Height (full level)
    real(DP), intent(in):: z_Exner (1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: r_Exner (0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in):: r_VelDiffCoef (0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(in):: r_TempDiffCoef (0:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: r_QMixDiffCoef (0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity

    real(DP), intent(out):: r_MomFluxX (0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward momentum flux
    real(DP), intent(out):: r_MomFluxY (0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP), intent(out):: r_HeatFlux (0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(out):: rf_QMixFlux(0:kmax, 1:ncmax)
                              ! 質量フラックス. 
                              ! Mass flux of compositions

    ! 作業変数
    ! Work variables
    !
    real(DP) :: r_VelTransCoef (0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP) :: r_TempTransCoef(0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP) :: r_QMixTransCoef(0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass (composition)

    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_my_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! 輸送係数の計算
    ! Calculate transfer coefficient
    !
    r_VelTransCoef (0)    = 0.0_DP
    r_VelTransCoef (kmax) = 0.0_DP
    r_TempTransCoef(0)    = 0.0_DP
    r_TempTransCoef(kmax) = 0.0_DP
    r_QMixTransCoef(0)    = 0.0_DP
    r_QMixTransCoef(kmax) = 0.0_DP

    do k = 1, kmax-1
      r_VelTransCoef(k) =                                                   &
        &                 r_VelDiffCoef(k)                                  &
        &                   * r_Press(k) / ( GasRDry * r_VirTemp(k) ) &
        &                   / ( z_Height(k+1) - z_Height(k) )

      r_TempTransCoef(k) =                                                  &
        &                 r_TempDiffCoef(k)                                 &
        &                   * r_Press(k) / ( GasRDry * r_VirTemp(k) ) &
        &                   / ( z_Height(k+1) - z_Height(k) )

      r_QMixTransCoef(k) =                                                  &
        &                 r_QMixDiffCoef(k)                                 &
        &                   * r_Press(k) / ( GasRDry * r_VirTemp(k) ) &
        &                   / ( z_Height(k+1) - z_Height(k) )
    end do

    ! フラックスの計算
    ! Calculate fluxes
    !
    r_MomFluxX(0)    = 0.0_DP
    r_MomFluxX(kmax) = 0.0_DP
    r_MomFluxY(0)    = 0.0_DP
    r_MomFluxY(kmax) = 0.0_DP
    r_HeatFlux(0)    = 0.0_DP
    r_HeatFlux(kmax) = 0.0_DP
    do n = 1, ncmax
      rf_QMixFlux(0,n)    = 0.0_DP
      rf_QMixFlux(kmax,n) = 0.0_DP
    end do

    do k = 1, kmax-1
      r_MomFluxX(k) =   &
        & - r_VelTransCoef(k) * ( z_U(k+1) - z_U(k) )

      r_MomFluxY(k) =   &
        & - r_VelTransCoef(k) * ( z_V(k+1) - z_V(k) )

      r_HeatFlux(k) = &
        & - CpDry * r_TempTransCoef(k) * r_Exner(k) &
        &   * (   z_Temp(k+1) / z_Exner(k+1) &
        &       - z_Temp(k)   / z_Exner(k)     )
    end do

    do n = 1, ncmax
      do k = 1, kmax-1
        rf_QMixFlux(k,n) = &
          & - r_QMixTransCoef(k) &
          &   * ( zf_QMix(k+1,n) - zf_QMix(k,n)  )
      end do
    end do


  end subroutine VDiffusionCalcFlux1D

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

  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
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

    ! 物理定数設定
    ! 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_Sh (0:imax-1, 1:jmax, 0:kmax)
                              ! $ S_h $ (温度, 比湿). 
                              ! $ S_h $ (temperature, specific humidity)
    real(DP):: xyr_Sm (0:imax-1, 1:jmax, 0:kmax)
                              ! $ S_m $ (運動量). 
                              ! $ 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

    real(DP):: xyr_TurKinEne(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Turbulent kinetic energy

    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_my_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! 定数計算
    ! Calculate constants
    !
    Gamma1 = ( 1.0_DP / 3.0_DP ) - ( 2.0_DP * MYConstA1 / MYConstB1 )
    Gamma2 =   ( MYConstB2 / MYConstB1 ) &
      &      + ( 6.0_DP * MYConstA1 / MYConstB1 )
    Alpha1 = 3.0_DP  * MYConstA2 * Gamma1
    Alpha2 = 3.0_DP  * MYConstA2 * ( Gamma1 + Gamma2 )
    Beta1  = MYConstA1 * MYConstB1 * ( Gamma1 - MYConstC1 )
    Beta2  = MYConstA1 &
      &      * (   MYConstB1 * ( Gamma1 - MYConstC1 ) &
      &          + 6.0_DP * MYConstA1 &
      &          + 3.0_DP * MYConstA2 )
    Beta3  = MYConstA2 * MYConstB1 * Gamma1
    Beta4  = MYConstA2 * (   MYConstB1 * ( Gamma1 + Gamma2 ) &
      &                       - 3.0_DP * MYConstA1 )
    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_Sh(:,:,kmax) = 0.0_DP
    xyr_Sm(:,:,kmax) = 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_Sh(i,j,k) = &
              &   (   Alpha1 - Alpha2 * xyr_FluxRiNum(i,j,k) ) &
              &   / ( 1.0_DP - 1.0_DP * xyr_FluxRiNum(i,j,k) )

            xyr_Sm(i,j,k) = &
              &   (   Beta1 - Beta2 * xyr_FluxRiNum(i,j,k) ) &
              &   / ( Beta3 - Beta4 * xyr_FluxRiNum(i,j,k) ) &
              &   * xyr_Sh(i,j,k)

            xyr_Sh(i,j,k) = max( xyr_Sh(i,j,k), ShMin )
            xyr_Sm(i,j,k) = max( xyr_Sm(i,j,k), SmMin )

          else

            xyr_Sh(i,j,k) = ShMin
            xyr_Sm(i,j,k) = SmMin

          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.0_DP + FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) / MixLengthMax )
    end do

    ! 拡散係数の算出
    ! Calculate diffusion constants
    !
    xyr_VelDiffCoef = &
      &   xyr_MixLength**2 &
      &   * xyr_DVelDz &
      &   * sqrt ( MYConstB1 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_Sm ) &
      &   * xyr_Sm

    xyr_TempDiffCoef = &
      &   xyr_MixLength ** 2 &
      &   * xyr_DVelDz &
      &   * sqrt ( MYConstB1 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_Sm ) &
      &   * xyr_Sh

    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 &
            &    )
        end do
      end do
    end do

    xyr_QMixDiffCoef = xyr_TempDiffCoef

    xyr_VelDiffCoef (:,:,0)    = 0.0_DP
    xyr_VelDiffCoef (:,:,kmax) = 0.0_DP
    xyr_TempDiffCoef(:,:,0)    = 0.0_DP
    xyr_TempDiffCoef(:,:,kmax) = 0.0_DP
    xyr_QMixDiffCoef(:,:,0)    = 0.0_DP
    xyr_QMixDiffCoef(:,:,kmax) = 0.0_DP


    ! Calculation of turbulent kinetic energy
    ! Turbulent kinetic energy is diagnosed.
    xyr_TurKinEne = MYConstB1 * xyr_MixLength**2              &
      & * ( 1.0_DP - xyr_FluxRiNum ) * xyr_Sm * xyr_DVelDz**2 &
      & / 2.0_DP
    xyr_TurKinEne(:,:,0   ) = 0.0_DP
    xyr_TurKinEne(:,:,kmax) = 0.0_DP

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'TurKinEne', xyr_TurKinEne )


  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_my_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)
    & xyr_Press, xyz_Exner, xyr_Exner,                         & ! (in)
    & xyr_VirTemp, xyz_Height,                                 & ! (in)
    & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef      & ! (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: &
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure
      & LatentHeat, &
                              ! $ L $ [J kg-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation
      & GasRDry
                              ! $ R $ [J kg-1 K-1]. 
                              ! 乾燥大気の気体定数. 
                              ! Gas constant of air

    ! 時刻管理
    ! 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):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (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(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual temperature (half level)
    real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
                              ! 高度 (整数レベル). 
                              ! Height (full level)

    real(DP), intent(in):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(in):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity

    ! 出力のための作業変数
    ! 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
    !
    real(DP) :: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP) :: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP) :: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass of constituents

    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_my_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 transfer coefficient
    !
    xyr_VelTransCoef (:,:,0)    = 0.0_DP
    xyr_VelTransCoef (:,:,kmax) = 0.0_DP
    xyr_TempTransCoef(:,:,0)    = 0.0_DP
    xyr_TempTransCoef(:,:,kmax) = 0.0_DP
    xyr_QMixTransCoef(:,:,0)    = 0.0_DP
    xyr_QMixTransCoef(:,:,kmax) = 0.0_DP

    do k = 1, kmax-1
      xyr_VelTransCoef(:,:,k) =                                     &
        &   xyr_VelDiffCoef(:,:,k)                                  &
        &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
        &     / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )

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

      xyr_QMixTransCoef(:,:,k) =                                    &
        &   xyr_QMixDiffCoef(:,:,k)                                 &
        &     * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) &
        &     / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
    end do


    ! 風速, 温度, 比湿フラックス補正
    ! 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.0_DP
    xyr_MomFluxXCor   (:,:,kmax) = 0.0_DP
    xyr_MomFluxYCor   (:,:,0)    = 0.0_DP
    xyr_MomFluxYCor   (:,:,kmax) = 0.0_DP
    xyr_HeatFluxCor   (:,:,0)    = 0.0_DP
    xyr_HeatFluxCor   (:,:,kmax) = 0.0_DP
    do n = 1, ncmax
      xyrf_QMixFluxCor(:,:,0,n)    = 0.0_DP
      xyrf_QMixFluxCor(:,:,kmax,n) = 0.0_DP
    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_my モジュールの初期化を行います. 
    ! NAMELIST#vdiffusion_my_nml の読み込みはこの手続きで行われます. 
    !
    ! "vdiffusion_my" module is initialized. 
    ! "NAMELIST#vdiffusion_my_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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & AxnameX, &
      & AxnameY, &
      & AxnameZ, &
      & AxnameR, &
      & AxnameT

    ! メッセージ制御
    ! Message control
    !
    use mpi_messagecntl, only : DoesOutputMPIMessage


    ! 宣言文 ; 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_my_nml/        &
      & FlagConstDiffCoef,              &
      & ConstDiffCoefM, ConstDiffCoefH, &
!
      & SquareVelMin, BulkRiNumMin,     &
!
      & MixLengthMax, ShMin, SmMin,      &
      & VelDiffCoefMin, TempDiffCoefMin, &
      & VelDiffCoefMax, TempDiffCoefMax, &
!
      & MYConstA1, MYConstB1, MYConstA2, MYConstB2, MYConstC1
          !
          ! デフォルト値については初期化手続 "vdiffusion_my#VDiffInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "vdiffusion_my#VDiffInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( vdiffusion_my_inited ) return


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

    SquareVelMin    =     0.1_DP
    BulkRiNumMin    = - 100.0_DP

    MixLengthMax    = 300.0_DP
    ShMin           =   0.0_DP
    SmMin           =   0.0_DP
    VelDiffCoefMin  =   0.1_DP
    TempDiffCoefMin =   0.1_DP
    VelDiffCoefMax  = 10000.0_DP
    TempDiffCoefMax = 10000.0_DP

    ! Parameters proposed by Mellor and Yamada (1982).
    !
    MYConstA1 =  0.92_DP
    MYConstB1 = 16.6_DP
    MYConstA2 =  0.74_DP
    MYConstB2 = 10.1_DP
    MYConstC1 =  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_my_nml,     &  ! (out)
        & iostat = iostat_nml )           ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
      if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = vdiffusion_my_nml )
    end if


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

    call HistoryAutoAddVariable( 'MomFluxX', &
      & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
      & 'eastward momentum flux', 'N m-2' )
    call HistoryAutoAddVariable( 'MomFluxY', &
      & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
      & 'northward momentum flux', 'N m-2' )
    call HistoryAutoAddVariable( 'HeatFlux', &
      & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
      & 'heat flux', 'W m-2' )
    call HistoryAutoAddVariable( 'QVapFlux', &
      & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
      & 'moisture flux', 'W m-2' )

    call HistoryAutoAddVariable( 'DUDtVDiff', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'tendency of zonal wind by vertical diffusion', 'm s-2' )
    call HistoryAutoAddVariable( 'DVDtVDiff', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'tendency of meridional wind by vertical diffusion', 'm s-2' )
    call HistoryAutoAddVariable( 'DTempDtVDiff', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'tendency of temperature by vertical diffusion', 'K s-1' )
    call HistoryAutoAddVariable( 'DQVapDtVDiff', &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'tendency of specific humidity by vertical diffusion', 's-1' )

    call HistoryAutoAddVariable( 'TurKinEne', &
      & (/ AxNameX, AxNameY, AxNameR, AxNameT /), &
      & 'turbulent kinetic energy', 'm2 s-2' )

    call HistoryAutoAddVariable( 'TKEPShear',  &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'turbulent kinetic energy production rate by shear', 'm2 s-3' )
    call HistoryAutoAddVariable( 'TKEPBuoy',   &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'turbulent kinetic energy production rate by buoyancy', 'm2 s-3' )
    call HistoryAutoAddVariable( 'TKEDiss',    &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'turbulent kinetic energy dissipation rate', 'm2 s-3' )
    call HistoryAutoAddVariable( 'MixLength',  &
      & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
      & 'mixing length', 'm' )

    ! 印字 ; 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, '  ShMin             = %f', d = (/ ShMin       /) )
    call MessageNotify( 'M', module_name, '  SmMin             = %f', d = (/ SmMin       /) )
    call MessageNotify( 'M', module_name, '  VelDiffCoefMin    = %f', d = (/ VelDiffCoefMin  /) )
    call MessageNotify( 'M', module_name, '  TempDiffCoefMin   = %f', d = (/ TempDiffCoefMin /) )
    call MessageNotify( 'M', module_name, '  VelDiffCoefMax    = %f', d = (/ VelDiffCoefMax  /) )
    call MessageNotify( 'M', module_name, '  TempDiffCoefMax   = %f', d = (/ TempDiffCoefMax /) )
    call MessageNotify( 'M', module_name, '  MYConstA1         = %f', d = (/ MYConstA1     /) )
    call MessageNotify( 'M', module_name, '  MYConstB1         = %f', d = (/ MYConstB1     /) )
    call MessageNotify( 'M', module_name, '  MYConstA2         = %f', d = (/ MYConstA2     /) )
    call MessageNotify( 'M', module_name, '  MYConstB2         = %f', d = (/ MYConstB2     /) )
    call MessageNotify( 'M', module_name, '  MYConstC1         = %f', d = (/ MYConstC1     /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    vdiffusion_my_inited = .true.

  end subroutine VDiffusionInit

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

end module vdiffusion_my
