!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine VelocityZ
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: velocityz.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!オイラー法を用いた Z 方向の速度の計算. 
!deepconv/arare では時間積分として He-VI 法を利用しているので, 
!Z 方向の速度は陰解法で解く. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * エクスナー関数の空間方向の離散化において, 2 次精度の離散化を陽に利用しているため, 本サブプログラムにおいても differentiate_center4 モジュールを指定することはできないので注意.
!
!== Future Plans
!
                                                                 !=end
subroutine VelocityZ(  &
     & sf_VelZ_in, &
     & DelTime,    &
     & sf_AdvDiffZ, fs_VelX, sf_VelZ, ss_Exner_n, ss_Exner_a, &
     & sf_VelZ_out )
                                                                 !=begin
  !== Dependency
  use dc_trace, only: BeginSub, EndSub  
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax
  use basicset, only: ss_CpBasicZ, ss_PotTempBasicZ
  use bcset,    only: sf_BC
  use arareset, only: alpha, beta
  use average,  only: sf_avr_ss
  use differentiate_center2,  only: ss_dx_fs, ss_dz_sf, sf_dz_ss

                                                                 !=end  
  !== 暗黙の型宣言禁止
  implicit none
                                                                 !=begin  
  !== Input
  real(8), intent(in)   :: DelTime
  real(8), intent(in)   :: sf_VelZ_in(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)   :: ss_Exner_n(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)   :: ss_Exner_a(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)   :: fs_VelX(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)   :: sf_VelZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)   :: sf_AdvDiffZ(DimXMin:DimXMax, DimZMin:DimZMax)

  !== Output
  real(8), intent(out)  :: sf_VelZ_out(DimXMin:DimXMax, DimZMin:DimZMax)
                                                                 !=end  
  !== Work
  real(8)               :: ss_DivVel(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)               :: sf_PressGradVelZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)               :: sf_TendVelZ(DimXMin:DimXMax, DimZMin:DimZMax)


  call BeginSub("VelocityZ", &
&           fmt="%c",        &
&           c1="Calculate VelZ_backward (sf_VelZ).")

  !=== 速度の収束
  ss_DivVel =  ss_dx_fs( fs_VelX ) + ss_dz_sf( sf_VelZ )
  

  !=== 速度 w の圧力勾配
  sf_PressGradVelZ =   &
       & - sf_avr_ss( ss_CpBasicZ * ss_PotTempBasicZ )      &
       &   * (                                              &
       &       beta * sf_dz_ss( ss_Exner_a )               &
       &       + (1.0d0 - beta) * sf_dz_ss( ss_Exner_n )   &
       &       - sf_dz_ss( alpha * ss_DivVel )           &
       &     )                                              


  !=== 微分係数の評価
  sf_TendVelZ = sf_PressGradVelZ + sf_AdvDiffZ

  
  !=== 数値積分
  sf_VelZ_out = sf_VelZ_in + DelTime * sf_TendVelZ


  !=== 境界条件
  call boundary(sf_BC, sf_VelZ_out)
  

  call EndSub("VelocityZ")
  
end subroutine VelocityZ
