Class vdiffusion_my
In: vdiffusion/vdiffusion_my.f90

鉛直拡散フラックス (Mellor and Yamada, 1974, 1982)

Vertical diffusion flux (Mellor and Yamada, 1974, 1982)

Note that Japanese and English are described in parallel.

鉛直拡散フラックスを計算します.

Vertical diffusion flux is calculated.

Procedures List

VDiffusion :鉛直拡散フラックスの計算
VDiffusionOutPut :フラックスの出力
———— :————
VDiffusion :Calculate vertical diffusion fluxes
VDiffusionOutPut :Output fluxes

NAMELIST

NAMELIST#vdiffusion_my_nml

Methods

BulkRiNumMin   ConstDiffCoefH   ConstDiffCoefM   FlagConstDiffCoef   MYConstA1   MYConstA2   MYConstB1   MYConstB2   MYConstC1   MixLengthMax   ShMin   SmMin   SquareVelMin   TempDiffCoefMax   TempDiffCoefMin   VDiffCoefficient   VDiffusion   VDiffusionExpTendency   VDiffusionInit   VDiffusionMY25   VDiffusionMY25Testing   VDiffusionOutPut   VDiffusionTransCoefAndFlux   VelDiffCoefMax   VelDiffCoefMin   module_name   vdiffusion_my_inited   version  

Included Modules

gridset composition dc_types dc_message constants timeset gtool_historyauto mpi_wrapper phy_implicit_utils namelist_util dc_iounit dc_string axesset

Public Instance methods

VDiffusion( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
Subroutine :
xyz_U(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ u $ . 東西風速. Eastward wind
xyz_V(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ v $ . 南北風速. Northward wind
xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) :real(DP), intent(in)
: $ q $ . 質量混合比. Mass mixing ratio
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ T $ . 温度. Temperature
xyr_Temp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{T} $ . 温度 (半整数レベル). Temperature (half level)
xyz_VirTemp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ T_v $ . 仮温度. Virtual temperature
xyr_VirTemp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{T}_v $ . 仮温度 (半整数レベル). Virtual temperature (half level)
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half level)
xy_SurfHeight(0:imax-1,1:jmax) :real(DP), intent(in)
: $ z_s $ . 地表面高度. Surface height.
xyz_Height(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: 高度 (整数レベル). Height (full level)
xyr_Height(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 高度 (半整数レベル). Height (half level)
xyz_Exner(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: Exner 関数 (整数レベル). Exner function (full level)
xyr_Exner(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: Exner 関数 (半整数レベル). Exner function (half level)
xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 東西方向運動量フラックス. Eastward momentum flux
xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 南北方向運動量フラックス. Northward momentum flux
xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 熱フラックス. Heat flux
xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) :real(DP), intent(out)
: 質量フラックス. Mass flux of compositions
xyr_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：運動量. Transfer coefficient: velocity
xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：温度. Transfer coefficient: temperature
xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：質量. Transfer coefficient: mass (composition)

鉛直拡散フラックスを計算します.

Vertical diffusion flux is calculated by use of MY2 model.

[Source]

  subroutine VDiffusion( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! Vertical diffusion flux is calculated by use of MY2 model.
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: FKarm, Grav
                              ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration

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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u $ .   東西風速. Eastward wind
    real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v $ .   南北風速. Northward wind
    real(DP), intent(in):: 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_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP), intent(out):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(out):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass (composition)

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

    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.0d0
      xyr_VelDiffCoef (:,:,1:kmax-1) = ConstDiffCoefM
      xyr_VelDiffCoef (:,:,kmax    ) = 0.0d0

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

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

    else

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

      do k = 1, kmax-1
        xyr_DVelDz(:,:,k) = sqrt( max( SquareVelMin , ( xyz_U(:,:,k+1) - xyz_U(:,:,k) )**2 + ( xyz_V(:,:,k+1) - xyz_V(:,:,k) )**2 ) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )

        xyr_BulkRiNum(:,:,k) = Grav / ( xyr_VirTemp(:,:,k) / xyr_Exner(:,:,k) ) * (   xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) - xyz_VirTemp(:,:,k  ) / xyz_Exner(:,:,k  ) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) / xyr_DVelDz(:,:,k)**2

        xyr_BulkRiNum(:,:,k) = max( xyr_BulkRiNum(:,:,k) , BulkRiNumMin )
      end do

      ! 拡散係数の計算
      ! Calculate diffusion coefficients
      !
      call VDiffCoefficient( xy_SurfHeight, xyr_Height, xyr_DVelDz, xyr_BulkRiNum, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef )

    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 VDiffusionTransCoefAndFlux( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_VirTemp, xyr_Press, xyz_Height, xyz_Exner, xyr_Exner, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )

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

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

時間変化率の計算を行います.

Calculate tendencies.

[Source]

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, 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
VDiffusionInit( )
Subroutine :

vdiffusion_my モジュールの初期化を行います. NAMELIST#vdiffusion_my_nml の読み込みはこの手続きで行われます.

"vdiffusion_my" module is initialized. "NAMELIST#vdiffusion_my_nml" is loaded in this procedure.

This procedure input/output NAMELIST#vdiffusion_my_nml .

[Source]

  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

    ! 宣言文 ; 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.0d0
    ConstDiffCoefH    = 0.0d0

    SquareVelMin    =     0.1
    BulkRiNumMin    = - 100.

    MixLengthMax    = 300.
    ShMin           =   0.
    SmMin           =   0.
    VelDiffCoefMin  =   0.1
    TempDiffCoefMin =   0.1
    VelDiffCoefMax  = 10000.
    TempDiffCoefMax = 10000.

    ! Parameters proposed by Mellor and Yamada (1982).
    !
    MYConstA1 =  0.92
    MYConstB1 = 16.6
    MYConstA2 =  0.74
    MYConstB2 = 10.1
    MYConstC1 =  0.08


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

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

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
      if ( iostat_nml == 0 ) 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
VDiffusionMY25( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyz_TurKinEne, xy_SurfMomFluxX, xy_SurfMomFluxY, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, xyz_DTurKinEneDt )
Subroutine :
xyz_U(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ u $ . 東西風速. Eastward wind
xyz_V(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ v $ . 南北風速. Northward wind
xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) :real(DP), intent(in)
: $ q $ . 質量混合比. Mass mixing ratio
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ T $ . 温度. Temperature
xyr_Temp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{T} $ . 温度 (半整数レベル). Temperature (half level)
xyz_VirTemp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ T_v $ . 仮温度. Virtual temperature
xyr_VirTemp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{T}_v $ . 仮温度 (半整数レベル). Virtual temperature (half level)
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half level)
xy_SurfHeight(0:imax-1,1:jmax) :real(DP), intent(in)
: $ z_s $ . 地表面高度. Surface height.
xyz_Height(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: 高度 (整数レベル). Height (full level)
xyr_Height(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 高度 (半整数レベル). Height (half level)
xyz_Exner(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: Exner 関数 (整数レベル). Exner function (full level)
xyr_Exner(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: Exner 関数 (半整数レベル). Exner function (half level)
xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: Turbulent kinetic energy (m2 s-2)
xy_SurfMomFluxX(0:imax-1, 1:jmax) :real(DP), intent(in)
: Eastward momentum flux at surface
xy_SurfMomFluxY(0:imax-1, 1:jmax) :real(DP), intent(in)
: Northward momentum flux at surface
xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 東西方向運動量フラックス. Eastward momentum flux
xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 南北方向運動量フラックス. Northward momentum flux
xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 熱フラックス. Heat flux
xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) :real(DP), intent(out)
: 質量フラックス. Mass flux of compositions
xyr_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：運動量. Transfer coefficient: velocity
xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：温度. Transfer coefficient: temperature
xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：質量. Transfer coefficient: mass (composition)
xyz_DTurKinEneDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
: Tendency of turbulent kinetic energy

鉛直拡散フラックスを計算します.

Vertical diffusion flux is calculated by use of MY2.5 model.

[Source]

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: FKarm, Grav, GasRDry, CpDry
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, 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_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP), intent(out):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(out):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass (composition)
    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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP) :: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP) :: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity
    real(DP) :: 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.0d-3 )**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 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

    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
    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 )


!!$    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 VDiffusionTransCoefAndFlux( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_VirTemp, xyr_Press, xyz_Height, xyz_Exner, xyr_Exner, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )


    ! 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 / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneNonZero**1.5
    xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneNonZero**0.5


    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. * 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. * 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. * 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, imax * jmax, kmax )

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

      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, a_FlagReCalcGlobal )
      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
VDiffusionOutPut( xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
Subroutine :
xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 東西方向運動量フラックス. Eastward wind flux
xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 南北方向運動量フラックス. Northward momentum flux
xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 熱フラックス. Heat flux
xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) :real(DP), intent(in)
: 質量フラックス. Mass flux of constituents
xyz_DUDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ DP{u}{t} $ . 東西風速変化. Eastward wind tendency
xyz_DVDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ DP{v}{t} $ . 南北風速変化. Northward wind tendency
xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ DP{T}{t} $ . 温度変化. Temperature tendency
xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) :real(DP), intent(in)
: $ DP{q}{t} $ . 混合比変化. Mass mixing ratio tendency
xyz_Exner(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: Exner 関数 (整数レベル). Exner function (full level)
xyr_Exner(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: Exner 関数 (半整数レベル). Exner function (half level)
xyr_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 輸送係数：運動量. Transfer coefficient: velocity
xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 輸送係数：温度. Transfer coefficient: temperature
xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 輸送係数：質量. Transfer coefficient: mass of constituents

フラックス (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.

[Source]

  subroutine VDiffusionOutPut( xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyz_DUDt,  xyz_DVDt,  xyz_DTempDt,  xyzf_DQMixDt, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
    !
    ! フラックス (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, LatentHeat
                              ! $ L $ [J kg-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation

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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward wind flux
    real(DP), intent(in):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 質量フラックス. 
                              ! Mass flux of constituents
    real(DP), intent(in):: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{u}{t} $ . 東西風速変化. 
                              ! Eastward wind tendency
    real(DP), intent(in):: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{v}{t} $ . 南北風速変化. 
                              ! Northward wind tendency
    real(DP), intent(in):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化. 
                              ! Temperature tendency
    real(DP), intent(in):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ \DP{q}{t} $ . 混合比変化. 
                              ! Mass mixing ratio tendency
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)
    real(DP), intent(in):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP), intent(in):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(in):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass of constituents

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


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

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. vdiffusion_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 )


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

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

          xyr_MomFluxXCor( i,j,k ) = xyr_MomFluxX( i,j,k ) + ( xyz_DUDt( i,j,k ) - xyz_DUDt( i,j,k+1 ) ) * xyr_VelTransCoef( i,j,k ) * DelTime

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

          xyr_HeatFluxCor( i,j,k ) = xyr_HeatFlux( i,j,k ) + (   xyz_DTempDt( i,j,k   ) / xyz_Exner( i,j,k   ) - xyz_DTempDt( i,j,k+1 ) / xyz_Exner( i,j,k+1 ) ) * CpDry * xyr_TempTransCoef( i,j,k ) * xyr_Exner( i,j,k ) * DelTime
        end do
      end do
    end do

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

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

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

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

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

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

  end subroutine VDiffusionOutPut

Private Instance methods

BulkRiNumMin
Variable :
BulkRiNumMin :real(DP), save
: バルク $ R_i $ 数最小値. Minimum value of bulk $ R_i $
ConstDiffCoefH
Variable :
ConstDiffCoefH :real(DP), save
: Diffusion coefficient for heat that is used in the case of FlagConstDiffCoef == .true.
ConstDiffCoefM
Variable :
ConstDiffCoefM :real(DP), save
: Diffusion coefficient for momentum that is used in the case of FlagConstDiffCoef == .true.
FlagConstDiffCoef
Variable :
FlagConstDiffCoef :logical , save
: Flag for use of constant diffusion coefficient
MYConstA1
Variable :
MYConstA1 :real(DP), save
MYConstA2
Variable :
MYConstA2 :real(DP), save
MYConstB1
Variable :
MYConstB1 :real(DP), save
MYConstB2
Variable :
MYConstB2 :real(DP), save
MYConstC1
Variable :
MYConstC1 :real(DP), save
MixLengthMax
Variable :
MixLengthMax :real(DP), save
: 最大混合距離. Maximum mixing length
ShMin
Variable :
ShMin :real(DP), save
: $ S_h $ 最小値. Minimum $ S_h $
SmMin
Variable :
SmMin :real(DP), save
: $ S_m $ 最小値. Minimum $ S_m $
SquareVelMin
Variable :
SquareVelMin :real(DP), save
: 風二乗差最小値. Minimum value of square of velocity
TempDiffCoefMax
Variable :
TempDiffCoefMax :real(DP), save
: $ T $ 拡散係数最大値. Maximum diffusion coefficient of $ T $
TempDiffCoefMin
Variable :
TempDiffCoefMin :real(DP), save
: $ T $ 拡散係数最小値. Minimum diffusion coefficient of $ T $
VDiffCoefficient( xy_SurfHeight, xyr_Height, xyr_DVelDz, xyr_BulkRiNum, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef )
Subroutine :
xy_SurfHeight(0:imax-1,1:jmax) :real(DP), intent(in)
: $ z_s $ . 地表面高度. Surface height.
xyr_Height(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 高度 (半整数レベル). Height (half level)
xyr_DVelDz(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ DD{|Dvect{v}|}{z} $
xyr_BulkRiNum(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: バルク $ R_i $ 数. Bulk $ R_i $
xyr_VelDiffCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 拡散係数：運動量. Diffusion coefficient: velocity
xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 拡散係数：温度. Transfer coefficient: temperature
xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 拡散係数：質量 Diffusion coefficient: mass of constituents

鉛直拡散フラックスを計算します.

Vertical diffusion flux is calculated.

[Source]

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

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimeN, 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. / 3. ) - ( 2. * MYConstA1 / MYConstB1 )
    Gamma2 =   ( MYConstB2 / MYConstB1 ) + ( 6. * MYConstA1 / MYConstB1 )
    Alpha1 = 3.  * MYConstA2 * Gamma1
    Alpha2 = 3.  * MYConstA2 * ( Gamma1 + Gamma2 )
    Beta1  = MYConstA1 * MYConstB1 * ( Gamma1 - MYConstC1 )
    Beta2  = MYConstA1 * (   MYConstB1 * ( Gamma1 - MYConstC1 ) + 6. * MYConstA1 + 3. * MYConstA2 )
    Beta3  = MYConstA2 * MYConstB1 * Gamma1
    Beta4  = MYConstA2 * (   MYConstB1 * ( Gamma1 + Gamma2 ) - 3. * 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. * Beta2 * Beta3 * xyr_BulkRiNum ) ) / ( 2. * Beta2 )

    ! $ \tilde{S_h} $ と $ \tilde{S_m} $ の算出
    ! Calculate $ \tilde{S_h} $ and $ \tilde{S_m} $
    !
    xyr_Sh(:,:,kmax) = 0.
    xyr_Sm(:,:,kmax) = 0.

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

          if ( xyr_FluxRiNum(i,j,k) < CrtlFluxRiNum ) then 

            xyr_Sh(i,j,k) = (   Alpha1 - Alpha2 * xyr_FluxRiNum(i,j,k) ) / ( 1. - 1. * 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. + FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) / MixLengthMax )
    end do

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

    xyr_TempDiffCoef = xyr_MixLength ** 2 * xyr_DVelDz * sqrt ( MYConstB1 * ( 1. - 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.
    xyr_VelDiffCoef (:,:,kmax) = 0.
    xyr_TempDiffCoef(:,:,0)    = 0.
    xyr_TempDiffCoef(:,:,kmax) = 0.
    xyr_QMixDiffCoef(:,:,0)    = 0.
    xyr_QMixDiffCoef(:,:,kmax) = 0.


    ! 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
VDiffusionMY25Testing( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyz_TurKinEne, xy_SurfMomFluxX, xy_SurfMomFluxY, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, xyz_DTurKinEneDt )
Subroutine :
xyz_U(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ u $ . 東西風速. Eastward wind
xyz_V(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ v $ . 南北風速. Northward wind
xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) :real(DP), intent(in)
: $ q $ . 質量混合比. Mass mixing ratio
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ T $ . 温度. Temperature
xyr_Temp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{T} $ . 温度 (半整数レベル). Temperature (half level)
xyz_VirTemp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ T_v $ . 仮温度. Virtual temperature
xyr_VirTemp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{T}_v $ . 仮温度 (半整数レベル). Virtual temperature (half level)
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half level)
xy_SurfHeight(0:imax-1,1:jmax) :real(DP), intent(in)
: $ z_s $ . 地表面高度. Surface height.
xyz_Height(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: 高度 (整数レベル). Height (full level)
xyr_Height(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 高度 (半整数レベル). Height (half level)
xyz_Exner(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: Exner 関数 (整数レベル). Exner function (full level)
xyr_Exner(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: Exner 関数 (半整数レベル). Exner function (half level)
xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: Turbulent kinetic energy (m2 s-2)
xy_SurfMomFluxX(0:imax-1, 1:jmax) :real(DP), intent(in)
: Eastward momentum flux at surface
xy_SurfMomFluxY(0:imax-1, 1:jmax) :real(DP), intent(in)
: Northward momentum flux at surface
xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 東西方向運動量フラックス. Eastward momentum flux
xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 南北方向運動量フラックス. Northward momentum flux
xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 熱フラックス. Heat flux
xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) :real(DP), intent(out)
: 質量フラックス. Mass flux of compositions
xyr_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：運動量. Transfer coefficient: velocity
xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：温度. Transfer coefficient: temperature
xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：質量. Transfer coefficient: mass (composition)
xyz_DTurKinEneDt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)
: Tendency of turbulent kinetic energy

鉛直拡散フラックスを計算します.

Vertical diffusion flux is calculated by use of MY2.5 model.

[Source]

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: FKarm, Grav, GasRDry, CpDry
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, 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_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP), intent(out):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(out):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass (composition)
    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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP) :: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP) :: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity
    real(DP) :: 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_DTurKinEneDtVDiff(0:imax-1, 1:jmax, 1:kmax)


    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.0d-3 )**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


    real(DP) :: xyz_TurKinEneTentative(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: DelTimeSubStep


    ! 実行文 ; 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 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

    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
    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 )


!!$    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 VDiffusionTransCoefAndFlux( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_VirTemp, xyr_Press, xyz_Height, xyz_Exner, xyr_Exner, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )


    ! 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



    ! Calculation of vertical mixing

    !
    ! 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. * 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. * 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. * 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, imax * jmax, kmax )

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

    xyz_DTurKinEneDtVDiff = xyz_DelTurKinEneLUVec / ( 2.0_DP * DelTime )


    ! Calculation tendency including source and dissipation terms


    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 / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneNonZero**1.5
    xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneNonZero**0.5


    ! 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


    xyz_TurKinEneTentative = xyz_TurKinEne

    DelTimeSubStep = 0.25_DP
    nloop = 2.0_DP * DelTime / DelTimeSubStep + 1
    DelTimeSubStep = 2.0_DP * DelTime / dble( nloop )

    nloop = kmax
    loop_solve : do iloop = 1, nloop

!!$      xyz_TurKinEneProShear =                              &
!!$        &   sqrt( 2.0_DP ) * xyz_MixLength * xyz_DVelDzSq  &
!!$        & * xyz_Sm                                         &
!!$        & * sqrt( xyz_TurKinEneTentative )
!!$      xyz_TurKinEneProBuoya =                              &
!!$        & - sqrt( 2.0_DP ) * xyz_MixLength * xyz_StatStab  &
!!$        & * xyz_Sh                                         &
!!$        & * sqrt( xyz_TurKinEneTentative )
      xyz_TurKinEneDiss = 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) * xyz_TurKinEneTentative**1.5

!!$      xyz_TurKinEneTentative =        &
!!$        &   xyz_TurKinEneTentative    &
!!$        & + (   xyz_DTurKinEneDtVDiff &
!!$        &     + xyz_TurKinEneProShear &
!!$        &     + xyz_TurKinEneProBuoya &
!!$        &     - xyz_TurKinEneDiss )   &
!!$        & * DelTimeSubStep
      xyz_TurKinEneTentative = xyz_TurKinEneTentative + (   xyz_TurKinEneProShear + xyz_TurKinEneProBuoya - xyz_TurKinEneDiss ) * DelTimeSubStep

      xyz_TurKinEneTentative = max( xyz_TurKinEneTentative, 0.0_DP )

!!$      ! 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

    xyz_DTurKinEneDt = ( xyz_TurKinEneTentative - xyz_TurKinEne ) / ( 2.0_DP * DelTime )





    !=========================================================================





    ! Calculate turbulent kinetic energy with offset
    !
    xyz_TurKinEneNonZero = xyz_TurKinEne + TurKinEneOffset
!!$    xyz_TurKinEneNonZero = xyz_TurKinEneTentative + TurKinEneOffset


    ! 混合距離の算出
    ! 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 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

    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
    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 )


!!$    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 VDiffusionTransCoefAndFlux( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_VirTemp, xyr_Press, xyz_Height, xyz_Exner, xyr_Exner, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )


    ! 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



    ! Calculation of vertical mixing

    !
    ! 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. * 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. * 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. * 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) ) - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav * xyz_DTurKinEneDt(:,:,k)
    end do


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

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

    xyz_DTurKinEneDt = xyz_DelTurKinEneLUVec / ( 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_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 VDiffusionMY25Testing
VDiffusionTransCoefAndFlux( xyz_U, xyz_V, xyzf_QMix, xyz_Temp, xyr_VirTemp, xyr_Press, xyz_Height, xyz_Exner, xyr_Exner, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )
Subroutine :
xyz_U(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ u $ . 東西風速. Eastward wind
xyz_V(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ v $ . 南北風速. Northward wind
xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) :real(DP), intent(in)
: $ q $ . 質量混合比. Mass mixing ratio
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: $ T $ . 温度. Temperature
xyr_VirTemp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{T}_v $ . 仮温度 (半整数レベル). Virtual temperature (half level)
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half level)
xyz_Height(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: 高度 (整数レベル). Height (full level)
xyz_Exner(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: Exner 関数 (整数レベル). Exner function (full level)
xyr_Exner(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: Exner 関数 (半整数レベル). Exner function (half level)
xyr_VelDiffCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 拡散係数：運動量. Diffusion coefficient: velocity
xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 拡散係数：温度. Transfer coefficient: temperature
xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: 拡散係数：比湿. Diffusion coefficient: specific humidity
xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 東西方向運動量フラックス. Eastward momentum flux
xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 南北方向運動量フラックス. Northward momentum flux
xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 熱フラックス. Heat flux
xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) :real(DP), intent(out)
: 質量フラックス. Mass flux of compositions
xyr_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：運動量. Transfer coefficient: velocity
xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：温度. Transfer coefficient: temperature
xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 輸送係数：質量. Transfer coefficient: mass (composition)

鉛直拡散フラックスを計算します.

Vertical diffusion flux is calculated.

[Source]

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: FKarm, GasRDry, CpDry
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure

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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u $ .   東西風速. Eastward wind
    real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v $ .   南北風速. Northward wind
    real(DP), intent(in):: 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
    real(DP), intent(out):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP), intent(out):: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(out):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
                              ! 輸送係数：質量. 
                              ! Transfer coefficient: mass (composition)

    ! 作業変数
    ! Work variables
    !
    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.
    xyr_VelTransCoef (:,:,kmax) = 0.
    xyr_TempTransCoef(:,:,0)    = 0.
    xyr_TempTransCoef(:,:,kmax) = 0.
    xyr_QMixTransCoef(:,:,0)    = 0.
    xyr_QMixTransCoef(:,:,kmax) = 0.

    do k = 1, kmax-1
!!$      xyr_VelTransCoef(:,:,k) = &
!!$        &                 xyr_VelDiffCoef(:,:,k) &
!!$        &                   * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) &
!!$        &                   / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
!!$
!!$      xyr_TempTransCoef(:,:,k) = &
!!$        &                 xyr_TempDiffCoef(:,:,k) &
!!$        &                   * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) &
!!$        &                   / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
!!$
!!$      xyr_QMixTransCoef(:,:,k) = &
!!$        &                 xyr_QMixDiffCoef(:,:,k) &
!!$        &                   * xyr_Press(:,:,k) / GasRDry / xyr_Temp(:,:,k) &
!!$        &                   / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )


      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.
    xyr_MomFluxX(:,:,kmax) = 0.
    xyr_MomFluxY(:,:,0)    = 0.
    xyr_MomFluxY(:,:,kmax) = 0.
    xyr_HeatFlux(:,:,0)    = 0.
    xyr_HeatFlux(:,:,kmax) = 0.
    do n = 1, ncmax
      xyrf_QMixFlux(:,:,0,n)    = 0.
      xyrf_QMixFlux(:,:,kmax,n) = 0.
    end do

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

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

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

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


  end subroutine VDiffusionTransCoefAndFlux
VelDiffCoefMax
Variable :
VelDiffCoefMax :real(DP), save
: $ Dvect{u} $ 拡散係数最大値. Maximum diffusion coefficient of $ Dvect{u} $
VelDiffCoefMin
Variable :
VelDiffCoefMin :real(DP), save
: $ Dvect{u} $ 拡散係数最小値. Minimum diffusion coefficient of $ Dvect{u} $
module_name
Constant :
module_name = ‘vdiffusion_my‘ :character(*), parameter
: モジュールの名称. Module name
vdiffusion_my_inited
Variable :
vdiffusion_my_inited = .false. :logical, save
: 初期設定フラグ. Initialization flag
version
Constant :
version = ’$Name: dcpam5-20130921 $’ // ’$Id: vdiffusion_my.f90,v 1.1 2013-09-16 12:01:19 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version