!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine EddyViscosity
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: eddyviscosity.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!1.5 $B<!$NMpN.%/%m!<%8%c!<$rMxMQ$7(B, 
!$B<!$N;~9o$N12G4@-78?t$H123H;678?t$r5a$a$k%5%V%k!<%A%s(B.
!
!== Error Handling
!
!== Known Bugs
!
!  * 2005/09/05 $B8=:_(B, $B6u5$Cf$N4%Ag@.J,$H6E=L@-J,$NJ,;RNL$rF1$8$K$7$F$"$k(B
!  * $BIbNO9`$NJ,;RNL8z2L$r9MN8$7$F$$$J$$(B
!
!== Note
!
!  * $BDj<02=$O(B CReSS $B$N%^%K%e%"%k$r;2>H$7$?(B
!  * $B0\N.9`$N7W;;$K$O(B ***_adv $B$H$J$kJQ?t$rMxMQ$9$k(B. 
!  * $B@8@.9`(B, $B3H;69`(B, $B>C;69`$N7W;;$K$O(B ***_dif $B$H$J$kJQ?t$rMxMQ$9$k(B. 
!
!== Future Plans
!
                                                                 !=end
subroutine EddyViscosity(                         &
     & ss_Km_in,                                  &
     & DelTime,                                   & 
     & fs_VelX_adv, sf_VelZ_adv, ss_Km_adv,       & !$B0\N.(B
     & fs_VelX_dif, sf_VelZ_dif, ss_PotTemp_dif,  & !$B@8@.(B, $B3H;6(B, $B>C;6(B
     & 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, DelX, DelZ
  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
  !== $B0EL[$N7?@k8@6X;_(B
  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_Km_tmp(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)
  real(8), parameter :: alpha = 1.0d-6 ! $B12G4@-78?t$r7h$a$k:]$KMQ$$$k78?t(B. $BCM$O(B CReSS $B;2>H(B. 


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


  !==== $BJQ?t=i4|2=(B
  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_Km_tmp  = 0.0d0
  ss_Temp_dif = 0.0d0
  ss_AA       = 0.0d0
  ss_MixRatSum_dif = 0.0d0


  !=== $B%9%+%i!<NL$N0\N.(B. $BCf?4:9J,(B. 
  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 ) ) 


  !=== $BIbNO$K$h$kMpN.%(%M%k%.!<@8@.(B
  !----- $BK\Ev$J$i(B, $B290L$rAjEv290L$KFI$_BX$($kI,MW$,$"$k$,(B, $B8=:_$O(B
  !----- $B4%Ag@.J,$H6E=L@.J,$NJ,;RNL$r0lDj$H$7$F$$$k$N$G(B, $B$=$NI,MW$O$J$$(B
  ss_Temp_dif = &
       &   ( ss_PotTempBasicZ + ss_PotTemp_dif ) &
       & * ( ss_ExnerBasicZ   + ss_Exner_dif   )
  ss_AA =  ( 1.0d0 + 1.61d0 * epsilon * LatentHeat * ss_MixRtV_dif &
       &      / ( GasRDry * ss_Temp_dif ) )                   &
       & / ( 1.0d0 + epsilon * ( LatentHeat ** 2.0d0 ) * ss_MixRtV_dif &
       &      / ( CpDry * GasRDry * ( ss_Temp_dif ** 2.0d0 ) ) ) &
       & / ( ss_PotTempBasicZ + ss_PotTemp_dif )
  ss_MixRatSum_dif = ss_MixRtC_dif + ss_MixRtR_dif
  
  
  where ( ss_MixRtC_dif + ss_MixRtR_dif > Min ) 
     ss_BuoyKm = &
          &   1.5d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) &
          & * ( &
          &       ss_AA * ss_avr_sf( sf_dz_ss( ss_PotTemp_dif ) ) &
          &     - ss_avr_sf( sf_dz_ss( ss_MixRatSum_dif ) )  &
         &   )
     
  elsewhere ( ss_MixRtC_dif + ss_MixRtR_dif < Min ) 
     ss_BuoyKm = &
          &   1.5d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 )     &
         & * (                                                        &
          &       ss_avr_sf( sf_dz_ss( ss_PotTemp_dif ) )              &
          &         / ( ss_PotTempBasicZ + ss_PotTemp_dif )            &
          &     + 0.61d0 * ss_avr_sf( sf_dz_ss( ss_MixRtV_dif ) ) &
          &   )
     
  end where

!$B0J2<(B, $BJ,;RNL8z2L$r=|$$$?>l9g$NIbNO$K$h$kMpN.%(%M%k%.!<@8@.(B
!     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 )
!
  
  !=== $BB.EY%7%"!<$K$h$kMpN.%(%M%k%.!<@8@.(B
  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

  
  !=== $BMpN.%(%M%k%.!<$N3H;6(B
  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
  

  !=== $BMpN.%(%M%k%.!<$N>C;6(B
  ss_DispKm = ( ss_Km_dif ** 2.0d0 ) / ( 2.0d0 * ( MixLen ** 2.0d0 ) )
  
  
  !=== $B9`$r$^$H$a$k(B
  ss_TendKm = - ss_AdvKm - ss_BuoyKm + ss_ShearKm + ss_DiffKm - ss_DispKm 

  
  !=== $B;~4V@QJ,(B
  ss_Km_tmp = ss_Km_in + DelTime * ss_TendKm


  !=== $B6-3&>r7o(B
  call boundary( ss_BC, ss_Km_out )
  

  !=== $BIi$NNN0h$KBP$9$k07$$(B (2005-10-05, CReSS $B$r;29M$K(B)
  if ( DelX >= DelZ ) then 
     ss_Km_out = max( ss_Km_tmp, alpha * ( DelZ **2.0d0 ) )
  else
     ss_Km_out = max( ss_Km_tmp, alpha * ( DelX **2.0d0 ) )
  end if
  
  
  !=== $B123H;6Dj?t$r7h$a$k(B
  ss_Kh_out = 3.0d0 * ss_Km_out

  
  call EndSub("EddyViscosity")
  
end subroutine EddyViscosity
