!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine DensityCloud
!
!   * Developer: Kitamori Taichi
!   * Version: $Id: densitycloud.f90,v 1.1.2.2 2005/12/21 05:47:10 kitamo Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!雲密度の時間発展を解く
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!   * フラックス形式の方程式を解いている
!   * 乱流拡散項は含まれていない 
!   * 北守修論によると, 雲粒落下の速さは十分小さいので, 雲粒の落下の効果は無視(山下, 080411)
!  
!== Future Plans
!
                                                                 !=end
subroutine DensityCloud(      & 
  & ss_DensCloud_in,                 & !(in) 雲の密度
  & DelTime,                         & !(in) leapfrog スキームの時間間隔
  & fs_VelX_nl,                      & !(in) 水平流速
  & sf_VelZ_nl,                      & !(in) 鉛直流速
  & ss_DensCloud_nl,                 & !(in) 雲の密度
  & ss_MassCond_nl,                  & !(inout) 凝結量
  & ss_DensCloud_out)                  !(out) 雲の密度

                                                                 !=begin  
  !== Dependancy
  use dc_trace, only: BeginSub, EndSub
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax
  use bcset,    only: ss_BC
  use average,  only: fs_avr_ss, sf_avr_ss
  use differentiate_center4, only: ss_dx_fs, ss_dz_sf
                                                                 !=end
  !== 暗黙の型宣言禁止
  implicit none
                                                                 !=begin
  !== Input
  real(8), intent(in)  :: DelTime       ! leapfrog スキームの時間間隔 [s]
  real(8), intent(in)  :: ss_DensCloud_in(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲の密度 [kg/m^3]
  real(8), intent(in)  :: fs_VelX_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 水平速度 [m/s]
  real(8), intent(in)  :: sf_VelZ_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 鉛直速度 [m/s]
  real(8), intent(in)  :: ss_DensCloud_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲の密度 [kg/m^3]

  !== Output
  real(8), intent(out)  :: ss_DensCloud_out(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲の密度 [kg/m^3]

  !== In/Out
  real(8), intent(inout):: ss_MassCond_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 凝結量 [kg/m^3 s]
                                                                 !=end  
  !== Work
  real(8)  :: ss_FluxDensCloud(DimXMin:DimXMax, DimZMin:DimZMax) 
                                        ! フラックスによる雲密度の時間変化率 [kg/m^3 s]
  real(8)  :: ss_SrcDensCloud(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 生成消滅による雲密度の時間変化率 [kg/m^3 s]

  call BeginSub("DensityCloud", &
       &        fmt="%c",              &
       &        c1="Time integration of density of cloud.")
  
  !=== フラックス項の計算. 4 次精度中心差分を用いて計算
  ss_FluxDensCloud  = &
    & - ss_dx_fs(fs_VelX_nl * fs_avr_ss(ss_DensCloud_nl)) &
    & - ss_dz_sf(sf_VelZ_nl * sf_avr_ss(ss_DensCloud_nl))
!  ss_FluxDensCloud = 0.0d0

  !=== 生成項の計算
  ss_SrcDensCloud = ss_MassCond_nl

  ss_DensCloud_out = ss_DensCloud_in  & 
    & + DelTime * (ss_FluxDensCloud + ss_SrcDensCloud)

  ! 雲の量が負になる場合は 0 にする. 
  ! 差し引き分だけ蒸発量を減らす
  ss_DensCloud_out = max(ss_DensCloud_out, 0.0d0)
  ss_MassCond_nl = max(ss_MassCond_nl,   &
    &                  - (ss_DensCloud_nl / DelTime  + ss_FluxDensCloud))

  !=== 境界条件
  call boundary(ss_BC, ss_DensCloud_out)

  call EndSub("DensityCloud")
  
end subroutine DensityCloud