!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine EddyViscosity
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: eddyviscosity.f90,v 1.2.4.1 2006/01/13 06:55:04 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!1.5 次の乱流クロージャーを利用し, 
!次の時刻の渦粘性係数と渦拡散係数を求めるサブルーチン.
!
!== Error Handling
!
!== Known Bugs
!
!  * 2005/09/05 現在, 空気中の乾燥成分と凝縮性分の分子量を同じにしてある
!  * 浮力項の分子量効果を考慮していない
!
!== Note
!
!  * 定式化は CReSS のマニュアルを参照した
!  * 移流項の計算には ***_adv となる変数を利用する. 
!  * 生成項, 拡散項, 消散項の計算には ***_dif となる変数を利用する. 
!
!== Future Plans
!
                                                                 !=end
subroutine EddyViscosity(                         &
     & ss_Km_in,                                  &
     & DelTime,                                   & 
     & fs_VelX_adv, sf_VelZ_adv, ss_Km_adv,       & !移流
     & fs_VelX_dif, sf_VelZ_dif, ss_PotTemp_dif,  & !生成, 拡散, 消散
     & ss_Km_dif, ss_Exner_dif,                   &
     & ss_MixRtV_dif, ss_MixRtC_dif,    &
     & ss_MixRtR_dif,                         &
     & ss_Km_out,   ss_Kh_out )
                                                                 !=begin
  !== Dependency
  use dc_trace, only: BeginSub, EndSub
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, MixLen
  use bcset,    only: ss_BC
  use physset,  only: Grav, CpDry, GasRDry, LatentHeat, epsilon
  use basicset, only: ss_PotTempBasicZ, ss_ExnerBasicZ
  use arareset, only: Cm, Min
  use average,  only: ss_avr_fs, ss_avr_ff, ss_avr_sf
  use differentiate_center4, only: fs_dx_ss, ss_dx_fs, ff_dx_sf, &
       &                   ss_dz_sf, sf_dz_ss, ff_dz_fs
                                                                 !=end
  !== 暗黙の型宣言禁止
  implicit none
                                                                 !=begin
  !== Input
  real(8), intent(in)  :: DelTime
  real(8), intent(in)  :: ss_Km_in(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: fs_VelX_adv(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: sf_VelZ_adv(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_Km_adv(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_Exner_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: fs_VelX_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: sf_VelZ_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_PotTemp_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_Km_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_MixRtV_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_MixRtC_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_MixRtR_dif(DimXMin:DimXMax, DimZMin:DimZMax)

  !== Output
  real(8), intent(out) :: ss_Km_out(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(out) :: ss_Kh_out(DimXMin:DimXMax, DimZMin:DimZMax)
                                                                 !=end
  !== Work
  real(8)  :: ss_AdvKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_BuoyKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_ShearKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_DiffKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_DispKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_TendKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_Temp_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_AA(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_MixRatSum_dif(DimXMin:DimXMax, DimZMin:DimZMax)

 
  call BeginSub("EddyViscosity", &
       &        fmt="%c",    &
       &        c1="calculation of eddy viscosity whith 2.5 order closure.")


  !==== 変数初期化
  ss_AdvKm   = 0.0d0
  ss_BuoyKm  = 0.0d0  
  ss_ShearKm = 0.0d0
  ss_DiffKm  = 0.0d0
  ss_DispKm  = 0.0d0
  ss_TendKm  = 0.0d0
  ss_Temp_dif = 0.0d0
  ss_AA       = 0.0d0
  ss_MixRatSum_dif = 0.0d0


  !=== スカラー量の移流. 中心差分. 
  ss_AdvKm = &
       &   ss_avr_fs( fs_VelX_adv * fs_dx_ss( ss_Km_adv ) ) &
       & + ss_avr_sf( sf_VelZ_adv * sf_dz_ss( ss_Km_adv ) ) 


  !=== 浮力による乱流エネルギー生成
  !----- 本当なら, 温位を仮温位に読み替える必要があるが, 現在は
  !----- 乾燥成分と凝縮成分の分子量を一定としているので, その必要はない
  ss_BuoyKm = &
    &    3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) &
    &       * ss_avr_sf( sf_dz_ss( ss_PotTemp_dif + ss_PotTempBasicZ ) )&
    &       / ( 2.0d0 * ss_PotTempBasicZ )

  !=== 速度シアーによる乱流エネルギー生成
  ss_ShearKm = &
    &  ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 )               &
    &    * (                                               &
    &          ( ss_dx_fs( fs_VelX_dif ) ) ** 2.0d0        &
    &       + ( ss_dz_sf( sf_VelZ_dif ) ) ** 2.0d0         &
    &       + ( (   ss_avr_ff( ff_dz_fs( fs_VelX_dif ) )   &
    &             + ss_avr_ff( ff_dx_sf( sf_VelZ_dif ) )   &
    &           ) ** 2.0d0 ) / 2.0d0                       &
    &      )                                               &
    & - ss_Km_dif                                          &
    &    * ( ss_dx_fs( fs_VelX_dif ) + ss_dz_sf( sf_VelZ_dif ) ) / 3.0d0

  
  !=== 乱流エネルギーの拡散
  ss_DiffKm =                                                 &
    &     ss_dx_fs( fs_dx_ss( ss_Km_dif ) ) / 2.0d0        &
    &   + ss_dz_sf( sf_dz_ss( ss_Km_dif ) ) / 2.0d0        &
    &   + ( ss_avr_fs( fs_dx_ss( ss_Km_dif ) ) ) ** 2.0d0  & 
    &   + ( ss_avr_sf( sf_dz_ss( ss_Km_dif ) ) ) ** 2.0d0


  !=== 乱流エネルギーの消散
  ss_DispKm = ( ss_Km_dif ** 2.0d0 ) / ( 2.0d0 * ( MixLen ** 2.0d0 ) )
  
  
  !=== 項をまとめる
  ss_TendKm = - ss_AdvKm - ss_BuoyKm + ss_ShearKm + ss_DiffKm - ss_DispKm 

  
  !=== 時間積分
  ss_Km_out = ss_Km_in + DelTime * ss_TendKm

  !=== 負の処理
  where (ss_Km_out < 0.0d0)
    ss_Km_out = 0.0d0
  end where
  
  !=== 境界条件
  call boundary( ss_BC, ss_Km_out )
    
  !=== 渦拡散定数を決める
  ss_Kh_out = 3.0d0 * ss_Km_out

  
  call EndSub("EddyViscosity")
  
end subroutine EddyViscosity
