!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine MassCondense
!
!   * Developer: KITAMORI Taichi 
!   * Version: $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!凝結/蒸発量を計算する.
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!拡散によって雲粒が成長する場合における成長方程式を解く. 
!蒸発量が存在する雲の量よりも多い場合, 蒸発量を存在する雲の量にする. 
!
!== Future Plans
!
                                                                 !=end
subroutine MassCondense(  &
  & xz_LatHeatPerMassNl,  &  !(in) 単位質量あたりの潜熱
  & xz_SatRatioNl, &         !(in) 過飽和度
  & xz_PotTempNl,     &      !(in) 温位
  & xz_ExnerNl,       &      !(in) エクスナー関数
  & xz_DensCloudNl,     &    !(in) 雲密度
  & xz_MassCondNl)           !(out) 凝結量(山下, 20080410, 蒸発量は入っているのか要確認)
                                                                 !=begin
  !==Dependency
  !(山下, 20080410, モジュールとモジュールの中の変数がちゃんと使えるようになっているのか要確認)
  use dc_trace, only: BeginSub, EndSub, DbgMessage
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax
!  use physset,  only: GasR,   &         ! 気体定数  
!    &                 Pi,     &         ! 円周率
!    &                 PressSfc          ! 地表面圧力
!  use physset,  only: Pi,     &         ! 円周率
!    &                 PressSfc          ! 地表面圧力
  use cloudset, only: DensIce, &        ! 固相物質の密度 
    &                 NumAerosol, &     ! 質量あたりの凝結核数密度
    &                 RadiAerosol, &    ! 凝結核半径
    &                 Kd,     &         ! 大気の拡散係数
!    &                 SatPressA, &      ! Antoine の式の係数 A 
!    &                 SatPressB,  &     ! Antoine の式の係数 B
    &                 SatRatioCr        ! 臨界飽和比
!  use basicset, only: ss_CpBasicZ, &       ! 定圧比熱
!    &                 ss_PotTempBasicZ,  & ! 温位基本場
!    &                 ss_ExnerBasicZ,  &   ! 無次元圧力
!    &                 ss_DensBasicZ    &   ! 密度
!    &                 GasRDry              ! 気体定数
!  use BasicSet, only: xz_PotTempBasicZ,  & ! 温位基本場
  use basicset, only: xz_PotTempBasicZ,  & ! 温位基本場
    &                 xz_ExnerBasicZ,  &   ! 無次元圧力
    &                 xz_DensBasicZ,   &   ! 密度
    &                 GasRDry,         &    ! 気体定数
    &                 PressSfc,        &    ! 地表面圧力
    &                 CpDry                 ! 定圧比熱
  use ChemData, only: ChemData_Cp,     &   ! 定圧比熱
    &                 ChemData_SVapPress_AntoineA, & ! Antoine の式の係数 B
    &                 ChemData_SVapPress_AntoineB    ! Antoine の式の係数 B
!  use basicset, only: PressSfc             ! 地表面圧力

                                                                 !=end
  !==暗黙の型宣言を禁止
  implicit none
  
  !==Input
  real(8), intent(in)   :: xz_LatHeatPerMassNl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                        ! 単位質量あたりの凝結熱 [J/kg]
  real(8), intent(in)   :: xz_SatRatioNl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 飽和比 [1]
  real(8), intent(in)   :: xz_PotTempNl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                        ! 温位 [K]
  real(8), intent(in)   :: xz_ExnerNl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                        ! 無次元圧力 [1]
  real(8), intent(in)   :: xz_DensCloudNl(DimXMin:DimXMax, DimZMin:DimZMax)   
                                        ! 雲の密度   [kg/m^3]


  !==Output
  real(8), intent(out)  :: xz_MassCondNl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 凝結量 [kg/m^3 s]
   
  !==Work
  real(8)               :: xz_Rh(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 熱拡散による抵抗係数
  real(8)               :: xz_RadiCloud(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲粒の半径 [m]
  real, parameter       :: Pi = 3.1415926535897932385d0
                                        ! 円周率
  integer               :: i,k          ! ループ変数

  call BeginSub("MassCondense", &
    &           fmt="%c",        &
    &           c1="Calculate latent heat per unit mass.")

  
  call DbgMessage( fmt="*** Debug Message [MassCondense] *** xz_Rh ")

!write(*,*) xz_PotTempNl(1,:)

!  xz_Rh = xz_LatHeatPerMassNl**2.0d0 / Kd
!  xz_Rh = 1.0 / (xz_ExnerBasicZ + xz_ExnerNl)** 2.0d0
!   xz_Rh = 1.0 / (xz_PotTempBasicZ + xz_PotTempNl)**2.0d0
!   xz_Rh = xz_PotTempBasicZ
!   xz_Rh = xz_PotTempNl

  xz_Rh =  & 
    & xz_LatHeatPerMassNl**2.0d0 &
    & / (  &
    &      Kd * GasRDry * (xz_ExnerBasicZ + xz_ExnerNl)** 2.0d0 &
    &      * (xz_PotTempBasicZ + xz_PotTempNl)**2.0d0  &
    &   )
      ! 山下, 20080410, 
      ! L**2 / (k * R * Exner**2 * theta**2) という値を計算している. 
      ! 参考文献 : 北守修論の付録 A に導出過程が記述されている. 
      ! 但し, 北守修論本文の Rh の式 (3.4) には誤植があるので注意されたい. 
      ! 教科書として 1 次参照元はどれなのか調べる必要がありそうだ. 

  call DbgMessage( fmt="*** Debug Message [MassCondense] *** xz_RadiCloud ")


  xz_RadiCloud =   &
    & (  &
    &      3.0d0 * xz_DensCloudNl &
    &      / (4.0d0 * Pi * DensIce * NumAerosol * xz_DensBasicZ) & 
    &      + RadiAerosol**3.0d0  &
    &   )**(1.0d0 / 3.0d0)

  do k= DimZMin, DimZMax
   do i = DimXMin, DimXMax
!      if (xz_DensCloudNl(i,k) == 0 .and. &        
!        & xz_SatRatioNl(i,k) <= SatRatioCr) then
!        ! 雲密度 0 のときは過飽和度を越えないと凝結しない
!        xz_MassCondNl(i,k) = 0.0d0
!      else 
!        ! 凝結/蒸発量を計算する
!        xz_MassCondNl(i,k) = &
!          &   4.0d0 * Pi * xz_RadiCloud(i,k) &
!          & * NumAerosol * xz_DensBasicZ(i,k)  &
!          & / xz_Rh(i,k)  &
!          & * (xz_SatRatioNl(i,k) - 1.0d0)


      if (xz_SatRatioNl(i,k) - SatRatioCr > epsilon(0.0e0)) then
! 飽和比が定めた臨界飽和比よりも大きい場合
! CO2 の飽和比の式が地球大気で用いられる飽和比の式と等価なものなのか確認する必要がある. 
!      if (xz_SatRatioNl(i,k) / SatRatioCr > 1.0d0) then
        xz_MassCondNl(i,k) = &
         &   4.0d0 * Pi * xz_RadiCloud(i,k) &
         & * NumAerosol * xz_DensBasicZ(i,k)  &
         & / xz_Rh(i,k)  &
         & * (xz_SatRatioNl(i,k) - 1.0d0) 
      else if (xz_DensCloudNl(i,k) > 0.0d0 ) then
! 臨界飽和比を超えていないが, 雲は存在する場合
! 飽和比が 1 以上のときに凝結, 1 以下のときに蒸発. 
        xz_MassCondNl(i,k) = &
          &   4.0d0 * Pi * xz_RadiCloud(i,k) &
          & * NumAerosol * xz_DensBasicZ(i,k)  &
          & / xz_Rh(i,k)  &
          & * (xz_SatRatioNl(i,k) - 1.0d0)
      else
! 臨界飽和比を超えておらず, 雲も存在しない場合         
        xz_MassCondNl(i,k) = 0.0d0

!      if (xz_DensCloudNl(i,k) == 0 .and. &        
!        & xz_SatRatioNl(i,k) <= SatRatioCr) then
!        ! 雲密度 0 のときは過飽和度を越えないと凝結しない
!        xz_MassCondNl(i,k) = 0.0d0
!      else 
!        ! 凝結/蒸発量を計算する
!        xz_MassCondNl(i,k) = &
!          &   4.0d0 * Pi * xz_RadiCloud(i,k) &
!          & * NumAerosol * xz_DensBasicZ(i,k)  &
!          & / xz_Rh(i,k)  &
!          & * (xz_SatRatioNl(i,k) - 1.0d0)

      end if
    end do
  end do

  call EndSub("MassCondense")

end subroutine MassCondense
