!= Module Turbulence_3d
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA, Masatsugu
! Version::   $Id: turbulence_kw1978.f90,v 1.1 2011-04-26 04:08:10 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 gridset,   only:imin,          &! x 方向の配列の下限
    &                 imax,           &! x 方向の配列の上限
    &                 jmin,           &! y 方向の配列の下限
    &                 jmax,           &! y 方向の配列の上限
    &                 kmin,           &! z 方向の配列の下限
    &                 kmax,           &! z 方向の配列の上限
    &                 nx,ny,nz,nf      ! 
  use axesset, only:  x_dx,              &! x 方向の格子点間隔
    &                 y_dy,              &! y 方向の格子点間隔
    &                 z_dz                ! z 方向の格子点間隔
  use constants, only: CpDry,            &!乾燥成分の比熱
    &                 MolWtDry,          &!乾燥成分の分子量
    &                 Grav,              &!重力加速度
    &                 NuH, NuV
  use basicset, only: xyz_PotTempBasicZ, &!基本場の温位
    &                 xyz_EffMolWtBasicZ, &!基本場の分子量効果
    &                 xyz_ExnerBasicZ,    &!基本場の無次元圧力
    &                 xyza_MixRtBasicZ     !基本場の混合比
  use axesset, only: &
    &                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 moistset, only: CondNum,          &!凝結過程の数
    &                 SpcWetID,         &!凝縮成分のID
    &                 MolWtWet,         &!凝縮成分の分子量
    &                 SpcWetSymbol,     &!凝縮成分の名前
    &                 IdxCG,            &!凝結過程(蒸気)の配列添え字
    &                 IdxCC,            &!凝結過程(雲)の配列添え字
    &                 GasNum,           &!気体の数
    &                 IdxG,             &!気体の配列添え字
    &                 IdxNH3,           &!NH3(蒸気)の配列添え字
    &                 IdxH2S             !H2S(蒸気)の配列添え字
  use ChemCalc, only: xyz_LatentHeat,   &!潜熱
    &                 ReactHeatNH4SH     !NH4SH の反応熱

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public Turbulence_KW1978_Init
  public Turbulence_KW1978_
  public Steady_Km_KW1978
  public Tendency_PotTemp_KW1978
  public Tendency_MixRt_KW1978
  public Tendency_VelX_KW1978
  public Tendency_VelY_KW1978
  public Tendency_VelZ_KW1978

  !変数定義
  real(DP), save, public:: KmMax = 800.0d0    !上限値 (デフォルト値は 800.0)
                                              !  中島健介(1994, 学位論文)参照
  real(DP), save        :: Cm     = 2.0d-1    !乱流エネルギー診断式の係数 
  real(DP), save        :: MixLen = 0.0d0     !平均混合距離
  real(DP), allocatable :: xyz_MixRtBasicZPerMolWt(:,:,:)
                              !基本場の混合比 / 分子量 の総和
  real(DP), allocatable :: xyz_MixRtBasicZ(:,:,:)
                              !基本場の混合比の総和
  save xyz_MixRtBasicZPerMolWt
  save xyz_MixRtBasicZ

contains

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

    use gtool_historyauto, only: HistoryAutoAddVariable
    use dc_iounit,         only: FileOpen

    !暗黙の型宣言禁止
    implicit none

    character(STRING), intent(in) :: cfgfile
    real(DP), allocatable :: xyza_MixRtBasicZPerMolWt(:,:,:,:) !基本場の混合比 / 分子量
    real(DP)              :: MinDelX
    real(DP)              :: MinDelY
    real(DP)              :: MinDelZ
    integer               :: unit
    integer               :: l

    !-----------------------------------------------------------------
    ! NAMELIST から情報を取得
    !
    NAMELIST /turbulence_kw1978_nml/ KmMax
    
    call FileOpen(unit, file=cfgfile, mode='r')
    read(unit, NML=turbulence_kw1978_nml)
    close(unit)

    !混合距離
    MinDelX = minval(x_dx)
    MinDelY = minval(y_dy)
    MinDelZ = minval(z_dz)

!    MixLen = ( MinDelX * MinDelY * MinDelZ ) ** (1.0d0/3.0d0)

    MixLen = min( MinDelZ, min( MinDelX, MinDelY ) )

    !-----------------------------------------------------------
    ! 作業配列の初期化. 
    !-----------------------------------------------------------
    allocate( &
      & xyza_MixRtBasicZPerMolWt(imin:imax,jmin:jmax, kmin:kmax, GasNum),  &
      & xyz_MixRtBasicZPerMolWt(imin:imax,jmin:jmax,kmin:kmax),           &
      & xyz_MixRtBasicZ(imin:imax,jmin:jmax,kmin:kmax)                    &
      & )
    
    xyza_MixRtBasicZPerMolWt = 0.0d0
    xyz_MixRtBasicZPerMolWt = 0.0d0
    xyz_MixRtBasicZ = 0.0d0

    do l = 1, GasNum
      xyza_MixRtBasicZPerMolWt(:,:,:,l) = &
        & xyza_MixRtBasicZ(:,:,:,IdxG(l)) / MolWtWet(IdxG(l))
    end do

    xyz_MixRtBasicZPerMolWt = sum(xyza_MixRtBasicZPerMolWt, 4) 
    xyz_MixRtBasicZ         = sum(xyza_MixRtBasicZ, 4) 

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

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

    do l = 1, nf
      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='double')
    end do
    
  end subroutine turbulence_init
  
!!!------------------------------------------------------------------------!!!
  subroutine Tendency_Km_KW1978(pyz_VelXBl, pyz_VelXNl, xqz_VelYBl, xqz_VelYNl, &
    &                           xyr_VelZBl, xyr_VelZNl, xyz_PotTempBl, xyz_ExnerBl, &
    &                           xyza_MixRtBl, xyz_KmBl, xyz_KmNl, xyz_KmTend)
    ! 乱流エネルギーの時間発展を解く. Klemp and Wilhelmson (1978) 
    !    

    use xyz_deriv_module, only: pyz_dx_xyz, xqz_dy_xyz, xyr_dz_xyz
    
    implicit none

    real(DP), intent(in)  :: pyz_VelXBl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: pyz_VelXNl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xqz_VelYBl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xqz_VelYNl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZBl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZNl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyz_PotTempBl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyz_ExnerBl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyza_MixRtBl(imin:imax, jmin:jmax, kmin:kmax, nf)
    real(DP), intent(in)  :: xyz_KmBl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyz_KmNl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(out) :: xyz_KmTend(imin:imax, jmin:jmax, kmin:kmax)
    
    real(DP) :: xyz_Adv(imin:imax, jmin:jmax, kmin:kmax)
    real(DP) :: xyz_Shr(imin:imax, jmin:jmax, kmin:kmax)
    real(DP) :: xyz_Diff(imin:imax, jmin:jmax, kmin:kmax)
    real(DP) :: xyz_Buoy(imin:imax, jmin:jmax, kmin:kmax)
    real(DP) :: xyz_Disp(imin:imax, jmin:jmax, kmin:kmax)

!    xyz_Adv  = xyz_AdvKm(xyz_KmNl, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) 
    xyz_Adv  =                                         &
      & - xyz_avr_pyz(pyz_VelXNl * pyz_dx_xyz(xyz_KmNl))  &
      & - xyz_avr_xqz(xqz_VelYNl * xqz_dy_xyz(xyz_KmNl))  &
      & - xyz_avr_xyr(xyr_VelZNl * xyr_dz_xyz(xyz_KmNl))    

!!!!!!!!!!! mosit buoy !!!!!!!!!!!
    write(*,*) "bug ( mosit buoy )"

!    if (sw_moist) then 
      xyz_Buoy = xyz_BuoyMoistKm(xyz_PotTempBl, xyz_ExnerBl, xyza_MixRtBl)       
!    else
!      xyz_Buoy = xyz_BuoyKm(xyz_PotTempBl)
!    end if

    xyz_Shr  = xyz_ShearKm(xyz_KmBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl) 
    xyz_Diff = xyz_NumDiffKm(xyz_KmBl)
    xyz_Disp = xyz_DispKm(xyz_KmBl) 

    xyz_KmTend = xyz_Adv + xyz_Buoy + xyz_Shr + xyz_Diff + xyz_Disp

  end subroutine Tendency_Km_KW1978
  
!!!------------------------------------------------------------------------!!!
  subroutine Steady_Km_KW1978(pyz_VelX, xqz_VelY, xyr_VelZ, xyz_PotTemp, xyz_Km, xyz_KmA)
    !
    ! 乱流エネルギーを定常とした場合の渦拡散係数, 渦粘性係数を求める. 
    ! この乱流パラメタリゼーションは, Mellor and Yamada (1974) の
    ! Level 1 Closure に対応するが, Level1 Closure に存在する 
    ! \bar{\theta^2} は無視されている. 

    !--- 暗黙の型宣言禁止
    implicit none

    real(DP),intent(in)  :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP),intent(in)  :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP),intent(in)  :: xyr_VelZ(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP),intent(in)  :: xyz_PotTemp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP),intent(in)  :: xyz_Km(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP),intent(out) :: xyz_KmA(imin:imax,jmin:jmax,kmin:kmax) 

    xyz_KmA = &
      & sqrt ( &
      &   max( 0.0d0, &
      &        ( xyz_ShearKm( xyz_Km, pyz_VelX, xqz_VelY, xyr_VelZ ) + &
      &          xyz_BuoyKm( xyz_PotTemp ) ) * MixLen ** 2.0d0 ))

    
  end subroutine Steady_Km_KW1978

!!!------------------------------------------------------------------------!!!
  subroutine Tendency_PotTemp_KW1978(Time, xyz_Km, xyz_Kh, xyz_PotTemp, xyz_Tend)
    !
    !温位の散逸加熱項
    !

    use xyz_deriv_module, only:  xyz_dx_pyz, pyz_dx_xyz, xyz_dy_xqz, xqz_dy_xyz, xyz_dz_xyr, xyr_dz_xyz

    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP),intent(in) :: Time
    real(DP),intent(in) :: xyz_Km(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数
    real(DP),intent(in) :: xyz_Kh(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数
    real(DP),intent(in) :: xyz_PotTemp(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量
    real(DP),intent(out):: xyz_Tend(imin:imax,jmin:jmax,kmin:kmax)
                                                    !
    real(DP)            :: xyz_PotTempAll(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: xyz_Turb(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: xyz_DispHeat(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流エネルギーの消散

    xyz_PotTempAll = xyz_PotTemp + xyz_PotTempBasicZ

    xyz_Turb =                                                   &
      &   xyz_dx_pyz( pyz_avr_xyz( xyz_Kh ) * pyz_dx_xyz( xyz_PotTempAll ) )  &
      & + xyz_dy_xqz( xqz_avr_xyz( xyz_Kh ) * xqz_dy_xyz( xyz_PotTempAll ) )  &
      & + xyz_dz_xyr( xyr_avr_xyz( xyz_Kh ) * xyr_dz_xyz( xyz_PotTempAll ) )

    xyz_DispHeat = (xyz_Km ** 3.0d0) &
      &        / (xyz_ExnerBasicZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0))

    xyz_Tend = xyz_Turb + xyz_DispHeat
    
    call HistoryAutoPut(Time, 'PotTempDisp', xyz_DispHeat(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(Time, 'PotTempTurb', xyz_Turb(1:nx, 1:ny, 1:nz))
    
  end subroutine Tendency_PotTemp_KW1978
  
!!!------------------------------------------------------------------------!!!
  subroutine Tendency_MixRt_KW1978(Time, xyza_Var, xyz_Kh, xyza_Tend)
    !
    ! x, z 方向に半格子ずれた点での乱流拡散
    !

    use xyz_deriv_module, only:  xyz_dx_pyz, pyz_dx_xyz, xyz_dy_xqz, xqz_dy_xyz, xyz_dz_xyr, xyr_dz_xyz
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP),intent(in) :: Time
    real(DP),intent(in) :: xyz_Kh(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP),intent(in) :: xyza_Var(imin:imax,jmin:jmax,kmin:kmax, nf)
                                                    !スカラー量
    real(DP),intent(out):: xyza_Tend(imin:imax,jmin:jmax,kmin:kmax, nf)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: xyza_MixRtAll(imin:imax,jmin:jmax,kmin:kmax, nf)
                                                    !スカラー量の水平乱流拡散
    integer             :: s

!    xyza_Turb = xyza_TurbScalar(xyza_MixRtBl + xyza_MixRtBasicZ, xyz_KhBl)  

    xyza_MixRtAll = xyza_MixRtAll + xyza_MixRtBasicZ
    
    do s = 1, nf
      xyza_Tend(:,:,:,s) = &
        &   xyz_dx_pyz( pyz_avr_xyz( xyz_Kh ) * pyz_dx_xyz( xyza_MixRtAll(:,:,:,s) ) )  &
        & + xyz_dy_xqz( xqz_avr_xyz( xyz_Kh ) * xqz_dy_xyz( xyza_MixRtAll(:,:,:,s) ) )  &
        & + xyz_dz_xyr( xyr_avr_xyz( xyz_Kh ) * xyr_dz_xyz( xyza_MixRtAll(:,:,:,s) ) )

      call HistoryAutoPut(Time, trim(SpcWetSymbol(s))//'_Turb', xyza_Tend(1:nx,1:ny,1:nz,s))
    end do
    
  end subroutine Tendency_MixRt_KW1978
  

!!!------------------------------------------------------------------------!!!
  subroutine Tendency_VelX_KW1978(pyz_VelX, xqz_VelY, xyr_VelZ, xyz_Km, pyz_TendVelX)
    !
    ! 水平速度に対する乱流拡散
    !

    use xyz_deriv_module, only:  pyz_dx_xyz, xyz_dx_pyz, pqz_dx_xqz, pyr_dx_xyr, &
      &                          pyz_dy_pqz, pqz_dy_pyz, &
      &                          pyz_dz_pyr, pyr_dz_pyz
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP),intent(in) :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平速度
    real(DP),intent(in) :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平速度
    real(DP),intent(in) :: xyr_VelZ(imin:imax,jmin:jmax,kmin:kmax)
                                                    !鉛直速度
    real(DP),intent(in) :: xyz_Km(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP), intent(out) :: pyz_TendVelX(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
      
    pyz_TendVelX = &
      &   2.0d0 * pyz_dx_xyz( xyz_Km * xyz_dx_pyz( pyz_VelX ) ) &
      & + pyz_dy_pqz(                                           &
      &       pqz_avr_xyz( xyz_Km ) * pqz_dx_xqz( xqz_VelY )    &
      &     + pqz_avr_xyz( xyz_Km ) * pqz_dy_pyz( pyz_VelX )    &
      &   )                                                     &
      & + pyz_dz_pyr(                                           &
      &       pyr_avr_xyz( xyz_Km ) * pyr_dx_xyr( xyr_VelZ )    &
      &     + pyr_avr_xyz( xyz_Km ) * pyr_dz_pyz( pyz_VelX )    &
      &   )                                                     &
      & - 2.0d0 * pyz_dx_xyz( ( xyz_Km ** 2.0d0 ) )             &
      &   / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )
    
  end subroutine Tendency_VelX_KW1978


!!!------------------------------------------------------------------------!!!
  subroutine Tendency_VelY_KW1978(pyz_VelX, xqz_VelY, xyr_VelZ, xyz_Km, xqz_TendVelY)
    !
    ! 水平速度に対する乱流拡散
    !

    use xyz_deriv_module, only: xqz_dx_pqz, pqz_dx_xqz, &
      &                         xqz_dy_xyz, xyz_dy_xqz, pqz_dy_pyz, xqr_dy_xyr, &
      &                         xqr_dz_xqz, xqz_dz_xqr
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP),intent(in) :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平速度
    real(DP),intent(in) :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平速度
    real(DP),intent(in) :: xyr_VelZ(imin:imax,jmin:jmax,kmin:kmax)
                                                    !鉛直速度
    real(DP),intent(in) :: xyz_Km(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP), intent(out) :: xqz_TendVelY(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
  
    xqz_TendVelY = &
      &   2.0d0 * xqz_dy_xyz( xyz_Km * xyz_dy_xqz( xqz_VelY ) ) &
      & + xqz_dx_pqz(                                           &
      &       pqz_avr_xyz( xyz_Km ) * pqz_dy_pyz( pyz_VelX )    &
      &     + pqz_avr_xyz( xyz_Km ) * pqz_dx_xqz( xqz_VelY )    &
      &   )                                                     &
      & + xqz_dz_xqr(                                           &
      &       xqr_avr_xyz( xyz_Km ) * xqr_dy_xyr( xyr_VelZ )    &
      &     + xqr_avr_xyz( xyz_Km ) * xqr_dz_xqz( xqz_VelY )    &
      &   )                                                     &
      & - 2.0d0 * xqz_dy_xyz( ( xyz_Km ** 2.0d0 ) )             &
      &   / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )
    
  end subroutine Tendency_VelY_KW1978

!!!------------------------------------------------------------------------!!!
  subroutine Tendency_VelZ_KW1978(pyz_VelX, xqz_VelY, xyr_VelZ, xyz_Km, xyr_TendVelZ)
    !
    !鉛直速度に対する乱流拡散
    !

    use xyz_deriv_module, only: xyr_dx_pyr, pyr_dx_xyr, &
      &                         xyr_dy_xqr, xqr_dy_xyr, &
      &                         xyr_dz_xyz, xyz_dz_xyr, pyr_dz_pyz, xqr_dz_xqz
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP),intent(in) :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平速度
    real(DP),intent(in) :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平速度
    real(DP),intent(in) :: xyr_VelZ(imin:imax,jmin:jmax,kmin:kmax)
                                                    !鉛直速度
    real(DP),intent(in) :: xyz_Km(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP),intent(out):: xyr_TendVelZ(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
  
    xyr_TendVelZ = &
      & + 2.0d0 * xyr_dz_xyz( xyz_Km * xyz_dz_xyr( xyr_VelZ ) )  &
      & + xyr_dx_pyr(                                            &
      &      pyr_avr_xyz( xyz_Km ) * pyr_dz_pyz( pyz_VelX )      &
      &    + pyr_avr_xyz( xyz_Km ) * pyr_dx_xyr( xyr_VelZ )      &
      &   )                                                      &
      & + xyr_dy_xqr(                                            &
      &      xqr_avr_xyz( xyz_Km ) * xqr_dz_xqz( xqz_VelY )      &
      &    + xqr_avr_xyz( xyz_Km ) * xqr_dy_xyr( xyr_VelZ )      &
      &   )                                                      &
      & - 2.0d0 * xyr_dz_xyz( ( xyz_Km ** 2.0d0 ) )              &
      &   / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )
    
  end subroutine Tendency_VelZ_KW1978
  
!!!------------------------------------------------------------------------!!!
  function xyz_ShearKm(xyz_Km, pyz_VelX, xqz_VelY, xyr_VelZ)
    !
    !速度シアーによる乱流エネルギー生成項
    !

    use xyz_deriv_module, only: xyz_dx_pyz, pyr_dx_xyr, pqz_dx_xqz, &
      &                         xyz_dy_xqz, xqr_dy_xyr, pqz_dy_pyz, &
      &                         xyz_dz_xyr, pyr_dz_pyz, xqr_dz_xqz
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP)            :: xyz_ShearKm(imin:imax,jmin:jmax,kmin:kmax)
                                                        !渦粘性係数の移流
    real(DP),intent(in) :: xyz_Km(imin:imax,jmin:jmax,kmin:kmax)
                                                        !渦粘性係数
    real(DP),intent(in) :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                                                        !水平速度
    real(DP),intent(in) :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                                                        !水平速度
    real(DP),intent(in) :: xyr_VelZ(imin:imax,jmin:jmax,kmin:kmax)
                                                        !鉛直速度
  
!    xyz_ShearKm = 0.0d0
    xyz_ShearKm = &
      &   ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 )                    &
      & * (                                                        &
      &      ( xyz_dx_pyz( pyz_VelX ) ) ** 2.0d0                   &
      &    + ( xyz_dy_xqz( xqz_VelY ) ) ** 2.0d0                   &
      &    + ( xyz_dz_xyr( xyr_VelZ ) ) ** 2.0d0                   &
      &    + 5.0d-1                                                &
      &      * (                                                   &
      &          (                                                 &
      &              xyz_avr_pyr( pyr_dz_pyz( pyz_VelX ) )         &
      &            + xyz_avr_pyr( pyr_dx_xyr( xyr_VelZ ) )         &
      &           ) ** 2.0d0                                       &
      &        + (                                                 &
      &              xyz_avr_xqr( xqr_dy_xyr( xyr_VelZ ) )         &
      &            + xyz_avr_xqr( xqr_dz_xqz( xqz_VelY ) )         &
      &           ) ** 2.0d0                                       &
      &        + (                                                 &
      &              xyz_avr_pqz( pqz_dx_xqz( xqz_VelY ) )         &
      &            + xyz_avr_pqz( pqz_dy_pyz( pyz_VelX ) )         &
      &           ) ** 2.0d0                                       &
      &        )                                                   &
      &   )                                                        &
      & - xyz_Km * (  xyz_dx_pyz( pyz_VelX ) &
      &             + xyz_dy_xqz( xqz_VelY ) &
      &             + xyz_dz_xyr( xyr_VelZ ) ) / 3.0d0
   
  end function xyz_ShearKm


!!!------------------------------------------------------------------------!!!
  function xyz_TurbKm(xyz_Km)
    !
    ! 乱流エネルギーの拡散項
    !

    use xyz_deriv_module, only: xyz_dx_pyz, pyz_dx_xyz, &
      &                         xyz_dy_xqz, xqz_dy_xyz, &
      &                         xyz_dz_xyr, xyr_dz_xyz 

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP)            :: xyz_TurbKm(imin:imax,jmin:jmax,kmin:kmax)
                                               ! 乱流エネルギーの拡散項
    real(DP),intent(in) :: xyz_Km(imin:imax,jmin:jmax,kmin:kmax)
                                               !渦粘性係数

!    xz_TurbKm = 0.0d0
    xyz_TurbKm =                                       &
      &   xyz_dx_pyz( pyz_dx_xyz( xyz_Km ** 2.0d0 ) )  &
      & + xyz_dy_xqz( xqz_dy_xyz( xyz_Km ** 2.0d0 ) )  &
      & + xyz_dz_xyr( xyr_dz_xyz( xyz_Km ** 2.0d0 ) )  &
      & + (xyz_avr_pyz( pyz_dx_xyz( xyz_Km ) ) ) ** 2.0d0 & 
      & + (xyz_avr_xqz( xqz_dy_xyz( xyz_Km ) ) ) ** 2.0d0 & 
      & + (xyz_avr_xyr( xyr_dz_xyz( xyz_Km ) ) ) ** 2.0d0
    
  end function xyz_TurbKm


!!!------------------------------------------------------------------------!!!
  function xyz_DispKm(xyz_Km)
    !
    !乱流エネルギーの消散項
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP)            :: xyz_DispKm(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流エネルギーの消散
    real(DP),intent(in) :: xyz_Km(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数

!    xz_DispKm = 0.0d0
    xyz_DispKm = - (xyz_Km ** 2.0d0) * 5.0d-1 / (MixLen ** 2.0d0)
    
  end function xyz_DispKm

!!!------------------------------------------------------------------------!!!
  function xyz_BuoyKm(xyz_PotTemp)
    !
    !乱流エネルギーの浮力項を計算. 乾燥大気版. 
    !

    use xyz_deriv_module, only: xyr_dz_xyz
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP),intent(in)  :: xyz_PotTemp(imin:imax,jmin:jmax,kmin:kmax)
                                               !温位擾乱
    real(DP)             :: xyz_BuoyKm(imin:imax,jmin:jmax,kmin:kmax)
                                               !浮力項

    !浮力項の計算
!    xz_BuoyKm = 0.0d0
    xyz_BuoyKm = &
      &  - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) &
      &       * xyz_avr_xyr( xyr_dz_xyz( xyz_PotTemp + xyz_PotTempBasicZ ) ) &
      &       / ( 2.0d0 * xyz_PotTempBasicZ )
    
  end function xyz_BuoyKm

!!!------------------------------------------------------------------------!!!
  function xyz_BuoyMoistKm(xyz_PotTemp, xyz_Exner, xyza_MixRt)
    !

    use xyz_deriv_module, only: xyr_dz_xyz

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: xyz_PotTemp(imin:imax,jmin:jmax,kmin:kmax)
                                               !温位
    real(DP), intent(in) :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)
                                               !無次元圧力
    real(DP), intent(in) :: xyza_MixRt(imin:imax,jmin:jmax,kmin:kmax, nf)
                                               !凝縮成分の混合比
    real(DP) :: xyza_MixRtAll(imin:imax,jmin:jmax,kmin:kmax, nf)
                                               !凝縮成分の混合比
    real(DP) :: xyza_MixRtAll2(imin:imax,jmin:jmax,kmin:kmax, nf)
                                               !凝縮成分の混合比
    real(DP) :: xyz_BuoyMoistKm(imin:imax,jmin:jmax,kmin:kmax)
                                               !
    real(DP) :: xyza_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_EffPotTemp(imin:imax,jmin:jmax,kmin:kmax)    
                                               !浮力に対する温度差の寄与
    real(DP) :: xyz_EffMolWt(imin:imax,jmin:jmax,kmin:kmax)    
                                               !浮力に対する分子量差の寄与
    real(DP) :: xyza_MixRtPerMolWt(imin:imax,jmin:jmax,kmin:kmax, GasNum)
                                               !混合比/分子量
    integer              :: s
    
    
    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xyz_TempAll = ( xyz_PotTemp + xyz_PotTempBasicZ ) * ( xyz_Exner + xyz_ExnerBasicZ )
    xyza_MixRtAll = xyza_MixRtBasicZ + xyza_MixRt
    xyza_LatentHeat = 0.0d0
    
    !作業配列の初期化. 気体のみ利用
    do s = 1, GasNum
      xyza_MixRtPerMolWt(:,:,:,s) = xyza_MixRt(:,:,:,IdxG(s)) / MolWtWet(IdxG(s))
    end do
    
    !温度の効果
    xyz_EffPotTemp = xyz_PotTemp / xyz_PotTempBasicZ 
    
    !分子量効果 + 引きづりの効果
    xyz_EffMolWt =                                           &
      & + sum(xyza_MixRtPerMolWt, 4)                         &
      &    / ( 1.0d0 / MolWtDry + xyz_MixRtBasicZPerMolWt )  &
      & - sum(xyza_MixRt, 4) / ( 1.0d0 + xyz_MixRtBasicZ )
    
    !蒸気が蒸発する場合の潜熱を計算
    !  分子量の部分はいつでも効くが潜熱は飽和していないと効かないので, 
    !  雲の混合比がゼロの時には, 潜熱の寄与はゼロとなるように調節している
    xyza_MixRtAll2 = xyza_MixRtAll
!    xzya_MixRtAll2(:,:,:,IdxNH3) = &
!      &  xyza_MixRtAll(:,:,:,IdxNH3) - xyza_MixRtAll(:,:,:,IdxH2S) 

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

    xyz_EffHeat = ( sum( xyza_LatentHeat, 4 ) * xyz_EffMolWtBasicZ &
!      &            + ReactHeatNH4SH * &
!      &            (xyza_MixRtAll(:,:,:,IdxH2S) + xyza_MixRtAll(:,:,:,IdxH2S)) &
      &          ) / ( CpDry * xyz_ExnerBasicZ ) 
    
    !乱流拡散係数の時間発展式の浮力項を決める
    xyz_BuoyMoistKm = &
      &  - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 )     &
      &    * xyz_avr_xyr(                                           &
      &        xyr_dz_xyz(                                          &
      &           xyz_EffHeat                                       &
      &         + xyz_PotTempBasicZ / xyz_EffMolWtBasicZ            &
      &           * ( 1.0d0 + xyz_EffPotTemp + xyz_EffMolWt )       &
      &         )                                                   &
      &      )                                                      &
      &    / ( 2.0d0 * xyz_PotTempBasicZ / xyz_EffMolWtBasicZ)                   

  end function xyz_BuoyMoistKm


!!!-------------------------------------------------------------------!!!
  function xyz_NumDiffKm(xyz_Scalar)
    !
    ! x, z 方向に半格子ずれた点での数値拡散項を評価
    !

    use xyz_deriv_module, only: pyz_dx_xyz, xyz_dx_pyz, xqz_dy_xyz, xyz_dy_xqz, &
      &                         xyr_dz_xyz, xyz_dz_xyr
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xyz_Scalar(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量
    real(8)             :: xyz_NumDiffKm(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平方向の数値拡散
    
    xyz_NumDiffKm =   &
      &    NuH * (xyz_dx_pyz(pyz_dx_xyz( xyz_Scalar ))) &
      &  + NuH * (xyz_dy_xqz(xqz_dy_xyz( xyz_Scalar ))) &
      &  + NuV * (xyz_dz_xyr(xyr_dz_xyz( xyz_Scalar ))) 
    
  end function xyz_NumDiffKm

end module Turbulence_KW1978


!      call HistoryAutoAddVariable(  &
!        & varname=trim(SpcWetSymbol(l))//'_Fall', & 
!        & dims=(/'x','y','z','t'/),     &
!        & longname='Fall Rain term of '          &
!        &           //trim(SpcWetSymbol(l))//' mixing ratio',  &
!        & units='kg.kg-1.s-1',    &
!        & xtype='double')
