!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine AdvectDiffuseZ
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: advectdiffusez.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!Z $BJ}8~$N1?F0J}Dx<0Cf$K8=$l$k0J2<$N9`$r(B
!$BD9$$;~4V%9%F%C%W(B DelTimeLong $B$K$*$$$FI>2A$9$k(B. 
!
!   * $B0\N.9`(B
!   * $BIbNO9`(B
!   * $BMpN.3H;69`(B
!   * $B?tCMG4@-9`(B
!
!== Error Handling
!
!== Known Bugs
!
!  * 2005-10-05 $B8=:_(B, $B6u5$Cf$N4%Ag@.J,$H6E=L@.J,$NJ,;RNL$OF1$8$H$7$F$$$k(B
!  * $BIbNO9`$NCf$NJ,;RNL8z2L$O9MN8$7$J$$(B. 
!
!== Note
!
!  * $B0\N.9`$N7W;;$K$O(B ***_adv $B$H$J$kJQ?t$rMxMQ$9$k(B. 
!  * $BMpN.3H;69`(B, $B?tCMG4@-9`$N7W;;$K$O(B ***_dif $B$H$J$kJQ?t$rMxMQ$9$k(B. 
!
!== Future Plans
!
                                                                 !=end
subroutine AdvectDiffuseZ(                       &
     & fs_VelX_adv, sf_VelZ_adv, ss_PotTemp_adv, & !$B0\N.(B, $BIbNO(B
     & fs_VelX_dif, sf_VelZ_dif, ss_Km_dif,      & !$BMpN.3H;6(B, $B?tCMG4@-(B
     & ss_MixRtV_adv, ss_MixRtC_adv, ss_MixRtR_adv, &
     & sf_Buoy, &
     & sf_AdvDiffZ )
                                                                 !=begin
  !== Dependency
  use dc_trace, only: BeginSub, EndSub
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, MixLen
  use physset,  only: Grav, epsilon
  use basicset, only: ss_PotTempBasicZ, ss_MixRtVBasicZ
  use arareset, only: Cm, NuH, NuV
  use average,  only: sf_avr_ff, sf_avr_ss, ss_avr_sf, ff_avr_fs, ff_avr_ss
  use differentiate_center4, only: ff_dx_sf, sf_dx_ff, &
       &                           ss_dz_sf, ff_dz_fs, sf_dz_ss 
                                                                 !=end  
  !== $B0EL[$N7?@k8@6X;_(B 
  implicit none
                                                                 !=begin
  !== Input
  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)  :: fs_VelX_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: sf_VelZ_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_Km_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_PotTemp_adv(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_MixRtV_adv(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_MixRtC_adv(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_MixRtR_adv(DimXMin:DimXMax, DimZMin:DimZMax)

  !== Output
  real(8), intent(out) :: sf_AdvDiffZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(out)  :: sf_Buoy(DimXMin:DimXMax, DimZMin:DimZMax)
                                                                 !=end    
  !=== Work
  real(8)  :: sf_AdvVelZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: sf_TurbVelZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: sf_NumDiffVelZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: sf_DelMixRtV_adv(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: sf_DelMixRatAll_adv(DimXMin:DimXMax, DimZMin:DimZMax)

  call BeginSub("AdvectDiffuseZ",            &
       &        fmt="%c",                    &
       &        c1="calculation of extra force term in z-direction.")

  !=== $B=i4|2=(B
  sf_AdvVelZ            = 0.0d0
  sf_TurbVelZ           = 0.0d0
  sf_NumDiffVelZ        = 0.0d0
  sf_DelMixRtV_adv = 0.0d0
  sf_DelMixRatAll_adv   = 0.0d0

  !=== $B0\N.9`(B
  sf_AdvVelZ = &
       &   sf_avr_ff( ff_avr_fs( fs_VelX_adv ) * ff_dx_sf( sf_VelZ_adv ) ) &
       & + sf_avr_ss( ss_avr_sf( sf_VelZ_adv ) * ss_dz_sf( sf_VelZ_adv ) ) 

  
  !=== $BIbNO9`(B 
  sf_DelMixRtV_adv = &
       & sf_avr_ss( ss_MixRtV_adv - ss_MixRtVBasicZ )  
  sf_DelMixRatAll_adv = &
       & sf_avr_ss(   &
       &              ss_MixRtV_adv - ss_MixRtVBasicZ &
       &            + ss_MixRtC_adv + ss_MixRtR_adv     )  
  sf_Buoy = Grav &
       & * (                                                          &
       &      sf_avr_ss( ss_PotTemp_adv  / ss_PotTempBasicZ )         &
       &    + sf_DelMixRtV_adv                                  &
       &       / ( epsilon +  sf_avr_ss( ss_MixRtVBasicZ ) )   &
       &    - sf_DelMixRatAll_adv                                    &
       &       / ( 1.0d0 + sf_avr_ss( ss_MixRtVBasicZ ) )      &
       &    )
!  sf_Buoy = Grav * ( &
!       &   sf_avr_ss( ss_PotTemp_adv  / ss_PotTempBasicZ ) )
  
  
  !=== $BMpN.3H;69`(B ($B;~9o$O(B t - $B&$(Bt $B$GI>2A(B)
  sf_TurbVelZ = &
       &   sf_dx_ff( ff_avr_ss( ss_Km_dif ) * ff_dz_fs( fs_VelX_dif ) )   &
       & + sf_dx_ff( ff_avr_ss( ss_Km_dif ) * ff_dx_sf( sf_VelZ_dif ) )   &
       & + 2.0d0 * sf_dz_ss( ss_Km_dif * ss_dz_sf( sf_VelZ_dif ) )  &
       & - 2.0d0 * sf_dz_ss( ( ss_Km_dif ** 2.0d0 ) ) &
       &         / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )


  !=== $B?tCMG4@-9`(B ($B;~9o$O(B t - $B&$(Bt $B$GI>2A(B)
  sf_NumDiffVelZ = &
       &   NuH * ( sf_dx_ff( ff_dx_sf( sf_VelZ_dif) ) ) &
       & + NuV * ( sf_dz_ss( ss_dz_sf( sf_VelZ_dif) ) ) 


  !=== $B30NO9`$H$7$F=PNO(B
  sf_AdvDiffZ = - sf_AdvVelZ + sf_TurbVelZ + sf_NumDiffVelZ + sf_Buoy 

  call EndSub("AdvectDiffuseZ")
  
end subroutine AdvectDiffuseZ
