!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine TurbulenceVelZ
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: turbulencevelz.f90,v 1.3 2005/02/02 09:53:48 odakker Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!Z 方向の運動方程式の粘性拡散項を評価する. 
!音波に関連しない項として扱うため長い時間ステップで評価. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
                                                                 !=end

subroutine TurbulenceVelZ(fs_VelX_n, sf_VelZ_n, ss_Km_n, sf_TurbVelZ_n)
                                                                 !=begin
  !== Dependency
  use dc_trace, only: BeginSub, EndSub
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, DelXZ
  use differentiate_center2,  only: sf_dx_ff, ff_dz_fs, ff_dx_sf, &
       &                    sf_dz_ss, ss_dz_sf
  use average,        only: ff_avr_ss
  use arareprm, only: Cm
                                                                 !=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_Km_n(DimXMin:DimXMax, DimZMin:DimZMax)   !渦粘性係数

  !== Output
  real(8), intent(out) :: sf_TurbVelZ_n(DimXMin:DimXMax, DimZMin:DimZMax) !粘性拡散項
                                                                 !=end  
  call BeginSub("TurbulenceVelZ", &
&              fmt="%c",       &
&              c1="Calculate terbulent diffusion term for VelZ (sf_TurbVelZ).")

  !--- 拡散項
  sf_TurbVelZ_n = &
       &   sf_dx_ff( ff_avr_ss( ss_Km_n ) * ff_dz_fs( fs_VelX_n ) )   &
       & + sf_dx_ff( ff_avr_ss( ss_Km_n ) * ff_dx_sf( sf_VelZ_n ) )   &
       & + 2.0d0 * sf_dz_ss( ss_Km_n * ss_dz_sf( sf_VelZ_n ) )  &
       & - 2.0d0 * sf_dz_ss( ( ss_Km_n ** 2.0d0 ) ) &
       &         / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( DelXZ ** 2.0d0 ) )

  write(*,*) "Dz ", maxval( abs( sf_TurbVelZ_n ) )

  call EndSub("TurbulenceVelZ")
    
end subroutine TurbulenceVelZ
