subroutine CloudT1993base( xyz_Press, xyr_Press, xyz_VirTemp, xyz_DQCloudWaterDtCum, xyz_MoistConvDetTend, xyz_OMG, xyz_MoistConvSubsidMassFlux, xyz_DTempDtPhy, xyz_Temp, xyz_QH2OVap, xyz_QCloudWater, xyz_CloudCover, xy_SurfRainFlux, xy_SurfSnowFlux )
! USE statements
!
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! 物理・数学定数設定
! Physical and mathematical constants settings
!
use constants0, only: PI ! $ \pi $ .
! 円周率. Circular constant
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, CpDry, GasRDry, GasRWet, LatentHeat
! $ L $ [J kg-1] .
! 凝結の潜熱.
! Latent heat of condensation
! 飽和比湿の算出
! Evaluate saturation specific humidity
!
use saturate, only: xyz_CalcQVapSat, xyz_CalcDQVapSatDTemp
! 大規模凝結 (非対流性凝結)
! Large scale condensation
!
use lscond, only: LScaleCond
! 雲関系ルーチン
! Cloud-related routines
!
use cloud_utils, only : CloudUtilsCalcPRCPKeyLLTemp
real(DP), intent(in ) :: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(in ) :: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_DQCloudWaterDtCum (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_MoistConvDetTend (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_OMG (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_MoistConvSubsidMassFlux(0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_DTempDtPhy (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(inout) :: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(inout) :: xyz_QH2OVap (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(inout) :: xyz_QCloudWater (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(inout) :: xyz_CloudCover (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(out) :: xy_SurfRainFlux (0:imax-1, 1:jmax)
real(DP), intent(out) :: xy_SurfSnowFlux (0:imax-1, 1:jmax)
real(DP) :: xyz_QCloudWaterB (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_QH2OVapSat (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DQH2OVapSatDPress(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DQH2OVapSatDTemp (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DQH2OVapSatDt (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_ZeroOneDQsDt (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DelCloudCoverStr(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactA (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactA1 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactA2 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactB (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactC (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactC1 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactC2 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactD (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactE (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactE1 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactE2 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_Rain (0:imax-1, 1:jmax)
real(DP) :: xyz_Rain (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_RainConvFactor (0:imax-1, 1:jmax)
real(DP) :: xy_FactCo (0:imax-1, 1:jmax)
real(DP) :: xy_FactBF (0:imax-1, 1:jmax)
real(DP) :: xy_QCloudWaterConvThreshold(0:imax-1, 1:jmax)
real(DP) :: xyz_DTempDtLSC (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DQVapDtLSC (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_RainLSC (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_RainLSC (0:imax-1, 1:jmax)
real(DP), parameter :: MixCoef = 1.0d-6
real(DP), parameter :: QCloudWaterEvapThreshold = 1.0d-10
!!$ real(DP), parameter :: RHThreshold = 0.999_DP
real(DP), parameter :: RHThreshold = 1.0_DP - 1.0d-10
! Values below are obtained from Sundqvist et al. (1989), but some of
! those are arbitrarily selected.
real(DP), parameter :: C0 = 1.0d-3
real(DP), parameter :: C1 = 300.0_DP
real(DP), parameter :: C2 = 0.5_DP
real(DP), parameter :: QCloudWaterConvThreshold = 4.0d-4
! Parameters for evaporation of rain
real(DP), parameter :: DensWater = 1.0d3
! rho_w
! Values below are from Kessler (1969)
real(DP), parameter :: RainFallVelFactor = 130.0d0
! K
real(DP), parameter :: MedianDiameterFactor = 3.67d0
! C'
real(DP), parameter :: RainDistFactor = 1.0d7
! N0
real(DP), parameter :: RainEvapRatUnitDiamFactor = 2.24d3
! C
real(DP), parameter :: H2OVapDiffCoef = 1.0d-5
! Kd
real(DP) :: Dens0
! rho_0
real(DP) :: V00
! V_{00}
real(DP) :: RainEvapFactor
real(DP) :: xyz_Dens (0:imax-1, 1:jmax, 1:kmax)
! rho
real(DP) :: xy_DensRain (0:imax-1, 1:jmax)
! (rho q_r)
real(DP) :: xy_RainArea (0:imax-1, 1:jmax)
! a_p
real(DP) :: xy_RainEvapArea (0:imax-1, 1:jmax)
! A = max( a_p - a, 0 )
real(DP) :: xyz_RainEvapRate(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_DelRain (0:imax-1, 1:jmax)
real(DP) :: DelQH2OVap
real(DP) :: RainFallVel
integer :: i
integer :: j
integer :: k
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. cloud_T1993base_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! Parameters for evaporation of rain
Dens0 = 1013.0d2 / ( GasRDry * 300.0_DP )
V00 = RainFallVelFactor * sqrt( MedianDiameterFactor ) / ( PI * DensWater * RainDistFactor )**(1.0d0/8.0d0)
RainEvapFactor = RainEvapRatUnitDiamFactor * 1.429624558860304d0 * H2OVapDiffCoef * RainDistFactor**(7.0d0/20.0d0) / ( PI * DensWater )**(13.0d0/20.0d0)
! Values for evaporation of rain
xyz_Dens = xyz_Press / ( GasRDry * xyz_VirTemp )
! Save cloud water amount
!
xyz_QCloudWaterB = xyz_QCloudWater
xyz_QH2OVapSat = xyz_CalcQVapSat( xyz_Temp, xyz_Press )
xyz_DQH2OVapSatDTemp = xyz_CalcDQVapSatDTemp( xyz_Temp, xyz_QH2OVapSat )
xyz_DQH2OVapSatDPress = xyz_QH2OVapSat / xyz_Press * ( LatentHeat * GasRDry * xyz_VirTemp - CpDry * GasRWet * xyz_Temp**2 ) / ( xyz_CloudCover * LatentHeat**2 * xyz_QH2OVapSat + CpDry * GasRWet * xyz_Temp**2 )
xyz_DQH2OVapSatDt = xyz_DQH2OVapSatDPress * ( xyz_OMG + Grav * xyz_MoistConvSubsidMassFlux ) + xyz_DQH2OVapSatDTemp / ( 1.0_DP + xyz_CloudCover * LatentHeat / CpDry * xyz_DQH2OVapSatDTemp ) * xyz_DTempDtPhy
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_DQH2OVapSatDt(i,j,k) < 0.0_DP ) then
xyz_ZeroOneDQsDt(i,j,k) = 1.0_DP
else
xyz_ZeroOneDQsDt(i,j,k) = 0.0_DP
end if
end do
end do
end do
! Rain at the surface
xy_Rain = 0.0_DP
! Rain area
xy_RainArea = 0.0_DP
k_loop : do k = kmax, 1, -1
! Evaporation of rain
!
if ( k == kmax ) then
xy_RainArea = 0.0_DP
else
do j = 1, jmax
do i = 0, imax-1
if ( xyz_Rain(i,j,k+1) > 0.0_DP ) then
if ( xyz_CloudCover(i,j,k+1) > xy_RainArea(i,j) ) then
xy_RainArea(i,j) = xyz_CloudCover(i,j,k+1)
end if
end if
end do
end do
end if
xy_DensRain = ( xy_Rain / ( xy_RainArea + 1.0d-10 ) / ( V00 * sqrt( Dens0 / xyz_Dens(:,:,k) ) ) )**(8.0d0/9.0d0)
xy_RainEvapArea = max( xy_RainArea - xyz_CloudCover(:,:,k), 0.0_DP )
xyz_RainEvapRate(:,:,k) = xyz_Dens(:,:,k) * xy_RainEvapArea * RainEvapFactor * max( xyz_QH2OVapSat(:,:,k) - xyz_QH2OVap(:,:,k), 0.0_DP ) * xy_DensRain**(13.0d0/20.0d0)
do j = 1, jmax
do i = 0, imax-1
RainFallVel = V00 * sqrt( Dens0 / xyz_Dens(i,j,k) ) * xy_DensRain(i,j)**(1.0d0/8.0d0)
if ( xy_RainArea(i,j) * RainFallVel * xyz_RainEvapRate(i,j,k) * 2.0_DP * DelTime > xy_Rain(i,j) ) then
xyz_RainEvapRate(i,j,k) = xy_Rain(i,j) / ( ( xy_RainArea(i,j) + 1.0d-10 ) * RainFallVel * 2.0_DP * DelTime )
xy_DelRain(i,j) = - xy_Rain(i,j)
else
xy_DelRain(i,j) = - xy_RainArea(i,j) * RainFallVel * xyz_RainEvapRate(i,j,k) * ( 2.0_DP * DelTime )
end if
xy_Rain(i,j) = xy_Rain(i,j) + xy_DelRain(i,j)
!!$ xyz_QH2OVap(i,j,k) = xyz_QH2OVap(i,j,k) &
!!$ & + xyz_RainEvapRate(i,j,k) * ( 2.0_DP * DelTime ) &
!!$ & / xyz_Dens(i,j,k)
!!$ xyz_Temp(i,j,k) = xyz_Temp(i,j,k) &
!!$ & - xyz_RainEvapRate(i,j,k) * ( 2.0_DP * DelTime ) &
!!$ & / xyz_Dens(i,j,k) &
!!$ & * LatentHeat / CpDry
! DelRain = DelQH2OVap * DelPress / Grav / ( 2.0_DP * DelTime )
DelQH2OVap = - xy_DelRain(i,j) * ( 2.0_DP * DelTime ) * Grav / ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) )
xyz_QH2OVap(i,j,k) = xyz_QH2OVap(i,j,k) + DelQH2OVap
xyz_Temp(i,j,k) = xyz_Temp(i,j,k) - DelQH2OVap * LatentHeat / CpDry
end do
end do
xyz_FactC1(:,:,k) = xyz_MoistConvDetTend(:,:,k)
do j = 1, jmax
do i = 0, imax-1
if ( xyz_QH2OVap(i,j,k) < RHThreshold * xyz_QH2OVapSat(i,j,k) ) then
xyz_FactC2(i,j,k) = - ( 1.0_DP - xyz_CloudCover(i,j,k) ) / ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) * xyz_DQH2OVapSatDt(i,j,k) * xyz_ZeroOneDQsDt(i,j,k)
else
!!$ xyz_FactC2(i,j,k) = 0.0_DP
xyz_FactC2(i,j,k) = 1.0_DP / ( 2.0_DP * DelTime )
end if
end do
end do
xyz_FactC(:,:,k) = xyz_FactC1(:,:,k) + xyz_FactC2(:,:,k)
do j = 1, jmax
do i = 0, imax-1
if ( xyz_QH2OVapSat(i,j,k) > xyz_QH2OVap(i,j,k) ) then
xyz_FactD(i,j,k) = 2.0_DP * xyz_CloudCover(i,j,k) * MixCoef * ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) / ( xyz_QCloudWater(i,j,k) + 1.0d-100 )
else
xyz_FactD(i,j,k) = 0.0_DP
end if
end do
end do
!
do j = 1, jmax
do i = 0, imax-1
if ( xyz_QH2OVap(i,j,k) < RHThreshold * xyz_QH2OVapSat(i,j,k) ) then
xyz_FactE1(i,j,k) = ( 1.0_DP - xyz_CloudCover(i,j,k) )**2 / ( 2.0_DP * ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) ) * xyz_DQH2OVapSatDt(i,j,k) * xyz_ZeroOneDQsDt(i,j,k)
else
xyz_FactE1(i,j,k) = 0.0_DP
end if
end do
end do
do j = 1, jmax
do i = 0, imax-1
if ( xyz_QH2OVapSat(i,j,k) > xyz_QH2OVap(i,j,k) ) then
xyz_FactE2(i,j,k) = + xyz_CloudCover(i,j,k)**2 * MixCoef * ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) / ( xyz_QCloudWater(i,j,k) + 1.0d-100 )
else
xyz_FactE2(i,j,k) = 0.0_DP
end if
end do
end do
xyz_FactE(:,:,k) = xyz_FactE1(:,:,k) + xyz_FactE2(:,:,k)
!
xyz_DelCloudCoverStr(:,:,k) = xyz_FactC2(:,:,k) * 2.0_DP * DelTime - xyz_FactC2(:,:,k) / ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) + 1.0d-100 ) * ( xyz_FactC(:,:,k) * 2.0_DP * DelTime + ( xyz_CloudCover(:,:,k) - xyz_FactC(:,:,k) / ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) + 1.0d-100 ) ) * ( 1.0_DP - exp( - ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) ) * 2.0_DP * DelTime ) ) )
!
xyz_FactA1(:,:,k) = xyz_DQCloudWaterDtCum(:,:,k)
xyz_FactA2(:,:,k) = - xyz_CloudCover(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) - xyz_DelCloudCoverStr(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * xyz_ZeroOneDQsDt(:,:,k) / 2.0_DP - xyz_CloudCover(:,:,k) * MixCoef * ( xyz_QH2OVapSat(:,:,k) - xyz_QH2OVap(:,:,k) )
! The value of xyz_FactA2 is checked, and is updated.
do j = 1, jmax
do i = 0, imax-1
if ( xyz_FactA2(i,j,k) * 2.0_DP * DelTime > xyz_QH2OVap(i,j,k) ) then
xyz_FactA2(i,j,k) = xyz_QH2OVap(i,j,k) / ( 2.0_DP * DelTime )
end if
end do
end do
xyz_FactA(:,:,k) = xyz_FactA1(:,:,k) + xyz_FactA2(:,:,k)
!!$ xy_RainConvFactor = 1.0_DP / ( CloudLifeTime + 1.0d-100 )
!
xy_FactCo = 1.0_DP + C1 * sqrt( xy_Rain )
xy_FactBF = 1.0_DP + C2 * sqrt( max( 268.0_DP - xyz_Temp(:,:,k), 0.0_DP ) )
xy_QCloudWaterConvThreshold = QCloudWaterConvThreshold / ( xy_FactCo * xy_FactBF )
xy_RainConvFactor = C0 * xy_FactCo * xy_FactBF * ( 1.0_DP - exp( - ( xyz_QCloudWater(:,:,k) / ( ( xyz_CloudCover(:,:,k) + 1.0d-100 ) * xy_QCloudWaterConvThreshold ) )**2 ) )
!
xyz_FactB(:,:,k) = xy_RainConvFactor
!
do j = 1, jmax
do i = 0, imax-1
if ( xyz_FactB(i,j,k) < 1.0_DP / 1.0d10 ) then
xyz_FactB(i,j,k) = 1.0_DP / 1.0d10
end if
end do
end do
! Values at next time step are calculated.
!
xyz_QCloudWater(:,:,k) = xyz_QCloudWater(:,:,k) * exp( - xyz_FactB(:,:,k) * 2.0_DP * DelTime ) + xyz_FactA(:,:,k) / xyz_FactB(:,:,k) * ( 1.0_DP - exp( - xyz_FactB(:,:,k) * 2.0_DP * DelTime ) )
! The value of cloud water amount is checked, and xyz_FactA
! is updated.
do j = 1, jmax
do i = 0, imax-1
if ( xyz_QCloudWater(i,j,k) < 0.0_DP ) then
xyz_FactA2(i,j,k) = - xyz_FactA1(i,j,k) - xyz_FactB(i,j,k) * xyz_QCloudWaterB(i,j,k) * exp( - xyz_FactB(i,j,k) * 2.0_DP * DelTime ) / ( 1.0_DP - exp( - xyz_FactB(i,j,k) * 2.0_DP * DelTime ) )
xyz_QCloudWater(i,j,k) = 0.0_DP
end if
end do
end do
xyz_FactA(:,:,k) = xyz_FactA1(:,:,k) + xyz_FactA2(:,:,k)
!
xyz_CloudCover(:,:,k) = xyz_CloudCover(:,:,k) * exp( - ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) ) * 2.0_DP * DelTime ) + ( xyz_FactC(:,:,k) + xyz_FactE(:,:,k) ) / ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) + 1.0d-100 ) * ( 1.0_DP - exp( - ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) ) * 2.0_DP * DelTime ) )
!
xyz_QH2OVap(:,:,k) = xyz_QH2OVap(:,:,k) - xyz_FactA2(:,:,k) * 2.0_DP * DelTime
!
xyz_Temp(:,:,k) = xyz_Temp(:,:,k) + xyz_FactA2(:,:,k) * 2.0_DP * DelTime * LatentHeat / CpDry
! Rain
!
xyz_Rain(:,:,k) = xyz_FactA(:,:,k) * 2.0_DP * DelTime + xyz_QCloudWaterB(:,:,k) * ( 1.0_DP - exp( - xyz_FactB(:,:,k) * 2.0_DP * DelTime ) ) - xyz_FactA(:,:,k) / xyz_FactB(:,:,k) * ( 1.0_DP - exp( - xyz_FactB(:,:,k) * 2.0_DP * DelTime ) )
xyz_Rain(:,:,k) = xyz_Rain(:,:,k) / ( 2.0_DP * DelTime )
! Rain at the surface
xy_Rain = xy_Rain + xyz_Rain(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
! Evaporation
!
do j = 1, jmax
do i = 0, imax-1
if ( xyz_QCloudWater(i,j,k) < QCloudWaterEvapThreshold ) then
xyz_CloudCover(i,j,k) = 0.0_DP
end if
end do
end do
! Cloud cover is restricted.
do j = 1, jmax
do i = 0, imax-1
if ( xyz_CloudCover(i,j,k) > 1.0_DP ) then
xyz_CloudCover(i,j,k) = 1.0_DP
else if ( xyz_CloudCover(i,j,k) < 0.0_DP ) then
xyz_CloudCover(i,j,k) = 0.0_DP
end if
end do
end do
! Check values
do j = 1, jmax
do i = 0, imax-1
if ( xyz_QCloudWater(i,j,k) < 0.0_DP ) then
write( 6, * ) 'QCloudWater is negative', i, j, k, xyz_QCloudWater(i,j,k)
end if
end do
end do
end do k_loop
xyz_DTempDtLSC = 0.0_DP
xyz_DQVapDtLSC = 0.0_DP
! 大規模凝結 (非対流性凝結) (Manabe, 1965)
! Large scale condensation (non-convective condensation) (Manabe, 1965)
!
call LScaleCond( xyz_Temp, xyz_QH2OVap, xyz_DTempDtLSC, xyz_DQVapDtLSC, xyz_Press, xyr_Press, xyz_RainLSC )
xy_RainLSC = 0.0_DP
do k = kmax, 1, -1
xy_RainLSC = xy_RainLSC + xyz_RainLSC(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
end do
xy_Rain = xy_Rain + xy_RainLSC
call CloudUtilsCalcPRCPKeyLLTemp( xyz_Temp, xy_Rain, xy_SurfRainFlux, xy_SurfSnowFlux )
end subroutine CloudT1993base