!= Module Turbulence_3d
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA, Masatsugu
! Version::   $Id: turbulence_kw1978.f90,v 1.18 2012-12-29 16:27:35 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!モデルの物理過程を計算するために必要となる関数群を束ねたモジュール
!具体的には以下の項を計算するための関数を格納する.  
!  * 乱流拡散項
!  * 乱流エネルギーの時間発展方程式に含まれる各項
!  * 散逸加熱項
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!  * 1.5 次のクロージャーを利用する
!  * 定式化は CReSS のマニュアルを参照した
!
!== Future Plans
!
!

module Turbulence_kw1978
  !
  !モデルの物理過程を計算するために必要となる関数群を束ねたモジュール
  !具体的には以下の項を計算するための関数を格納する.  
  !  * 乱流拡散項
  !  * 乱流エネルギーの時間発展方程式に含まれる各項
  !  * 散逸加熱項
  !

  !モジュール読み込み 
  use dc_types, only: DP, STRING
  use dc_iounit,  only: FileOpen
  use dc_message, only: MessageNotify
  use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut

  use mpi_wrapper,only: myrank
  use timeset, only:  DelTimeLong, TimeN
  use gridset, only:  imin,           &! x 方向の配列の下限
    &                 imax,           &! x 方向の配列の上限
    &                 jmin,           &! y 方向の配列の下限
    &                 jmax,           &! y 方向の配列の上限
    &                 kmin,           &! z 方向の配列の下限
    &                 kmax,           &! z 方向の配列の上限
    &                 nx,ny,nz,ncmax      ! 
  use basicset, only: xyz_PTempBZ,    &!基本場の温位
    &                 xyzf_QMixBZ,    &!基本場の混合比
    &                 xyz_QMixBZ,    &!基本場の混合比
    &                 xyz_QMixBZPerMolWt,    &!基本場の混合比
    &                 xyz_EffMolWtBZ, &!基本場の混合比
    &                 xyz_VelSoundBZ,  &! 音速
    &                 xyz_VPTempBZ,    &! 仮温位
    &                 xyz_ExnerBZ      !エクスナー関数の基本場
  use constants,only: Grav,          &
    &                 MolWtDry,      & 
    &                 CpDry
  use composition, only: MolWtWet,      &
    &                 SpcWetSymbol,      &
    &                 SpcWetID,      &
    &                 CondNum,          &!凝結過程の数
    &                 IdxCG,            &!凝結過程(蒸気)の配列添え字
    &                 IdxCC,            &!凝結過程(雲)の配列添え字
!    &                 IdxNH3,           &!NH3(蒸気)の配列添え字
!    &                 IdxH2S,           &!H2S(蒸気)の配列添え字
    &                 GasNum,           &!気体の数
    &                 IdxG               !気体の配列添え字
  use ChemCalc, only: xyz_LatentHeat     !潜熱
!    &                 ReactHeatNH4SH     !NH4SH の反応熱
  use axesset, only: &  
    &                 x_dx,           &! x 方向の格子点間隔
    &                 y_dy,           &! y 方向の格子点間隔
    &                 z_dz,           &! z 方向の格子点間隔
    &                 xyz_avr_pyz, xyr_avr_pyr, xqz_avr_pqz, &
    &                 pyz_avr_xyz, pyr_avr_xyr, pqz_avr_xqz, &
    &                 xyz_avr_xqz, pyz_avr_pqz, xyr_avr_xqr, &
    &                 xqz_avr_xyz, pqz_avr_pyz, xqr_avr_xyr, &
    &                 xyz_avr_xyr, pyz_avr_pyr, xqz_avr_xqr, &
    &                 xyr_avr_xyz, pyr_avr_pyz, xqr_avr_xqz, &
    &                 pqz_avr_xyz, pyr_avr_xyz, xqr_avr_xyz, &
    &                 xyz_avr_pqz, xyz_avr_pyr, xyz_avr_xqr
  use xyz_deriv_c4_module, only: pyz_c4dx_xyz, xqz_c4dy_xyz, xyr_c4dz_xyz
  use xyz_deriv_module, only: xyz_dx_pyz, xyz_dy_xqz, xyz_dz_xyr, &
    &                         pyz_dx_xyz, xqz_dy_xyz, xyr_dz_xyz, &
    &                         pqz_dx_xqz, pqz_dy_pyz, pyz_dy_pqz, &
    &                         pyr_dx_xyr, pyz_dz_pyr, pyr_dz_pyz, &
    &                         pyr_dz_pyz, pyr_dx_xyr, xyr_dx_pyr, &
    &                         xqr_dz_xqz, xqr_dy_xyr, xyr_dy_xqr, &
    &                         pqz_dy_pyz, pqz_dx_xqz, xqz_dx_pqz, &
    &                         xqr_dy_xyr, xqr_dz_xqz, xqz_dz_xqr  
  use namelist_util, only: namelist_filename
  use setmargin, only: SetMargin_xyz
  use DExnerDt,  only: xyz_DExnerDt_xyz

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public Turbulence_KW1978_Init
  public Turbulence_KW1978_Forcing

  !変数定義
  real(DP), save :: Cm     = 2.0d-1    !乱流エネルギー診断式の係数 
  real(DP), save :: MixLen = 0.0d0     !平均混合距離
  real(DP), save :: KmMax  = 0.0d0     !最大値
  real(DP), save :: ConstKm = 0.0d0    !定数値
  real(DP), save :: ConstKh = 0.0d0    !定数値
  logical,  save :: FlagConstTurbCoef
  logical,  save :: FlagDExnerDtTurb =.true.

  character(*), parameter:: module_name = 'turbulence_kw1978'


contains

!!!------------------------------------------------------------------------!!!
  subroutine turbulence_kw1978_init
    !
    ! Turbulence モジュールの初期化ルーチン
    ! 

    !暗黙の型宣言禁止
    implicit none

    real(DP):: MinDelX
    real(DP):: MinDelY
    real(DP):: MinDelZ
    integer :: l
    integer :: unit

    ! NAMELIST から情報を取得
    NAMELIST /turbulence_kw1978_nml/ &
      & Cm, KmMax,                   &
      & FlagConstTurbCoef, ConstKm, ConstKh, &
      & FlagDExnerDtTurb

    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=turbulence_kw1978_nml)
    close(unit)

    ! 混合距離
    ! 2 次元計算の場合には DelY に依存しないようにするために if 文を利用.
    ! 
    MinDelX = minval(x_dx)
    MinDelY = minval(y_dy)
    MinDelZ = minval(z_dz)

    if (jmin == jmax) then 
!       MixLen = min( MinDelZ, MinDelX)
       MixLen = sqrt(MinDelX * MinDelZ) 
    else
!       MixLen = min( MinDelZ, min( MinDelX, MinDelY ) )
       MixLen = ( MinDelZ * MinDelX * MinDelY ) ** (1.0 / 3.0)
    end if

    
    ! KmMax が設定されていない場合. 
    ! 安定性解析では, dt / l**2 < 0.5 を満たす必要がある. ここでは 0.1 にしておいた. 
    !
    if (KmMax == 0.0d0) then 
       KmMax = 0.1 * (MixLen ** 2.0d0) / (DelTimeLong * 2.0d0) !LeapFrog
    end if

    ! Output
    !
    if (myrank == 0) then 
      call MessageNotify( "M", &
        & module_name, "Cm = %f", d=(/Cm/))
      call MessageNotify( "M", &
        & module_name, "KmMax = %f", d=(/KmMax/))
      call MessageNotify( "M", &
        & module_name, "MixLen= %f", d=(/MixLen/))
      call MessageNotify( "M", &
        & module_name, "FlagConstTurbCoef= %b", l=(/ FlagConstTurbCoef /))
      call MessageNotify( "M", &
        & module_name, "ConstKm= %f", d=(/ConstKm/))
      call MessageNotify( "M", &
        & module_name, "ConstKh= %f", d=(/ConstKh/))
      call MessageNotify( "M", &
        & module_name, "FlagDExnerDtTurb= %b", l=(/ FlagDExnerDtTurb /))
    end if

    call HistoryAutoAddVariable(  &
      & varname='VelXTurb',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Turbulence term of velocity (x)', &
      & units='K.s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelYTurb',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Turbulence term of velocity (y)', &
      & units='K.s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelZTurb',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Turbulence term of velocity (z)', &
      & units='K.s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='PTempTurb',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Turbulence term of potential temperature', &
      & units='K.s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='PTempDisp',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Dissipation term of potential temperature', &
      & units='K.s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='ExnerDisp',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Dissipation term of exner function', &
      & units='K.s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='CDensTurb',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Turbulence term of cloud density', &
      & units='kg.m-3.s-1',    &
      & xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_Turb', & 
        & dims=(/'x','y','z','t'/),     &
        & longname='Turbulence term of '          &
        &           //trim(SpcWetSymbol(l))//' mixing ratio',  &
        & units='kg.kg-1.s-1',    &
        & xtype='float')
    end do

    call HistoryAutoAddVariable(  &
      & varname='KmDiff',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Diffusion term of Km', &
      & units='s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='KmBuoyT',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Buoyancy term of Km (temperature)', &
      & units='s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='KmBuoyM',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Buoyancy term of Km (composition)', &
      & units='s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='KmShear',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Shear term of Km', &
      & units='s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='KmDisp',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Dissipation term of Km', &
      & units='s-1',    &
      & xtype='float')


  end subroutine turbulence_kw1978_init

!!!------------------------------------------------------------------------!!!
  subroutine turbulence_KW1978_forcing(        &
    & pyz_VelXBl, xqz_VelYBl,   xyr_VelZBl,    &
    & xyz_PTempBl, xyz_ExnerBl, xyzf_QMixBl,   &
    & xyz_KmBl,    xyz_KhBl,    xyz_CDensBl,   &
    & pyz_DVelXDt, xqz_DVelYDt,  xyr_DVelZDt,  &
    & xyz_DPTempDt,xyz_DExnerDt, xyzf_DQMixDt, &
    & xyz_DKmDt,   xyz_DCDensDt,               &
    & xyz_KmAl, xyz_KhAl &
    )

    implicit none

    real(DP),intent(in) :: pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平速度
    real(DP),intent(in) :: xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平速度
    real(DP),intent(in) :: xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !鉛直速度
    real(DP),intent(in) :: xyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !温位
    real(DP),intent(in) :: xyz_ExnerBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !無次元圧力
    real(DP),intent(in) :: xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                                    !凝縮成分の混合比
    real(DP),intent(in) :: xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP),intent(in) :: xyz_KhBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP),intent(in) :: xyz_CDensBl(imin:imax,jmin:jmax,kmin:kmax)

    real(DP),intent(inout):: pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP),intent(inout):: xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP),intent(inout):: xyr_DVelZDt(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP),intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax)

    real(DP),intent(inout):: xyz_DExnerDt(imin:imax,jmin:jmax,kmin:kmax)

    real(DP),intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax)

    real(DP),intent(inout):: xyz_DKmDt(imin:imax,jmin:jmax,kmin:kmax)

    real(DP),intent(inout) :: xyz_DCDensDt(imin:imax,jmin:jmax,kmin:kmax)
    
    real(DP),intent(out):: xyz_KmAl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP),intent(out):: xyz_KhAl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP)            :: xyz_Buoy(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数の
    real(DP)            :: xyz_BuoyT(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数の
    real(DP)            :: xyz_BuoyM(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数の
    real(DP)            :: xyz_Shear(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数の
    real(DP)            :: xyz_Diff(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数の
    real(DP)            :: xyz_Disp(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流エネルギーの消散
    real(DP)            :: xyz_DispPI(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流エネルギーの消散
    real(DP)            :: xyz_DispHeat(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流エネルギーの消散
    real(DP)            :: xyz_Turb(imin:imax,jmin:jmax,kmin:kmax)
                                                    !
    real(DP)            :: xyzf_Turb(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: pyz_Turb(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xqz_Turb(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xyr_Turb(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: pyz_DVelXDt0(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: xqz_DVelYDt0(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: xyr_DVelZDt0(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xyz_DExnerDt0(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, ncmax)

    real(DP)            :: xyz_DKmDt0(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xyz_DCDensDt0(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xyz_PTempBlAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: xyzf_QMixBlAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    integer             :: f


    !----------------------------------
    ! 初期化
    !
    xyz_PTempBlAll = xyz_PTempBl + xyz_PTempBZ
    xyzf_QMixBlAll  = xyzf_QMixBl + xyzf_QMixBZ 
    pyz_DVelXDt0 = pyz_DVelXDt
    xqz_DVelYDt0 = xqz_DVelYDt
    xyr_DVelZDt0 = xyr_DVelZDt
    xyz_DPTempDt0 = xyz_DPTempDt
    xyz_DExnerDt0 = xyz_DExnerDt
    xyzf_DQMixDt0 = xyzf_DQMixDt
    xyz_DKmDt0    = xyz_DKmDt
    xyz_DCDensDt0 = xyz_DCDensDt
    
    !----------------------------------
    ! 拡散係数の時間発展 (エネルギー方程式を Km の式に変形したもの)
    !
    if ( .not. FlagConstTurbCoef ) then
    ! Buoyancy term
    !
    xyz_Buoy = xyz_BuoyMoistKm(xyz_PTempBl, xyz_ExnerBl, xyzf_QMixBl) 
    xyz_BuoyT = &
        &  - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) &
        &       * xyz_avr_xyr( xyr_dz_xyz( xyz_PTempBlAll ) ) &
        &       / ( 2.0d0 * xyz_PTempBZ )
    xyz_BuoyM = xyz_Buoy - xyz_BuoyT


    xyz_Shear = &
      &   ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 )                    &
      & * (                                                        &
      &      ( xyz_dx_pyz( pyz_VelXBl ) ) ** 2.0d0                 &
      &    + ( xyz_dy_xqz( xqz_VelYBl ) ) ** 2.0d0                 &
      &    + ( xyz_dz_xyr( xyr_VelZBl ) ) ** 2.0d0                 &
      &    + 5.0d-1                                                &
      &      * (                                                   &
      &          (                                                 &
      &              xyz_avr_pyr( pyr_dz_pyz( pyz_VelXBl ) )       &
      &            + xyz_avr_pyr( pyr_dx_xyr( xyr_VelZBl ) )       &
      &           ) ** 2.0d0                                       &
      &        + (                                                 &
      &              xyz_avr_xqr( xqr_dy_xyr( xyr_VelZBl ) )       &
      &            + xyz_avr_xqr( xqr_dz_xqz( xqz_VelYBl ) )       &
      &           ) ** 2.0d0                                       &
      &        + (                                                 &
      &              xyz_avr_pqz( pqz_dx_xqz( xqz_VelYBl ) )       &
      &            + xyz_avr_pqz( pqz_dy_pyz( pyz_VelXBl ) )       &
      &           ) ** 2.0d0                                       &
      &        )                                                   &
      &   )                                                        &
      & - xyz_KmBl * (  xyz_dx_pyz( pyz_VelXBl )                   &
      &               + xyz_dy_xqz( xqz_VelYBl )                   &
      &               + xyz_dz_xyr( xyr_VelZBl ) ) / 3.0d0

    xyz_Diff =                                                     &
      &   5.0d-1                                                   &
      &   * (                                                      &
      &         xyz_dx_pyz(pyz_dx_xyz(xyz_KmBl ** 2.0d0))          &
      &       + xyz_dy_xqz(xqz_dy_xyz(xyz_KmBl ** 2.0d0))          &
      &       + xyz_dz_xyr(xyr_dz_xyz(xyz_KmBl ** 2.0d0))          &
      &     )                                                      &
      & + (                                                        &
      &      (xyz_avr_pyz(pyz_dx_xyz(xyz_KmBl))) ** 2.0d0          &
      &    + (xyz_avr_xqz(xqz_dy_xyz(xyz_KmBl))) ** 2.0d0          &
      &    + (xyz_avr_xyr(xyr_dz_xyz(xyz_KmBl))) ** 2.0d0          &
      &   )
    
    ! t - \Delta t で評価
    !
    xyz_Disp = - (xyz_KmBl ** 2.0d0) * 5.0d-1 / (MixLen ** 2.0d0)

    ! tendency
    !
    xyz_DKmDt = xyz_DKmDt0 + xyz_Buoy + xyz_Shear + xyz_Diff + xyz_Disp

    call HistoryAutoPut(TimeN, 'KmBuoyT', xyz_BuoyT(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'KmBuoyM', xyz_BuoyM(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'KmShear', xyz_Shear(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'KmDiff', xyz_Diff(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'KmDisp', xyz_Disp(1:nx,1:ny,1:nz))    
   
    ! 時間積分
    !
    xyz_KmAl = xyz_KmBl + (2.0d0 * DelTimeLong) * xyz_DKmDt

    ! 上限値の設定
    !
    xyz_KmAl = max( 0.0d0, min( xyz_KmAl, KmMax ) )

    ! Kh の設定
    !
    xyz_KhAl = 3.0d0 * xyz_KmAl

    else 
      xyz_KmAl = ConstKm
      xyz_KhAl = ConstKh
    end if
    
    !--------------------------------
    ! 雲密度の tendency
    !
    xyz_Turb =                                                              &
      &   xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyz_CDensBl ) ) &
      & + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyz_CDensBl ) ) &
      & + xyz_dz_xyr( xyr_avr_xyz( xyz_KhBl ) * xyr_dz_xyz( xyz_CDensBl ) )   

    xyz_DCDensDt = xyz_DCDensDt0 + xyz_Turb

    ! output
    !
    call HistoryAutoPut(TimeN, 'CDensTurb', xyz_Turb(1:nx,1:ny,1:nz))    

    !--------------------------------
    ! 温位の tendency
    !
    xyz_Turb =                                                                  &
      &   xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyz_PTempBlAll ) )  &
      & + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyz_PTempBlAll ) )  &
      & + xyz_dz_xyr( xyr_avr_xyz( xyz_KhBl ) * xyr_dz_xyz( xyz_PTempBlAll ) )  
    
    xyz_DispHeat = (xyz_KmBl ** 3.0d0) &
      & / (xyz_ExnerBZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0))

    xyz_DPTempDt = xyz_DPTempDt0 + xyz_Turb + xyz_DispHeat

    call HistoryAutoPut(TimeN, 'PTempDisp', xyz_DispHeat(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'PTempTurb', xyz_Turb(1:nx, 1:ny, 1:nz))

    !--------------------------------
    ! 混合比の tendency
    !
    do f = 1, ncmax    
      xyzf_Turb(:,:,:,f) =                                                     &
        &   xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyzf_QMixBlAll(:,:,:,f) ) )  &
        & + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyzf_QMixBlAll(:,:,:,f) ) )  &
        & + xyz_dz_xyr( xyr_avr_xyz( xyz_KhBl ) * xyr_dz_xyz( xyzf_QMixBlAll(:,:,:,f) ) )    
      
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_Turb', xyzf_Turb(1:nx,1:ny,1:nz,f))
    end do

    xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Turb
    
    !--------------------------------
    ! 各速度成分の tendency
    !
    pyz_Turb = &
      &   2.0d0 * pyz_dx_xyz( xyz_KmBl * xyz_dx_pyz( pyz_VelXBl ) )&
      & + pyz_dy_pqz(                                              &
      &       pqz_avr_xyz( xyz_KmBl ) * pqz_dx_xqz( xqz_VelYBl )   &
      &     + pqz_avr_xyz( xyz_KmBl ) * pqz_dy_pyz( pyz_VelXBl )   &
      &   )                                                        &
      & + pyz_dz_pyr(                                              &
      &       pyr_avr_xyz( xyz_KmBl ) * pyr_dx_xyr( xyr_VelZBl )   &
      &     + pyr_avr_xyz( xyz_KmBl ) * pyr_dz_pyz( pyz_VelXBl )   &
      &   )                                                        &
      & - 2.0d0 * pyz_dx_xyz( ( xyz_KmBl ** 2.0d0 ) )              &
      &   / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )

    pyz_DVelXDt = pyz_DVelXDt0 + pyz_Turb

    call HistoryAutoPut(TimeN, 'VelXTurb', pyz_Turb(1:nx, 1:ny, 1:nz))

    xqz_Turb = &
      &   2.0d0 * xqz_dy_xyz( xyz_KmBl * xyz_dy_xqz( xqz_VelYBl ) ) &
      & + xqz_dx_pqz(                                             &
      &       pqz_avr_xyz( xyz_KmBl ) * pqz_dy_pyz( pyz_VelXBl )  &
      &     + pqz_avr_xyz( xyz_KmBl ) * pqz_dx_xqz( xqz_VelYBl )  &
      &   )                                                       &
      & + xqz_dz_xqr(                                             &
      &       xqr_avr_xyz( xyz_KmBl ) * xqr_dy_xyr( xyr_VelZBl )  &
      &     + xqr_avr_xyz( xyz_KmBl ) * xqr_dz_xqz( xqz_VelYBl )  &
      &   )                                                       &
      & - 2.0d0 * xqz_dy_xyz( ( xyz_KmBl ** 2.0d0 ) )             &
      &   / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )
    
    xqz_DVelYDt = xqz_DVelYDt0 + xqz_Turb

    call HistoryAutoPut(TimeN, 'VelYTurb', xqz_Turb(1:nx, 1:ny, 1:nz))

    xyr_Turb = &
      & + 2.0d0 * xyr_dz_xyz( xyz_KmBl * xyz_dz_xyr( xyr_VelZBl ) ) &
      & + xyr_dx_pyr(                                               &
      &      pyr_avr_xyz( xyz_KmBl ) * pyr_dz_pyz( pyz_VelXBl )     &
      &    + pyr_avr_xyz( xyz_KmBl ) * pyr_dx_xyr( xyr_VelZBl )     &
      &   )                                                         &
      & + xyr_dy_xqr(                                               &
      &      xqr_avr_xyz( xyz_KmBl ) * xqr_dz_xqz( xqz_VelYBl )     &
      &    + xqr_avr_xyz( xyz_KmBl ) * xqr_dy_xyr( xyr_VelZBl )     &
      &   )                                                         & 
      & - 2.0d0 * xyr_dz_xyz( ( xyz_KmBl ** 2.0d0 ) )               &
      &   / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )

    xyr_DVelZDt = xyr_DVelZDt0 + xyr_Turb

    call HistoryAutoPut(TimeN, 'VelZTurb', xyr_Turb(1:nx, 1:ny, 1:nz))

    !--------------------
    ! Exner function
    !
    if ( FlagDExnerDtTurb ) then
      xyz_DispPI = xyz_DExnerDt_xyz( xyz_DispHeat )
    else 
      xyz_DispPi = 0.0d0
    end if
    xyz_DExnerDt = xyz_DExnerDt0 + xyz_DispPI

    call HistoryAutoPut(TimeN, 'ExnerDisp', xyz_DispPI(1:nx, 1:ny, 1:nz))

    ! Set Margin
    !
    call SetMargin_xyz(xyz_KmAl)
    call SetMargin_xyz(xyz_KhAl)

  end subroutine Turbulence_KW1978_forcing


!!!------------------------------------------------------------------------!!!
  function xyz_BuoyMoistKm(xyz_PTemp, xyz_Exner, xyzf_QMix)
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax)
                                               !温位
    real(DP), intent(in) :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)
                                               !無次元圧力
    real(DP), intent(in) :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                               !凝縮成分の混合比
    real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                               !凝縮成分の混合比
    real(DP) :: xyzf_QMixAll2(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                               !凝縮成分の混合比
    real(DP) :: xyz_BuoyMoistKm(imin:imax,jmin:jmax,kmin:kmax)
                                               !
    real(DP) :: xyzf_LatentHeat(imin:imax,jmin:jmax,kmin:kmax, GasNum)
                                               !潜熱
    real(DP) :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax)
                                               !温度
    real(DP) :: xyz_EffHeat(imin:imax,jmin:jmax,kmin:kmax)
                                               !
    real(DP) :: xyz_EffPTemp(imin:imax,jmin:jmax,kmin:kmax)    
                                               !浮力に対する温度差の寄与
    real(DP) :: xyz_EffMolWt(imin:imax,jmin:jmax,kmin:kmax)    
                                               !浮力に対する分子量差の寄与
    real(DP) :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, GasNum)
                                               !混合比/分子量
    integer              :: s
    
    
    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xyz_TempAll = ( xyz_PTemp + xyz_PTempBZ ) * ( xyz_Exner + xyz_ExnerBZ )
    xyzf_QMixAll = xyzf_QMixBZ + xyzf_QMix
    xyzf_LatentHeat = 0.0d0
    
    !作業配列の初期化. 気体のみ利用
    do s = 1, GasNum
      xyzf_QMixPerMolWt(:,:,:,s) = xyzf_QMix(:,:,:,IdxG(s)) / MolWtWet(IdxG(s))
    end do
    
    !温度の効果
    xyz_EffPTemp = xyz_PTemp / xyz_PTempBZ 
    
    !分子量効果 + 引きづりの効果
    xyz_EffMolWt =                                           &
      & + sum(xyzf_QMixPerMolWt, 4)                         &
      &    / ( 1.0d0 / MolWtDry + xyz_QMixBZPerMolWt )  &
      & - sum(xyzf_QMix, 4) / ( 1.0d0 + xyz_QMixBZ )
    
    !蒸気が蒸発する場合の潜熱を計算
    !  分子量の部分はいつでも効くが潜熱は飽和していないと効かないので, 
    !  雲の混合比がゼロの時には, 潜熱の寄与はゼロとなるように調節している
    xyzf_QMixAll2 = xyzf_QMixAll
!    xzya_QMixAll2(:,:,:,IdxNH3) = &
!      &  xyzf_QMixAll(:,:,:,IdxNH3) - xyzf_QMixAll(:,:,:,IdxH2S) 

    do s = 1, CondNum
      xyzf_LatentHeat(:,:,:,s) =                             &
        & xyz_LatentHeat( SpcWetID(IdxCC(s)), xyz_TempAll )  &
        &  * xyzf_QMixAll2(:,:,:,IdxCG(s))                    &
        &  * ( 5.0d-1 + sign( 5.0d-1, (xyzf_QMixAll2(:,:,:,IdxCC(s)) - 1.0d-4) ) )
    end do

    xyz_EffHeat = ( sum( xyzf_LatentHeat, 4 ) * xyz_EffMolWtBZ &
!      &            + ReactHeatNH4SH * &
!      &            (xyzf_QMixAll(:,:,:,IdxH2S) + xyzf_QMixAll(:,:,:,IdxH2S)) &
      &          ) / ( CpDry * xyz_ExnerBZ ) 

    
    !乱流拡散係数の時間発展式の浮力項を決める
    xyz_BuoyMoistKm = &
      &  - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 )     &
      &    * xyz_avr_xyr(                                           &
      &        xyr_dz_xyz(                                          &
      &           xyz_EffHeat                                       &
      &         + xyz_PTempBZ / xyz_EffMolWtBZ                      &
      &           * ( 1.0d0 + xyz_EffPTemp + xyz_EffMolWt )         &
      &         )                                                   &
      &      )                                                      &
      &    / ( 2.0d0 * xyz_PTempBZ / xyz_EffMolWtBZ)                   

  end function xyz_BuoyMoistKm

end module Turbulence_kw1978
