!= Module Turbulence
!
! Authors::   SUGIYAMA Ko-ichiro
! Version::   $Id: turbulence.f90,v 1.10 2006/11/22 15:40:54 sugiyama Exp $
! Tag Name::  $Name: arare4-20061224 $
! 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, pz_avr_pr, &
    &                 pr_avr_xr, pr_avr_pz, pr_avr_xz, &
    &                 xr_avr_pr, 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 StoreSet, only: StoreDisp, StoreTurb
  use StoreSet2, only: Store2Turb


  !暗黙の型宣言禁止
  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

  !変数定義
  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 StoreTurb( 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 Store2Turb( 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_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 StoreDisp( 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
  
end module Turbulence
