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


subroutine TurbulenceVelX(fs_VelX_n, sf_VelZ_n, ss_Km_n, fs_TurbVelX_n)
                                                                 !=begin  
  !== Dependency 
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, DelXZ
  use differentiate,  only: fs_dx_ss, ss_dx_fs, fs_dz_ff, ff_dx_sf, &
       &                    fs_dz_ff, ff_dz_fs, fs_dx_ss
  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) :: fs_TurbVelX_n(DimXMin:DimXMax, DimZMin:DimZMax) !粘性拡散項
                                                                 !=end


  !--- 速度 u に対する拡散項 
  fs_TurbVelX_n = &
       &   2.0d0 * fs_dx_ss( ss_Km_n * ss_dx_fs( fs_VelX_n ) ) &
       & + fs_dz_ff( ff_avr_ss( ss_Km_n ) * ff_dx_sf( sf_VelZ_n ) )  &
       & + fs_dz_ff( ff_avr_ss( ss_Km_n ) * ff_dz_fs( fs_VelX_n ) )  &
       & - 2.0d0 * fs_dx_ss( ( ss_Km_n ** 2.0d0 ) ) &
       &   / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( DelXZ ** 2.0d0 ) )   

  write(*,*) "Dx ", maxval( abs( fs_TurbVelX_n ) )

  
end subroutine TurbulenceVelX
