!= Module Turbulence_3d
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA, Masatsugu
! Version::   $Id: turbulence_3d.f90,v 1.3 2008-06-04 04:51:28 odakker 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_3d
  !
  !モデルの物理過程を計算するために必要となる関数群を束ねたモジュール
  !具体的には以下の項を計算するための関数を格納する.  
  !  * 乱流拡散項
  !  * 乱流エネルギーの時間発展方程式に含まれる各項
  !  * 散逸加熱項
  !

  !モジュール読み込み 
  use dc_types, only : DP

  use gridset, only:  DimXMin,           &! x 方向の配列の下限
    &                 DimXMax,           &! x 方向の配列の上限
    &                 DimYMin,           &! y 方向の配列の下限
    &                 DimYMax,           &! y 方向の配列の上限
    &                 DimZMin,           &! z 方向の配列の下限
    &                 DimZMax,           &! z 方向の配列の上限
    &                 SpcNum              ! z 方向の配列の上限
  use axesset, only:  x_dx,              &! x 方向の格子点間隔
    &                 y_dy,              &! y 方向の格子点間隔
    &                 z_dz                ! z 方向の格子点間隔
  use constants, only: CpDry,          &!乾燥成分の比熱
    &                  Grav                !重力加速度
  use basicset_3d, only: xyz_PotTempBasicZ, &!基本場の温位
    &                    xyz_ExnerBasicZ     !エクスナー関数の基本場
  use xyz_base_module, 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 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 StorePotTemp_3d, only: StorePotTempDisp, StorePotTempTurb
!  use StoreMixRt,   only: StoreMixRtTurb

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public Turbulence_Init
  public xyz_TurbScalar
  public xyza_TurbScalar
  public pyz_TurbVelX
  public xqz_TurbVelY
  public xyr_TurbVelZ
  public xyz_ShearKm
  public xyz_DispKm
  public xyz_DispHeat
  public xyz_BuoyKm
  public EddyViscosity
  public xyz_TurbKm


  !変数定義
  real(DP) :: Cm     = 2.0d-1    !乱流エネルギー診断式の係数 
  real(DP) :: MixLen = 0.0d0     !平均混合距離

  !値を保存
  save Cm, MixLen

contains

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

    !暗黙の型宣言禁止
    implicit none

    real(DP):: MinDelX
    real(DP):: MinDelY
    real(DP):: MinDelZ

    !混合距離
    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 ) )

  end subroutine turbulence_init
  
!!!------------------------------------------------------------------------!!!
  function xyz_TurbScalar(xyz_Var, xyz_Kh)
    !
    ! x, y, z 方向に半格子ずれた点での乱流拡散
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP),intent(in) :: xyz_Var(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !スカラー量
    real(DP),intent(in) :: xyz_Kh(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !乱流拡散係数
    real(DP)            :: xyz_TurbScalar(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !スカラー量の水平乱流拡散
  
!    xyz_TurbScalar = 0.0d0  !初期化
    xyz_TurbScalar =                                                   &
      &   xyz_dx_pyz( pyz_avr_xyz( xyz_Kh ) * pyz_dx_xyz( xyz_Var ) )  &
      & + xyz_dy_xqz( xqz_avr_xyz( xyz_Kh ) * xqz_dy_xyz( xyz_Var ) )  &
      & + xyz_dz_xyr( xyr_avr_xyz( xyz_Kh ) * xyr_dz_xyz( xyz_Var ) )

!    call StorePotTempTurb( xyz_TurbScalar )    
    
  end function xyz_TurbScalar


!!!------------------------------------------------------------------------!!!
  function xyza_TurbScalar(xyza_Var, xyz_Kh)
    !
    ! x, z 方向に半格子ずれた点での乱流拡散
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP),intent(in) :: xyza_Var(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                                    !スカラー量
    real(DP),intent(in) :: xyz_Kh(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !乱流拡散係数
    real(DP)            :: xyza_TurbScalar(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                                    !スカラー量の水平乱流拡散
    integer             :: s
  
    do s = 1, SpcNum
      xyza_TurbScalar(:,:,:,s) = xyz_TurbScalar( xyza_Var(:,:,:,s), xyz_Kh )  
    end do

!    call StoreMixRtTurb( xyza_TurbScalar )    
    
  end function xyza_TurbScalar


!!!------------------------------------------------------------------------!!!
  function pyz_TurbVelX(xyz_Km, pyz_VelX, xqz_VelY, xyr_VelZ)
    !
    ! 水平速度に対する乱流拡散
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP),intent(in) :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !水平速度
    real(DP),intent(in) :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !水平速度
    real(DP),intent(in) :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !鉛直速度
    real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !乱流拡散係数
    real(DP)            :: pyz_TurbVelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !スカラー量の水平乱流拡散
  
!    pyz_TurbVelX = 0.0d0  !初期化  
    pyz_TurbVelX = &
      &   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 function pyz_TurbVelX


!!!------------------------------------------------------------------------!!!
  function xqz_TurbVelY(xyz_Km, pyz_VelX, xqz_VelY, xyr_VelZ)
    !
    ! 水平速度に対する乱流拡散
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP),intent(in) :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !水平速度
    real(DP),intent(in) :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !水平速度
    real(DP),intent(in) :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !鉛直速度
    real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !乱流拡散係数
    real(DP)            :: xqz_TurbVelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !スカラー量の水平乱流拡散
  
!    xqz_TurbVelY = 0.0d0  !初期化  
    xqz_TurbVelY = &
      &   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 function xqz_TurbVelY


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

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP),intent(in) :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !水平速度
    real(DP),intent(in) :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !水平速度
    real(DP),intent(in) :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !鉛直速度
    real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !乱流拡散係数
    real(DP)            :: xyr_TurbVelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !スカラー量の水平乱流拡散
  
!    xyr_TurbVelZ = 0.0d0  !初期化
    xyr_TurbVelZ = &
      & + 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 function xyr_TurbVelZ

!!!------------------------------------------------------------------------!!!
  function xyz_ShearKm(xyz_Km, pyz_VelX, xqz_VelY, xyr_VelZ)
    !
    !速度シアーによる乱流エネルギー生成項
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP)            :: xyz_ShearKm(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                        !渦粘性係数の移流
    real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                        !渦粘性係数
    real(DP),intent(in) :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                        !水平速度
    real(DP),intent(in) :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                        !水平速度
    real(DP),intent(in) :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                        !鉛直速度
  
!    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)
    !
    ! 乱流エネルギーの拡散項
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP)            :: xyz_TurbKm(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                               ! 乱流エネルギーの拡散項
    real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                               !渦粘性係数

!    xz_TurbKm = 0.0d0
    xyz_TurbKm =                                               &
      &   xyz_dx_pyz( pyz_dx_xyz( xyz_Km ** 2.0d0 ) ) * 5.0d-1 &
      & + xyz_dy_xqz( xqz_dy_xyz( xyz_Km ** 2.0d0 ) ) * 5.0d-1 &
      & + xyz_dz_xyr( xyr_dz_xyz( xyz_Km ** 2.0d0 ) ) * 5.0d-1 &
      & + (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(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !乱流エネルギーの消散
    real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !渦粘性係数

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


!!!------------------------------------------------------------------------!!!
  function xyz_DispHeat(xyz_Km)
    !
    !温位の散逸加熱項
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP)            :: xyz_DispHeat(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !乱流エネルギーの消散
    real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !渦粘性係数

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

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

    !値の保管
!    call StorePotTempDisp( xyz_DispHeat )
    
  end function xyz_DispHeat


!!!------------------------------------------------------------------------!!!
  function xyz_BuoyKm(xyz_PotTemp)
    !
    !乱流エネルギーの浮力項を計算. 乾燥大気版. 
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP),intent(in)  :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                               !温位擾乱
    real(DP)             :: xyz_BuoyKm(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                               !浮力項

    !浮力項の計算
!    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

!!!------------------------------------------------------------------------!!!
  subroutine EddyViscosity(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(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) 
    real(DP),intent(in)  :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) 
    real(DP),intent(in)  :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) 
    real(DP),intent(in)  :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
    real(DP),intent(in)  :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) 
    real(DP),intent(out) :: xyz_KmA(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) 

    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 EddyViscosity

end module Turbulence_3d
