!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine DensityCloud
!
!   * Developer: Kitamori Taichi
!   * Version: $Id: densitycloud_falltest3.f90,v 1.2 2011-12-15 03:42:15 yamasita Exp $ 
!   * Tag Name: $Name: arare4-20111215 $
!   * Change History: 
!
!== Overview 
!
! 雲密度の時間発展を解く
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!   * フラックス形式の方程式を解いている
!   * 雲密度落下計算用テストプログラム.
!     * 乱流拡散項, 数値粘性項, 移流項は陽にゼロを代入.
!     * 雲粒半径は雲密度のみに依存するものとする. 
!     * Stokes 則の係数部分は雲粒半径のみに依存するものとする. 
!     * Knudsen 数において, 平均自由行程は定数であるとする. 
!  
!== Future Plans
!
                                                                 !=end
subroutine DensityCloud(      &
  & xz_DensCloudNs,                 & !(in) 雲の密度
  & DelTimeShort,                   & !(in) 時間ステップ
  & xz_DensCloudWorkAs)               !(out) 雲の密度(負の密度の調整前の値)

                                                                 !=begin  
  !== Dependancy
  use dc_trace, only: BeginSub, EndSub, DbgMessage
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax
  use average,  only: xr_avr_xz
  use differentiate_center4, only: xz_dz_xr
  use basicset, only: Grav, PressSfc, GasRDry, CpDry, &
               &      xz_ExnerBasicZ, xz_PotTempBasicZ, &
               &      xz_DensBasicZ
  use cloudset, only: DensIce, &        ! 固相物質の密度
    &                 NumAerosol, &     ! 質量あたりの凝結核数密度
    &                 RadiAerosol       ! 凝結核半径

  use StoreDensCloud, only: StoreDensCloudFall
                                                                 !=end
  !== 暗黙の型宣言禁止
  implicit none
                                                                 !=begin
  !== Input
  real(8), intent(in)  :: DelTimeShort  ! タイムステップ [s]
  real(8), intent(in)  :: xz_DensCloudNs(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲の密度 [kg/m^3]
  !== Output
  real(8), intent(out)  :: xz_DensCloudWorkAs(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲の密度(負の密度調整前) [kg/m^3]

                                                                 !=end  
  !== Work

  real(8)  :: xz_TendDensCloudNl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 長い時間ステップで評価する項 [kg/m^3 s]
  real(8)  :: xz_RadiCloudNs(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲粒の半径 [m]
  real(8)  :: xz_FallDensCloud(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲粒落下による雲密度の時間変化率 [kg/m^3 s]
  real(8)  :: xz_VelZTerm(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲粒の終端速度
  real(8)  :: xz_VisCoffCO2(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! CO2 の粘性係数
  real(8)  :: xz_MeanFreePath(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! CO2 の平均自由行程
  real(8)  :: xz_KnudsenNum(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲粒に対する Knudsen 数
  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 :: VisCoffRef = 1.47d-5
                                        ! 粘性係数の参照値
  real, parameter :: TempRef = 2.93d2
                                        ! 粘性係数の参照温度

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

  !=== 粘性係数の計算
  xz_VisCoffCO2 = &
    & VisCoffRef * (TempRef + SutherlandConst) &
    & / ( xz_PotTempBasicZ * xz_ExnerBasicZ  &
    &     + SutherlandConst ) &
    & * ( xz_PotTempBasicZ * xz_ExnerBasicZ  &
    &     / TempRef )**1.5d0

  !=== 平均自由行程の計算
!  xz_MeanFreePath = &
!    & BoltzMannConst * &
!    & xz_PotTempBasicZ * xz_ExnerBasicZ  &
!    & / ( 2.0**0.5d0 * Pi * (CO2Diam)**2.0d0 * PressSfc &
!    &    * (xz_ExnerBasicZ)**(CpDry/GasRDry) )
  xz_MeanFreePath = 1.0d-5 ! テスト用

  !=== 雲粒半径の計算
!  xz_RadiCloudNs =   &
!    & (  &
!    &      3.0d0 * xz_DensCloudNs &
!    &      / (4.0d0 * Pi * DensIce * NumAerosol * xz_DensBasicZ) &
!    &      + RadiAerosol**3.0d0  &
!    &   )**(1.0d0 / 3.0d0)
  xz_RadiCloudNs = &
    & ( 3.0d-11 * xz_DensCloudNs + RadiAerosol**3.0d0 )**(1.0d0 / 3.0d0)

  !=== Knudsen 数の計算
  xz_KnudsenNum = &
    & xz_MeanFreePath / xz_RadiCloudNs

  !=== 終端速度の計算
!  xz_VelZTerm = &
!    & ( 1.0d0 + 4.0d0 * xz_KnudsenNum / 3.0d0 ) &
!    & * 2.0d0 * xz_RadiCloudNs * xz_RadiCloudNs * Grav * DensIce  &
!    & / (9.0d0 * xz_VisCoffCO2)
  xz_VelZTerm = &
    & ( 1.0d0 + 4.0d0 * xz_KnudsenNum / 3.0d0 ) &
    & * 2.0d8 * xz_RadiCloudNs**2.0d0

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

  call StoreDensCloudFall( xz_FallDensCloud )

  xz_DensCloudWorkAs = xz_DensCloudNs  &
    & + DelTimeShort * xz_FallDensCloud

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