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

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

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_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: specific humidity

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

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_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! 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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数:運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数:温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数:比湿. 
                              ! Diffusion coefficient: specific humidity

    ! 作業変数
    ! Work variables
    !
    real(DP) :: xyr_DVelDz (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \DD{|\Dvect{v}|}{z} $
    real(DP) :: xyr_BulkRiNum (0:imax-1, 1:jmax, 0:kmax)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $

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

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

    ! 実行文 ; Executable statement
    !

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


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


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

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

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

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

    else

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

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

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

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

      ! 拡散係数の計算
      ! Calculate diffusion coefficients
      !
      call VDiffCoefficient( xy_SurfHeight, 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 VDiffusionCalcFlux( 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 )

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

  end subroutine VDiffusion
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
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.0_DP
    ConstDiffCoefH    = 0.0_DP

    SquareVelMin    =     0.1_DP
    BulkRiNumMin    = - 100.0_DP

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

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


    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, 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
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_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)
: 拡散係数:温度. Diffusion coefficient: temperature
xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 拡散係数:比湿. Diffusion coefficient: specific humidity
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_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, 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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数:運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数:温度. 
                              ! Diffusion coefficient: temperature
    real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数:比湿. 
                              ! Diffusion coefficient: specific humidity

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

    ! 作業変数
    ! Work variables
    !

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    ! 実行文 ; Executable statement
    !

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


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


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

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

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

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


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


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

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


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


    ! 拡散係数の計算
    ! Calculation of diffusion coefficient
    !
    xyz_VelDiffCoef  = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sm
    xyz_TempDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sh
    !
    do k = 0, kmax
      if ( ( k == 0 ) .or. ( k == kmax ) ) then
        xyr_VelDiffCoef (:,:,k) = 0.0_DP
        xyr_TempDiffCoef(:,:,k) = 0.0_DP
      else
        xyr_VelDiffCoef (:,:,k) = ( xyz_VelDiffCoef (:,:,k) + xyz_VelDiffCoef (:,:,k+1) ) / 2.0_DP
        xyr_TempDiffCoef(:,:,k) = ( xyz_TempDiffCoef(:,:,k) + xyz_TempDiffCoef(:,:,k+1) ) / 2.0_DP
      end if
    end do
    !
    do k = 1, kmax-1
      do j = 1, jmax
        do i = 0, imax-1
          xyr_VelDiffCoef(i,j,k) = max( min( xyr_VelDiffCoef(i,j,k), VelDiffCoefMax ), VelDiffCoefMin )
          xyr_TempDiffCoef(i,j,k) = max( min( xyr_TempDiffCoef(i,j,k), TempDiffCoefMax ), TempDiffCoefMin )
        end do
      end do
    end do
    !
    xyr_QMixDiffCoef      = xyr_TempDiffCoef


    ! 輸送係数とフラックスの計算
    ! Calculate transfer coefficient and flux
    !
    call VDiffusionCalcFlux( xyz_U, xyz_V, xyzf_QMix, 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 )


    ! Calculate tendency of turbulent kinetic energy

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

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

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


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


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


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


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


    nloop = kmax
    loop_solve : do iloop = 1, nloop

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

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


      !
      ! Solve simultaneous linear equations by use of LU decomposition technique
      !
      xyza_TurKinEneLUMtx = xyza_TurKinEneMtx
      !
      call PhyImplLUDecomp3( xyza_TurKinEneLUMtx, 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
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_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)
: 拡散係数:温度. Diffusion coefficient: temperature
xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
: 拡散係数:比湿. Diffusion coefficient: specific humidity
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 VDiffusionMY251DWrapper3D( 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_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, 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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u $ .   東西風速. Eastward wind
    real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v $ .   南北風速. Northward wind
    real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q $ .   質量混合比. Mass mixing ratio
    real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .   温度. Temperature
    real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T} $ . 温度 (半整数レベル). 
                              ! Temperature (half level)
    real(DP), intent(in):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T_v $ .   仮温度. Virtual temperature
    real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual temperature (half level)
    real(DP), intent(in):: xyr_Press  (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xy_SurfHeight (0:imax-1,1:jmax)
                              ! $ z_s $ . 地表面高度. 
                              ! Surface height. 
    real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
                              ! 高度 (整数レベル). 
                              ! Height (full level)
    real(DP), intent(in):: xyr_Height (0:imax-1, 1:jmax, 0:kmax)
                              ! 高度 (半整数レベル). 
                              ! Height (half level)
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in):: xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! Turbulent kinetic energy (m2 s-2)
    real(DP), intent(in):: xy_SurfMomFluxX (0:imax-1, 1:jmax)
                              ! 
                              ! Eastward momentum flux at surface
    real(DP), intent(in):: xy_SurfMomFluxY (0:imax-1, 1:jmax)
                              ! 
                              ! Northward momentum flux at surface

    real(DP), intent(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
                              ! 東西方向運動量フラックス. 
                              ! Eastward momentum flux
    real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
    real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 熱フラックス. 
                              ! Heat flux
    real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
                              ! 質量フラックス. 
                              ! Mass flux of compositions
    real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数:運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数:温度. 
                              ! Diffusion coefficient: temperature
    real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
                              ! 拡散係数:比湿. 
                              ! Diffusion coefficient: specific humidity
    real(DP), intent(out):: xyz_DTurKinEneDt (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Tendency of turbulent kinetic energy

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


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


    ! 実行文 ; Executable statement
    !

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


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


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

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


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

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

      end do
    end do



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



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

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

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

    ! 作業変数
    ! Work variables
    !

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

    ! 実行文 ; Executable statement
    !

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


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


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


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

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

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

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

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

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


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

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


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


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


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

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


    ! 拡散係数の計算
    ! Calculation of diffusion coefficient
    !
    xyz_VelDiffCoef  = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sm
    xyz_TempDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sh
    !
    do k = 0, kmax
      if ( ( k == 0 ) .or. ( k == kmax ) ) then
        xyr_VelDiffCoef (:,:,k) = 0.0_DP
        xyr_TempDiffCoef(:,:,k) = 0.0_DP
      else
        xyr_VelDiffCoef (:,:,k) = ( xyz_VelDiffCoef (:,:,k) + xyz_VelDiffCoef (:,:,k+1) ) / 2.0_DP
        xyr_TempDiffCoef(:,:,k) = ( xyz_TempDiffCoef(:,:,k) + xyz_TempDiffCoef(:,:,k+1) ) / 2.0_DP
      end if
    end do
    !
    do k = 1, kmax-1
      do j = 1, jmax
        do i = 0, imax-1
          xyr_VelDiffCoef(i,j,k) = max( min( xyr_VelDiffCoef(i,j,k), VelDiffCoefMax ), VelDiffCoefMin )
          xyr_TempDiffCoef(i,j,k) = max( min( xyr_TempDiffCoef(i,j,k), TempDiffCoefMax ), TempDiffCoefMin )
        end do
      end do
    end do
    !
    xyr_QMixDiffCoef      = xyr_TempDiffCoef


    ! 輸送係数とフラックスの計算
    ! Calculate transfer coefficient and flux
    !
    call VDiffusionCalcFlux( xyz_U, xyz_V, xyzf_QMix, 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 )


    ! Calculate tendency of turbulent kinetic energy

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

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

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


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

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


    !
    ! Solve simultaneous linear equations by use of LU decomposition technique
    !
    xyza_TurKinEneLUMtx = xyza_TurKinEneMtx
    !
    call PhyImplLUDecomp3( xyza_TurKinEneLUMtx, 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 )

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


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



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

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

    call HistoryAutoPut( TimeN, 'MixLength' , xyz_MixLength )



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

  end subroutine VDiffusionMY25GBT94
Subroutine :
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
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (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_VirTemp(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
: $ hat{T}_v $ . 仮温度 (半整数レベル). Virtual temperature (half level)
xyz_Height(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
: 高度 (整数レベル). Height (full 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, 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, xyr_Press, xyz_Exner, xyr_Exner, xyr_VirTemp, xyz_Height, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef )
    !
    ! フラックス (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, GasRDry
                              ! $ R $ [J kg-1 K-1]. 
                              ! 乾燥大気の気体定数. 
                              ! Gas constant of air

    ! 時刻管理
    ! 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):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

    real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual temperature (half level)
    real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
                              ! 高度 (整数レベル). 
                              ! Height (full level)

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

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


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

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

    ! 実行文 ; Executable statement
    !

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


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


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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

  end subroutine VDiffusionOutPut

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 $
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.0_DP / 3.0_DP ) - ( 2.0_DP * MYConstA1 / MYConstB1 )
    Gamma2 =   ( MYConstB2 / MYConstB1 ) + ( 6.0_DP * MYConstA1 / MYConstB1 )
    Alpha1 = 3.0_DP  * MYConstA2 * Gamma1
    Alpha2 = 3.0_DP  * MYConstA2 * ( Gamma1 + Gamma2 )
    Beta1  = MYConstA1 * MYConstB1 * ( Gamma1 - MYConstC1 )
    Beta2  = MYConstA1 * (   MYConstB1 * ( Gamma1 - MYConstC1 ) + 6.0_DP * MYConstA1 + 3.0_DP * MYConstA2 )
    Beta3  = MYConstA2 * MYConstB1 * Gamma1
    Beta4  = MYConstA2 * (   MYConstB1 * ( Gamma1 + Gamma2 ) - 3.0_DP * MYConstA1 )
    CrtlFluxRiNum = Gamma1 / ( Gamma1 + Gamma2 )

    ! フラックス $ R_i $ 数の算出
    ! Calculate flux $ R_i $ number
    !
    xyr_FluxRiNum = (   Beta1 + Beta4 * xyr_BulkRiNum - sqrt(   ( Beta1 + Beta4 * xyr_BulkRiNum )**2 - 4.0_DP * Beta2 * Beta3 * xyr_BulkRiNum ) ) / ( 2.0_DP * Beta2 )

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

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

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

            xyr_Sh(i,j,k) = (   Alpha1 - Alpha2 * xyr_FluxRiNum(i,j,k) ) / ( 1.0_DP - 1.0_DP * xyr_FluxRiNum(i,j,k) )

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

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

          else

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

          end if

        end do
      end do
    end do


    ! 混合距離の算出
    ! Calculate mixing length
    !
    do k = 0, kmax
      xyr_MixLength(:,:,k) = FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) / (1.0_DP + FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) / MixLengthMax )
    end do

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

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

    do k = 0, kmax-1
      do i = 0, imax-1
        do j = 1, jmax
          xyr_VelDiffCoef(i,j,k) = max( min( xyr_VelDiffCoef(i,j,k), VelDiffCoefMax ), VelDiffCoefMin )
          xyr_TempDiffCoef(i,j,k) = max( min( xyr_TempDiffCoef(i,j,k), TempDiffCoefMax ), TempDiffCoefMin )
        end do
      end do
    end do

    xyr_QMixDiffCoef = xyr_TempDiffCoef

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


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

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


  end subroutine VDiffCoefficient
Subroutine :
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

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

Vertical diffusion flux is calculated.

[Source]

  subroutine VDiffusionCalcFlux( 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 )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! 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

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

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

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

    ! 実行文 ; Executable statement
    !

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


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

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

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

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

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

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

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

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

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


  end subroutine VDiffusionCalcFlux
Subroutine :
z_U(1:kmax) :real(DP), intent(in)
: $ u $ . 東西風速. Eastward wind
z_V(1:kmax) :real(DP), intent(in)
: $ v $ . 南北風速. Northward wind
zf_QMix(1:kmax, 1:ncmax) :real(DP), intent(in)
: $ q $ . 質量混合比. Mass mixing ratio
z_Temp(1:kmax) :real(DP), intent(in)
: $ T $ . 温度. Temperature
r_VirTemp(0:kmax) :real(DP), intent(in)
: $ hat{T}_v $ . 仮温度 (半整数レベル). Virtual temperature (half level)
r_Press(0:kmax) :real(DP), intent(in)
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half level)
z_Height(1:kmax) :real(DP), intent(in)
: 高度 (整数レベル). Height (full level)
z_Exner(1:kmax) :real(DP), intent(in)
: Exner 関数 (整数レベル). Exner function (full level)
r_Exner(0:kmax) :real(DP), intent(in)
: Exner 関数 (半整数レベル). Exner function (half level)
r_VelDiffCoef(0:kmax) :real(DP), intent(in)
: 拡散係数:運動量. Diffusion coefficient: velocity
r_TempDiffCoef(0:kmax) :real(DP), intent(in)
: 拡散係数:温度. Transfer coefficient: temperature
r_QMixDiffCoef(0:kmax) :real(DP), intent(in)
: 拡散係数:比湿. Diffusion coefficient: specific humidity
r_MomFluxX(0:kmax) :real(DP), intent(out)
: 東西方向運動量フラックス. Eastward momentum flux
r_MomFluxY(0:kmax) :real(DP), intent(out)
: 南北方向運動量フラックス. Northward momentum flux
r_HeatFlux(0:kmax) :real(DP), intent(out)
: 熱フラックス. Heat flux
rf_QMixFlux(0:kmax, 1:ncmax) :real(DP), intent(out)
: 質量フラックス. Mass flux of compositions

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

Vertical diffusion flux is calculated.

[Source]

  subroutine VDiffusionCalcFlux1D( z_U, z_V, zf_QMix, z_Temp, r_VirTemp, r_Press, z_Height, z_Exner, r_Exner, r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef, r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux )
    !
    ! 鉛直拡散フラックスを計算します. 
    !
    ! 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):: z_U (1:kmax)
                              ! $ u $ .   東西風速. Eastward wind
    real(DP), intent(in):: z_V (1:kmax)
                              ! $ v $ .   南北風速. Northward wind
    real(DP), intent(in):: zf_QMix(1:kmax, 1:ncmax)
                              ! $ q $ .   質量混合比. Mass mixing ratio
    real(DP), intent(in):: z_Temp (1:kmax)
                              ! $ T $ .   温度. Temperature
    real(DP), intent(in):: r_VirTemp (0:kmax)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル). 
                              ! Virtual temperature (half level)
    real(DP), intent(in):: r_Press  (0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: z_Height (1:kmax)
                              ! 高度 (整数レベル). 
                              ! Height (full level)
    real(DP), intent(in):: z_Exner (1:kmax)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP), intent(in):: r_Exner (0:kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)

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

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

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

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

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

    ! 実行文 ; Executable statement
    !

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


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

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

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

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

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

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

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

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

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


  end subroutine VDiffusionCalcFlux1D
Subroutine :
z_U(1:kmax) :real(DP), intent(in)
: $ u $ . 東西風速. Eastward wind
z_V(1:kmax) :real(DP), intent(in)
: $ v $ . 南北風速. Northward wind
zf_QMix(1:kmax, 1:ncmax) :real(DP), intent(in)
: $ q $ . 質量混合比. Mass mixing ratio
z_Temp(1:kmax) :real(DP), intent(in)
: $ T $ . 温度. Temperature
r_Temp(0:kmax) :real(DP), intent(in)
: $ hat{T} $ . 温度 (半整数レベル). Temperature (half level)
z_VirTemp(1:kmax) :real(DP), intent(in)
: $ T_v $ . 仮温度. Virtual temperature
r_VirTemp(0:kmax) :real(DP), intent(in)
: $ hat{T}_v $ . 仮温度 (半整数レベル). Virtual temperature (half level)
r_Press(0:kmax) :real(DP), intent(in)
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half level)
SurfHeight :real(DP), intent(in)
: $ z_s $ . 地表面高度. Surface height.
z_Height(1:kmax) :real(DP), intent(in)
: 高度 (整数レベル). Height (full level)
r_Height(0:kmax) :real(DP), intent(in)
: 高度 (半整数レベル). Height (half level)
z_Exner(1:kmax) :real(DP), intent(in)
: Exner 関数 (整数レベル). Exner function (full level)
r_Exner(0:kmax) :real(DP), intent(in)
: Exner 関数 (半整数レベル). Exner function (half level)
z_TurKinEne(1:kmax) :real(DP), intent(in)
: Turbulent kinetic energy (m2 s-2)
SurfMomFluxX :real(DP), intent(in)
: Eastward momentum flux at surface
SurfMomFluxY :real(DP), intent(in)
: Northward momentum flux at surface
r_MomFluxX(0:kmax) :real(DP), intent(out)
: 東西方向運動量フラックス. Eastward momentum flux
r_MomFluxY(0:kmax) :real(DP), intent(out)
: 南北方向運動量フラックス. Northward momentum flux
r_HeatFlux(0:kmax) :real(DP), intent(out)
: 熱フラックス. Heat flux
rf_QMixFlux(0:kmax, 1:ncmax) :real(DP), intent(out)
: 質量フラックス. Mass flux of compositions
r_VelDiffCoef(0:kmax) :real(DP), intent(out)
: 拡散係数:運動量. Diffusion coefficient: velocity
r_TempDiffCoef(0:kmax) :real(DP), intent(out)
: 拡散係数:温度. Diffusion coefficient: temperature
r_QMixDiffCoef(0:kmax) :real(DP), intent(out)
: 拡散係数:比湿. Diffusion coefficient: specific humidity
z_DTurKinEneDt(1:kmax) :real(DP), intent(out)
: Tendency of turbulent kinetic energy

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

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

[Source]

  subroutine VDiffusionMY251D( z_U, z_V, zf_QMix, z_Temp, r_Temp, z_VirTemp, r_VirTemp, r_Press, SurfHeight, z_Height, r_Height, z_Exner, r_Exner, z_TurKinEne, SurfMomFluxX, SurfMomFluxY, r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux, r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef, z_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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

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

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

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

    ! 作業変数
    ! Work variables
    !

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    ! 実行文 ; Executable statement
    !

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


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


    ! Calculate turbulent kinetic energy with offset
    !
    z_TurKinEneNonZero = z_TurKinEne + TurKinEneOffset

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

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

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


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


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

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



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


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


    ! Calculate tendency of turbulent kinetic energy

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

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

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


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


    nloop = kmax
    loop_solve : do iloop = 1, nloop

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

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


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

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

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


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

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

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

    end do loop_solve



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

  end subroutine VDiffusionMY251D
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-20150213 $’ // ’$Id: vdiffusion_my.f90,v 1.5 2015/01/29 12:09:17 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version