!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine EddyViscosity
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: eddyviscosity.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!2.5 次の乱流クロージャーを利用し, 
!次の時刻の渦粘性係数と渦拡散係数を求めるサブルーチン.
!
!== Error Handling
!
!== Known Bugs
!
!== 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_Km_out, ss_Kh_out )
                                                                 !=begin
  !== Dependency
  use dc_trace, only: BeginSub, EndSub
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, DelXZ
  use bcset,    only: ss_BC
  use physset,  only: Grav
  use basicset, only: ss_PotTempBasicZ
  use arareset, only: Cm
  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)  :: 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)

  !== 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)
  

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


  !=== スカラー量の移流. 中心差分. 
  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 ) * ( DelXZ ** 2.0d0 ) &
       &       * ss_avr_sf( sf_dz_ss( ss_PotTemp_dif ) ) &
       &       / ( 2.0d0 * ss_PotTempBasicZ )

  
  !=== 速度シアーによる乱流エネルギー生成
  ss_ShearKm = &
       &  ( Cm ** 2.0d0 ) * ( DelXZ ** 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 * ( DelXZ ** 2.0d0 ) )
  
  
  !=== 項をまとめる
  ss_TendKm = - ss_AdvKm  - ss_BuoyKm  +  ss_ShearKm &
       &      + ss_DiffKm - ss_DispKm 

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


  !=== 境界条件
  call boundary( ss_BC, ss_Km_out )
  

  !=== 負の領域を零に固定
  where ( ss_Km_out < 0.0d0 )
     ss_Km_out = 0.0d0
  end where

  
  !=== 渦拡散定数を決める
  ss_Kh_out = 3.0d0 * ss_Km_out


  call EndSub("EddyViscosity")
  
end subroutine EddyViscosity
