!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2004. All rights reserved. !--------------------------------------------------------------------- !=begin != Subroutine MassCondense ! ! * Developer: KITAMORI Taichi, YAMASHITA Tatsuya ! * Version: $ ! * Tag Name: $Name: arare4-20120911 $ ! * Change History: ! !== Overview ! !凝結/蒸発量を計算する. ! !== Error Handling ! !== Known Bugs ! !== Note ! ! 拡散によって雲粒が成長する場合における成長方程式を解く. ! 蒸発量が存在する雲の量よりも多い場合, 蒸発量を存在する雲の量にする. ! 雲密度を時間発展(凝結部分)を計算. ! 雲が非常に小さくなったときに蒸発が起きなくなってしまっていたので, ! 条件分岐を書き直した. ! !== Future Plans ! !=end subroutine MassCondense( & & DelTime, & !(in) 時間間隔 & xz_LatHeatPerMassNl, & !(in) 単位質量あたりの潜熱 & xz_SatRatioNl, & !(in) 過飽和度 & xz_PotTempNl, & !(in) 温位 & xz_ExnerNl, & !(in) エクスナー関数 & xz_DensCloudNl, & !(inout) 雲密度 & xz_MassCondNl) !(out) 凝結量 !=begin !==Dependency use dc_trace, only: BeginSub, EndSub, DbgMessage use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax 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 StoreDensCloud, only: StoreDensCloudCond !=end !==暗黙の型宣言を禁止 implicit none !==Input real(8), intent(in) :: DelTime ! 時間ステップ [s] 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(inout) :: 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 = 5.0d-5 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_SatRatioNl(i,k) - SatRatioCr > epsilon(0.0d0)) then ! 飽和比が定めた臨界飽和比よりも大きい場合 ! CO2 の飽和比の式が地球大気で用いられる飽和比の式と等価なものなのか確認する必要がある. xz_MassCondNl(i,k) = & & max( - xz_DensCloudNl(i,k) / DelTime, & & 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) = & & max( - xz_DensCloudNl(i,k) / DelTime, & & 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 .and. xz_SatRatioNl(i,k) < 1.0d0 ) then ! 雲密度が閾値未満で, 飽和比 1.0 未満である場合に蒸発. xz_MassCondNl(i,k) = & & max( - xz_DensCloudNl(i,k) / DelTime, & & 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 ! xz_MassCondNl(i,k) = & ! & max( - xz_DensCloudNl(i,k) / DelTime, & ! & ( 0.5d0 * & ! & ( 1.0d0 + sign(1.0d0, xz_DensCloudNl(i,k) - epsilon(0.0d0)) ) & ! & ) * & ! & ( 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) ) end if ! call StoreDensCloudCond( xz_MassCondNl(i,k) ) xz_DensCloudNl(i,k) = & & xz_DensCloudNl(i,k) + DelTime * xz_MassCondNl(i,k) end do end do call StoreDensCloudCond( xz_MassCondNl ) call EndSub("MassCondense") end subroutine MassCondense