!---------------------------------------------------------------------
!     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/01/28 12:41:21 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!乱流エネルギーを定常とした場合の渦拡散係数, 渦粘性係数を求める. 
!この乱流パラメタリゼーションは, Mellor and Yamada (1974) の
!Level 1 Closure に対応するが, Level1 Closure に存在する 
!\bar{\theta^2} は無視されている. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!乱流エネルギーが負になったら, 零を入れる.
!
!== Future Plans
!
                                                                 !=end

subroutine EddyViscosity(fs_VelX_n, sf_VelZ_n, ss_PotTemp_n, &
     &                   ss_Km_n, ss_Kh_n)
                                                                 !=begin
  !== Dependency
  use timeset,  only: DelTimeLong
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, DelX
  use average,  only: ss_avr_ff, ss_avr_sf
  use arareprm, only: Cm, Grav
  use differentiate, only: ss_dx_fs, sf_dz_ss, ff_dx_sf, ff_dz_fs, ss_dz_sf
  use basicset, only: ss_PotTempBasicZ
  use bcset,    only: ss_BC
                                                                 !=end  
  !--- 暗黙の型宣言禁止
  implicit none
                                                                 !=begin  
  !== Input
  real(8), intent(in)  :: &
       & fs_VelX_n(DimXMin:DimXMax, DimZMin:DimZMax) !速度 u
  real(8), intent(in)  :: &
       & sf_VelZ_n(DimXMin:DimXMax, DimZMin:DimZMax) !速度 w
  real(8), intent(in)  :: &
       & ss_PotTemp_n(DimXMin:DimXMax, DimZMin:DimZMax) !温位

  !== Output
  real(8), intent(out) :: &
       & ss_Km_n(DimXMin:DimXMax, DimZMin:DimZMax) !渦粘性係数
  real(8), intent(out) :: &
       & ss_Kh_n(DimXMin:DimXMax, DimZMin:DimZMax) !渦拡散係数
                                                                 !=end  
  !--- 内部変数
  real(8)  :: ss_Km2_n(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: Max
  

  !--- 渦粘性係数
  ss_Km2_n = &
       & ( Cm ** 2.0d0 ) * ( DelX ** 4.0d0 ) &
       &   * (                                & 
       &        2.0d0 * ( ss_dx_fs( fs_VelX_n ) ** 2.0d0 )  &
       &      + 2.0d0 * ( ss_dz_sf( sf_VelZ_n ) ** 2.0d0 )  &
       &      + ss_avr_ff(                                  &
       &          ( ff_dx_sf( sf_VelZ_n ) + ff_dz_fs( fs_VelX_n ) )** 2.0d0 &
       &        ) &
       &      - 3.0d0 * Grav * ss_avr_sf( &
       &          ( sf_dz_ss( ss_PotTemp_n ) / ss_PotTempBasicZ )  &
       &        ) &
       &      )

  !--- 値の調整. deepconv 参照
  where ( ss_Km2_n < 0.0d0 ) 
     ss_Km2_n = 0.0d0
  end where

  Max = ( DelX ** 2.0d0 ) / ( 2.0d0 * DelTimeLong )
  where ( ss_Km2_n > (Max ** 2.0d0) ) 
     ss_Km2_n = Max ** 2.0d0
  end where

  !--- 出力
  ss_Km_n = sqrt( ss_Km2_n )
  ss_Kh_n = 3.0d0 * ss_Km_n 

!  ss_Km_n = 0.0d0
!  ss_Kh_n = 0.0d0

  !--- 境界条件
  call boundary( ss_BC, ss_Km_n )
  call boundary( ss_BC, ss_Kh_n )
  
end subroutine EddyViscosity
