!= Module Turbulence_3d
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA, Masatsugu
! Version::   $Id: turbulence_steady.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

  use gridset_3d, 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, only: xyz_PotTempBasicZ, &!基本場の温位
    &                 xyz_ExnerBasicZ     !エクスナー関数の基本場
  use xyz_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  

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public Turbulence_Init
  public EddyViscosity

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

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_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_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


end module Turbulence_KW1978
