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

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

  !モジュール読み込み 
  use gridset, only:  DimXMin,           &! x 方向の配列の下限
    &                 DimXMax,           &! x 方向の配列の上限
    &                 DimZMin,           &! z 方向の配列の下限
    &                 DimZMax,           &! z 方向の配列の上限
    &                 SpcNum,            &! z 方向の配列の上限
    &                 DelX,              &! x 方向の格子点間隔
    &                 DelZ                ! z 方向の格子点間隔
  use basicset, only: CpDry,             &!乾燥成分の比熱
    &                 Grav,              &!重力加速度
    &                 xz_PotTempBasicZ,  &!基本場の温位
    &                 xz_EffMolWtBasicZ, &!基本場の分子量の効果
    &                 xz_ExnerBasicZ      !エクスナー関数の基本場
  use average,  only: xz_avr_pz, xz_avr_xr, xz_avr_pr, &
    &                 pz_avr_xz, &
    &                 pr_avr_xz, &
    &                 xr_avr_xz
!  use differentiate_center2, only: xz_dx_pz, xz_dz_xr, &
!    &                              xr_dx_pr, xr_dz_xz, &
!    &                              pz_dx_xz, pz_dz_pr, &
!    &                              pr_dx_xr, pr_dz_pz
  use differentiate_center4, only: xz_dx_pz, xz_dz_xr, &
    &                              xr_dx_pr, xr_dz_xz, &
    &                              pz_dx_xz, pz_dz_pr, &
    &                              pr_dx_xr, pr_dz_pz
  use StorePotTemp, only: StorePotTempDisp, StorePotTempTurb
  use StoreMixRt,   only: StoreMixRtTurb


  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public Turbulence_Init
  public xz_TurbScalar
  public xza_TurbScalar
  public pz_TurbVelX
  public xr_TurbVelZ
  public xz_ShearKm
  public xz_DispKm
  public xz_DispHeat
  public xz_BuoyKm
  public EddyViscosity


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

  !値を保存
  save Cm, MixLen

contains

!!!------------------------------------------------------------------------!!!
  subroutine turbulence_init()
    !
    ! Turbulence モジュールの初期化ルーチン
    ! 
    
    !暗黙の型宣言禁止
    implicit none

    !混合距離
!    MixLen = sqrt(DelX * DelZ) 
    MixLen = min(DelX,  DelZ) 

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

    !変数定義
    real(8), intent(in) :: xz_Var(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !スカラー量
    real(8), intent(in) :: xz_Kh(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !乱流拡散係数
    real(8)             :: xz_TurbScalar(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !スカラー量の水平乱流拡散
  
!    xz_TurbScalar = 0.0d0  !初期化
    xz_TurbScalar =                                          &
      &   xz_dx_pz(pz_avr_xz(xz_Kh) * pz_dx_xz(xz_Var))  &
      & + xz_dz_xr(xr_avr_xz(xz_Kh) * xr_dz_xz(xz_Var))

    call StorePotTempTurb( xz_TurbScalar )    
    
  end function xz_TurbScalar


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

    !変数定義
    real(8), intent(in) :: xza_Var(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                    !スカラー量
    real(8), intent(in) :: xz_Kh(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !乱流拡散係数
    real(8)             :: xza_TurbScalar(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                    !スカラー量の水平乱流拡散
    integer             :: s
  
    do s = 1, SpcNum
      xza_TurbScalar(:,:,s) =                                      &
        &   xz_dx_pz(pz_avr_xz(xz_Kh) * pz_dx_xz(xza_Var(:,:,s)))  &
        & + xz_dz_xr(xr_avr_xz(xz_Kh) * xr_dz_xz(xza_Var(:,:,s)))
    end do

    call StoreMixRtTurb( xza_TurbScalar )    
    
  end function xza_TurbScalar


!!!------------------------------------------------------------------------!!!
  function pz_TurbVelX(xz_Km, pz_VelX, xr_VelZ)
    !
    ! 水平速度に対する乱流拡散
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in) :: pz_VelX(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !水平速度
    real(8), intent(in) :: xr_VelZ(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !鉛直速度
    real(8), intent(in) :: xz_Km(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !乱流拡散係数
    real(8)             :: pz_TurbVelX(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !スカラー量の水平乱流拡散
  
!    pz_TurbVelX = 0.0d0  !初期化  
    pz_TurbVelX = &
      &   2.0d0 * pz_dx_xz( xz_Km * xz_dx_pz( pz_VelX ) )     &
      & + pz_dz_pr(                                           &
      &       pr_avr_xz( xz_Km ) * pr_dx_xr( xr_VelZ )        &
      &     + pr_avr_xz( xz_Km ) * pr_dz_pz( pz_VelX )        &
      &   )                                                   &
      & - 2.0d0 * pz_dx_xz( ( xz_Km ** 2.0d0 ) )              &
      &   / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )
    
  end function pz_TurbVelX


!!!------------------------------------------------------------------------!!!
  function xr_TurbVelZ(xz_Km, pz_VelX, xr_VelZ)
    !
    !鉛直速度に対する乱流拡散
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in) :: pz_VelX(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !水平速度
    real(8), intent(in) :: xr_VelZ(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !鉛直速度
    real(8), intent(in) :: xz_Km(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !乱流拡散係数
    real(8)             :: xr_TurbVelZ(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !スカラー量の水平乱流拡散
  
!    xr_TurbVelZ = 0.0d0  !初期化
    xr_TurbVelZ = &
      & + 2.0d0 * xr_dz_xz( xz_Km * xz_dz_xr( xr_VelZ ) )        &
      & + xr_dx_pr(                                              &
      &      pr_avr_xz( xz_Km ) * pr_dz_pz( pz_VelX )            &
      &    + pr_avr_xz( xz_Km ) * pr_dx_xr( xr_VelZ )            &
      &   )                                                      &
      & - 2.0d0 * xr_dz_xz( ( xz_Km ** 2.0d0 ) )                 &
      &   / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )
    
  end function xr_TurbVelZ


!!!------------------------------------------------------------------------!!!
  function xz_ShearKm(xz_Km, pz_VelX, xr_VelZ)
    !
    !速度シアーによる乱流エネルギー生成項
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8)             :: xz_ShearKm(DimXMin:DimXMax, DimZMin:DimZMax)
                                                        !渦粘性係数の移流
    real(8), intent(in) :: xz_Km(DimXMin:DimXMax, DimZMin:DimZMax)
                                                        !渦粘性係数
    real(8), intent(in) :: pz_VelX(DimXMin:DimXMax, DimZMin:DimZMax)
                                                        !水平速度
    real(8), intent(in) :: xr_VelZ(DimXMin:DimXMax, DimZMin:DimZMax)
                                                        !鉛直速度
  
!    xz_ShearKm = 0.0d0
    xz_ShearKm = &
      &   ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 )                    &
      & * (                                                        &
      &      (xz_dx_pz(pz_VelX)) ** 2.0d0                          &
      &    + (xz_dz_xr(xr_VelZ)) ** 2.0d0                          &
      &    + 5.0d-1                                                &
      &      * (                                                   &
      &          (                                                 &
      &              xz_avr_pr(pr_dz_pz(pz_VelX))                  &
      &            + xz_avr_pr(pr_dx_xr(xr_VelZ))                  &
      &           ) ** 2.0d0                                       &
      &        )                                                   &
      &   )                                                        &
      & - xz_Km * (xz_dx_pz(pz_VelX) + xz_dz_xr(xr_VelZ)) / 3.0d0
    
  end function xz_ShearKm


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

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

!    xz_TurbKm = 0.0d0
    xz_TurbKm =                                      &
!      &     xz_dx_pz(pz_dx_xz(xz_Km)) * 5.0d-1       &
!      &   + xz_dz_xr(xr_dz_xz(xz_Km)) * 5.0d-1       &
      &     xz_dx_pz(pz_dx_xz(xz_Km ** 2.0d0)) * 5.0d-1  &
      &   + xz_dz_xr(xr_dz_xz(xz_Km ** 2.0d0)) * 5.0d-1  &
      &   + (xz_avr_pz(pz_dx_xz(xz_Km))) ** 2.0d0        & 
      &   + (xz_avr_xr(xr_dz_xz(xz_Km))) ** 2.0d0
    
  end function xz_TurbKm


!!!------------------------------------------------------------------------!!!
  function xz_DispKm(xz_Km)
    !
    !乱流エネルギーの消散項
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8)             :: xz_DispKm(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !乱流エネルギーの消散
    real(8), intent(in) :: xz_Km(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !渦粘性係数

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


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

!    xz_DispHeat = 0.0d0
    xz_DispHeat = (xz_Km ** 3.0d0) * xz_EffMolWtBasicZ &
      &        / (xz_ExnerBasicZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0))

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


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

    !浮力項の計算
!    xz_BuoyKm = 0.0d0
    xz_BuoyKm = &
      &  - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) &
      &       * xz_avr_xr( xr_dz_xz( xz_PotTemp + xz_PotTempBasicZ ) ) &
      &       / ( 2.0d0 * xz_PotTempBasicZ )
    
  end function xz_BuoyKm

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

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

    real(8), intent(in)  :: pz_VelX(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8), intent(in)  :: xr_VelZ(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8), intent(in)  :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in)  :: xz_Km(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8), intent(out) :: xz_KmA(DimXMin:DimXMax, DimZMin:DimZMax) 

    xz_KmA = &
      & sqrt ( &
      &   max( 0.0d0, &
      &        ( xz_ShearKm(xz_Km, pz_VelX, xr_VelZ) + &
      &          xz_BuoyKm(xz_PotTemp) ) * MixLen ** 2.0d0 ) )

  end subroutine EddyViscosity

end module Turbulence
