!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine DensityCloud
!
!   * Developer: Kitamori Taichi
!   * Version: $Id: densitycloud_precipitation.f90,v 1.1 2009-03-05 05:39:42 yamasita Exp $ 
!   * Tag Name: $Name: arare4-20100305 $
!   * Change History: 
!
!== Overview 
!
!雲密度の時間発展を解く
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!   * フラックス形式の方程式を解いている
!   * 乱流拡散項, 数値粘性項は含まれていない 
!   * 雲粒の重力落下項を考慮する
!  
!== Future Plans
!
                                                                 !=end
subroutine DensityCloud(      & 
  & xz_DensCloud_in,                 & !(in) 雲の密度
  & DelTime,                         & !(in) leapfrog スキームの時間間隔
  & pz_VelXNl,                      & !(in) 水平流速
  & xr_VelZNl,                      & !(in) 鉛直流速
  & xz_PotTempNl,                   & !(in) 温位擾乱
  & xz_ExnerNl,                     & !(in) エクスナー擾乱
  & xz_DensCloudNl,                 & !(in) 雲の密度
  & xz_RadiCloud,                   & !(in) 雲粒半径
  & xz_MasscondNl,                  & !(inout) 凝結量
  & xz_DensCloud_out)                  !(out) 雲の密度

                                                                 !=begin  
  !== Dependancy
  use dc_trace, only: BeginSub, EndSub, DbgMessage
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax
!  use bcset,    only: xz_BC
  use average,  only: pz_avr_xz, xr_avr_xz
  use differentiate_center4, only: xz_dx_pz, xz_dz_xr
  use basicset, only: Grav, PressSfc, GasRDry, CpDry, &
               &      xz_ExnerBasicZ, xz_PotTempBasicZ
  use cloudset, only: DensIce

                                                                 !=end
  !== 暗黙の型宣言禁止
  implicit none
                                                                 !=begin
  !== Input
  real(8), intent(in)  :: DelTime       ! leapfrog スキームの時間間隔 [s]
  real(8), intent(in)  :: xz_DensCloud_in(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲の密度 [kg/m^3]
  real(8), intent(in)  :: pz_VelXNl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 水平速度 [m/s]
  real(8), intent(in)  :: xr_VelZNl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 鉛直速度 [m/s]
  real(8), intent(in)  :: xz_DensCloudNl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲の密度 [kg/m^3]
  real(8), intent(in)  :: xz_RadiCloud(DimXMin:DimXMax, DimZMin:DimZMax)  
                                        ! 雲粒の半径 [m]
  real(8), intent(in)  :: xz_ExnerNl(DimXMin:DimXMax, DimZMin:DimZMax)  
                                        ! エクスナー関数の擾乱成分 [1]
  real(8), intent(in)  :: xz_PotTempNl(DimXMin:DimXMax, DimZMin:DimZMax)  
                                        ! 温位の擾乱成分 [K]

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

  !== In/Out
  real(8), intent(inout):: xz_MasscondNl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 凝結量 [kg/m^3 s]
                                                                 !=end  
  !== Work
  real(8)  :: xz_FluxDensCloud(DimXMin:DimXMax, DimZMin:DimZMax) 
                                        ! フラックスによる雲密度の時間変化率 [kg/m^3 s]
  real(8)  :: xz_SrcDensCloud(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 生成消滅による雲密度の時間変化率 [kg/m^3 s]
  real(8)  :: xz_PrecipDensCloud(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 落下による雲密度の時間変化率 [kg/m^3 s]
  real(8)  :: xz_TermVel(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲粒の終端速度
  real(8)  :: xz_VisCoffCO2(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! CO2 の粘性係数
  real(8)  :: xz_MeanFreePath(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! CO2 の平均自由行程
  real, parameter :: Pi = 3.1415926535897932385d0
                                        ! 円周率
  real, parameter :: BoltzmannConst = 1.38065d-23
                                        ! ボルツマン定数
  real, parameter :: CO2Diam = 3.30d-10
                                        ! CO2 分子の kinetic diameter
  real, parameter :: SutherlandConst = 2.40d2
                                        ! CO2 分子のサザランド定数
  real, parameter :: VisCoffConst = 1.47d-5
                                        ! 1 気圧, 20 度でのCO2 分子の粘性係数


  call BeginSub("DensityCloud", &
       &        fmt="%c",              &
       &        c1="Time integration of density of cloud.")
  
!  write(*,*) xz_DensCloudNl(1,:)

  !=== 粘性係数の計算
  xz_VisCoffCO2 = &
    & VisCoffConst * (293.15d0 + SutherlandConst) &
    & / ( (xz_PotTempBasicZ + xz_PotTempNl)*(xz_ExnerBasicZ + xz_ExnerNl) &
    &     + SutherlandConst ) &
    & * ( (xz_PotTempBasicZ + xz_PotTempNl)*(xz_ExnerBasicZ + xz_ExnerNl) &
    &     / 293.15d0 )**1.5d0
  !=== 平均自由行程の計算
  xz_MeanFreePath = &
    & BoltzMannConst * &
    & (xz_PotTempBasicZ + xz_PotTempNl)*(xz_ExnerBasicZ + xz_ExnerNl) &
    & / ( 2.0**0.5d0 * Pi * (CO2Diam)**2.0d0 * PressSfc &
    &    * (xz_ExnerBasicZ + xz_ExnerNl)**(CpDry/GasRDry) )

  !=== 終端速度の計算
  xz_TermVel = &
    & 2.0d0 * Grav * DensIce * xz_DensCloudNl  &
    & * (xz_RadiCloud + 1.255d0 * xz_MeanFreePath ) &
    & / (9.0d0 * xz_VisCoffCO2)

  !=== 落下項の計算. 4 次精度中心差分を用いて計算
  xz_PrecipDensCloud  = &
    & xz_dz_xr(xr_avr_xz(xz_TermVel) * xr_avr_xz(xz_DensCloudNl))


  !=== フラックス項の計算. 4 次精度中心差分を用いて計算
  xz_FluxDensCloud  = &
    & - xz_dx_pz(pz_VelXNl * pz_avr_xz(xz_DensCloudNl)) &
    & - xz_dz_xr(xr_VelZNl * xr_avr_xz(xz_DensCloudNl))


  call DbgMessage( fmt="*** Debug Message [DensityCloud] *** xz_SrcDensCloud ")

  !=== 生成項の計算
  xz_SrcDensCloud = xz_MasscondNl

  call DbgMessage( fmt="*** Debug Message [DensityCloud] *** xz_DensCloud_out ")

  xz_DensCloud_out = xz_DensCloud_in  & 
    & + DelTime * (xz_PrecipDensCloud + xz_FluxDensCloud + xz_SrcDensCloud)

  call DbgMessage( fmt="*** Debug Message [DensityCloud] *** xz_DensCloud_out(adjustment) ")

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

!  !=== 境界条件
!  call boundary(xz_BC, xz_DensCloud_out)
!
  call EndSub("DensityCloud")
  
end subroutine DensityCloud
