!---------------------------------------------------------------------
!     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(  &
  & ss_LatHeatPerMass_nl,  &  !(in)
  & ss_SupSatRatio_nl, &      !(in)
  & ss_PotTemp_nl,     &      !(in)
  & ss_Exner_nl,       &      !(in)
  & ss_DensCloud_nl,     &    !(in)
  & ss_MassCond_nl)           !(out)
                                                                 !=begin
  !==Dependency
  use dc_trace, only: BeginSub, EndSub  
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax
  use physset,  only: GasR,   &         ! 気体定数  
    &                 Pi,     &         ! 円周率
    &                 PressSfc          ! 地表面圧力
  use cloudset, only: DensIce, &        ! 固相物質の密度 
    &                 NumAerosol, &     ! 質量あたりの凝結核数密度
    &                 RadiAerosol, &    ! 凝結核半径
    &                 Kd,     &         ! 大気の拡散係数
    &                 SatPressA, &      ! Antoine の式の係数 A 
    &                 SatPressB,  &     ! Antoine の式の係数 B
    &                 SupSatRatioCr     ! 臨界過飽和度
  use basicset, only: ss_CpBasicZ, &       ! 定圧比熱
    &                 ss_PotTempBasicZ,  & ! 温位基本場
    &                 ss_ExnerBasicZ,  &   ! 無次元圧力
    &                 ss_DensBasicZ        ! 密度
                                                                 !=end
  !==暗黙の型宣言を禁止
  implicit none
  
  !==Input
  real(8), intent(in)   :: ss_LatHeatPerMass_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                        ! 単位質量あたりの凝結熱 [J/kg]
  real(8), intent(in)   :: ss_SupSatRatio_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 過飽和度 [1]
  real(8), intent(in)   :: ss_PotTemp_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                        ! 温位 [K]
  real(8), intent(in)   :: ss_Exner_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                        ! 無次元圧力 [1]
  real(8), intent(in)   :: ss_DensCloud_nl(DimXMin:DimXMax, DimZMin:DimZMax)   
                                        ! 雲の密度   [kg/m^3]

  !==Output
  real(8), intent(out)  :: ss_MassCond_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 凝結量 [kg/m^3 s]
   
  !==Work
  real(8)               :: ss_Rh(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 熱拡散による抵抗係数
  real(8)               :: ss_RadiCloud(DimXMin:DimXMax, DimZMin:DimZMax)
                                        ! 雲粒の半径 [m]
  integer               :: i,k          ! ループ変数

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


!  ss_SupSatRatio =  &
!    &  PressSfc * (ss_ExnerBasicZ + ss_Exner_nl)**(ss_CpBasicZ / GasR) &
!    &  / (  exp(  &
!    &             SatPressA &
!    &           - SatPressB  &
!    &             / (ss_ExnerBasicZ + ss_Exner_nl) &
!    &             / (ss_PotTempBasicZ + ss_PotTemp_nl) &
!    &          )  &
!    &    ) &
!    &  - 1.0d0

  ss_Rh =  & 
    & ss_LatHeatPerMass_nl**2.0d0 &
    & / (  &
    &      Kd * GasR * (ss_ExnerBasicZ + ss_Exner_nl)** 2.0d0 &
    &      * (ss_PotTempBasicZ + ss_PotTemp_nl)**2.0d0  &
    &   ) 

  ss_RadiCloud =   &
    & (  &
    &      3.0d0 * ss_DensCloud_nl &
    &      / (4.0d0 * Pi * DensIce * NumAerosol * ss_DensBasicZ) & 
    &      + RadiAerosol**3.0d0  &
    &   )**(1.0d0 / 3.0d0) 

  do k= DimZMin, DimZMax
    do i = DimXMin, DimXMax
      if (ss_DensCloud_nl(i,k) == 0 .and. &        
        & ss_SupSatRatio_nl(i,k) < SupSatRatioCr) then
        ! 雲密度 0 のときは過飽和度を越えないと凝結しない
        ss_MassCond_nl(i,k) = 0.0d0
      else
        ! 凝結/蒸発量を計算する
        ss_MassCond_nl(i,k) = &
          &   4.0d0 * Pi * ss_RadiCloud(i,k) &
          & * NumAerosol * ss_DensBasicZ(i,k)  &
          & / ss_Rh(i,k)  &
          & * ss_SupSatRatio_nl(i,k)
      end if
    end do
  end do

  call EndSub("MassCondense")

end subroutine MassCondense
