subroutine MassCondense( xz_LatHeatPerMassNl, xz_SatRatioNl, xz_PotTempNl, xz_ExnerNl, xz_DensCloudNl, xz_MassCondNl) !(out) 凝結量
!=begin
!==Dependency
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, 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, 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
! 円周率
! real(8) :: DensCloudThreshold = 1.0d-5
! real(8) :: DensCloudThreshold = 5.0d-6
real(8) :: DensCloudThreshold = 1.0d-6
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 * 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
! if (xz_SatRatioNl(i,k) / SatRatioCr > 1.0d0) then
! 飽和比が定めた臨界飽和比よりも大きい場合
! CO2 の飽和比の式が地球大気で用いられる飽和比の式と等価なものなのか確認する必要がある.
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) > DensCloudThreshold ) 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
! 臨界飽和比を超えておらず, 雲密度が閾値未満の場合
! 雲密度がゼロではなく, 飽和比 1.0 未満のときに限って蒸発.
! それ以外の場合には何も起きない.
! xz_MassCondNl(i,k) = 0.0d0
xz_MassCondNl(i,k) = ( -0.5d0 * ( -1.0d0 + sign(1.0d0, - xz_DensCloudNl(i,k)) ) ) * ( -0.5d0 * ( -1.0d0 + sign(1.0d0, xz_SatRatioNl(i,k) - 1.0d0) ) ) * 4.0d0 * Pi * xz_RadiCloud(i,k) * NumAerosol * xz_DensBasicZ(i,k) / xz_Rh(i,k) * (xz_SatRatioNl(i,k) - 1.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