| Class | cloud_T1993base |
| In: |
radiation/cloud_T1993base.f90
|
Note that Japanese and English are described in parallel.
簡単雲モデルによる雲の計算.
In this module, the amount of cloud is calculated by use of a simple cloud model.
| !$ ! RadiationFluxDennouAGCM : | 放射フラックスの計算 |
| !$ ! ———— : | ———— |
| !$ ! RadiationFluxDennouAGCM : | Calculate radiation flux |
| Subroutine : | |
| 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(in ) |
| 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(inout) |
| xy_SurfRainFlux(0:imax-1, 1:jmax) : | real(DP), intent(out) |
| xy_SurfSnowFlux(0:imax-1, 1:jmax) : | real(DP), intent(out) |
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
! 簡単雲モデル
! Simple cloud model
!
use cloud_simple, only : CloudSimpleCalcPRCPKeyLLTemp
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_ZeroOneCloudProd(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_ZeroOneCloudLoss(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_QCloudWatConvThreshold(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) :: RainConvFactor0
real(DP), parameter :: C1 = 300.0_DP
real(DP), parameter :: C2 = 0.5_DP
! 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
real(DP) :: xyz_Sigma (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudProdRHThreshold(0:imax-1, 1:jmax, 1:kmax)
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
! set zero-one flag
do k = 1, kmax
xyz_Sigma(:,:,k) = xyz_Press(:,:,k) / xyr_Press(:,:,0)
end do
xyz_CloudProdRHThreshold = RHThresholdCrtl + ( 1.0_DP - xyz_QH2OVap / xyz_QH2OVapSat ) * ( max( ( xyz_Sigma - RHThresholdSigmaMin ), 0.0_DP ) / ( xyz_Sigma - RHThresholdSigmaMin ) )**RHThresholdOrd
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_DQH2OVapSatDt(i,j,k) < 0.0_DP ) then
if ( xyz_QH2OVap(i,j,k) >= xyz_CloudProdRHThreshold(i,j,k) * xyz_QH2OVapSat(i,j,k) ) then
xyz_ZeroOneCloudProd(i,j,k) = 1.0_DP
else
xyz_ZeroOneCloudProd(i,j,k) = 0.0_DP
end if
xyz_ZeroOneCloudLoss(i,j,k) = 1.0_DP
else
xyz_ZeroOneCloudProd(i,j,k) = 0.0_DP
xyz_ZeroOneCloudLoss(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) * 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_ZeroOneCloudProd(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_ZeroOneCloudProd(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_FactE1(:,:,k) * 2.0_DP * DelTime
!
xyz_FactA1(:,:,k) = xyz_DQCloudWaterDtCum(:,:,k)
xyz_FactA2(:,:,k) = - xyz_CloudCover(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * xyz_ZeroOneCloudProd(:,:,k) - xyz_DelCloudCoverStr(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * xyz_ZeroOneCloudProd(:,:,k) / 2.0_DP - xyz_CloudCover(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * ( 1.0_DP - xyz_ZeroOneCloudLoss(:,:,k) ) - 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 ) * ( 1.0_DP - 1.0d-14 )
end if
end do
end do
xyz_FactA(:,:,k) = xyz_FactA1(:,:,k) + xyz_FactA2(:,:,k)
!!$ xy_RainConvFactor = 1.0_DP / ( CloudLifeTime0 + 1.0d-100 )
!
xy_FactCo = 1.0_DP + C1 * sqrt( xy_Rain )
! Factor for Bergeron-Findeisen effect
! Sundqvist et al. (1989)
!!$ xy_FactBF = 1.0_DP &
!!$ & + C2 * sqrt( max( 268.0_DP - xyz_Temp(:,:,k), 0.0_DP ) )
! An original equation following that by IFS CY38r1
! (p.82 of http://www.ecmwf.int/research/ifsdocs/CY38r1/IFSPart4.pdf)
!!$ xy_FactBF = 1.0_DP &
!!$ & + C2 * sqrt( min( &
!!$ & max( 268.0_DP - xyz_Temp(:,:,k), 0.0_DP ), &
!!$ & max( 268.0_DP - TempBFEffectSat, 0.0_DP ) &
!!$ & ) &
!!$ & )
! Constant (no effect)
xy_FactBF = 1.0_DP
!
RainConvFactor0 = 1.0_DP / ( CloudLifeTime0 + 1.0d-100 )
!!$ xy_QCloudWatConvThreshold = &
!!$ & QCloudWatEffConv0 / ( xy_FactCo * xy_FactBF )
xy_QCloudWatConvThreshold = ( QCloudWatEffConv0 + 1.0e-100_DP ) / ( xy_FactCo * xy_FactBF )
xy_RainConvFactor = RainConvFactor0 * xy_FactCo * xy_FactBF * ( 1.0_DP - exp( - ( xyz_QCloudWater(:,:,k) / ( ( xyz_CloudCover(:,:,k) + 1.0d-100 ) * xy_QCloudWatConvThreshold ) )**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_QH2OVap(i,j,k) < 0.0_DP ) then
write( 6, * ) 'QH2OVap is negative', i, j, k, xyz_QH2OVap(i,j,k)
end if
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_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 CloudSimpleCalcPRCPKeyLLTemp( xyz_Temp, xy_Rain, xy_SurfRainFlux, xy_SurfSnowFlux )
end subroutine CloudT1993base
| Subroutine : | |
| ArgFlagSnow : | logical, intent(in) |
This procedure input/output NAMELIST#cloud_T1993base_nml .
subroutine CloudT1993baseInit( ArgFlagSnow )
! ファイル入出力補助
! File I/O support
!
use dc_iounit, only: FileOpen
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoAddVariable
! 飽和比湿の算出
! Evaluate saturation specific humidity
!
use saturate, only: SaturateInit
! 大規模凝結 (非対流性凝結)
! Large scale condensation (non-convective condensation)
!
use lscond, only : LScaleCondInit
! 宣言文 ; Declaration statements
!
logical, intent(in) :: ArgFlagSnow
integer:: unit_nml ! NAMELIST ファイルオープン用装置番号.
! Unit number for NAMELIST file open
integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT.
! IOSTAT of NAMELIST read
! NAMELIST 変数群
! NAMELIST group name
!
namelist /cloud_T1993base_nml/ RHThresholdCrtl, RHThresholdSigmaMin, RHThresholdOrd, CloudLifeTime0, CloudWatLifeTime0, CloudIceLifeTime0, QCloudWatEffConv0, QCloudIceEffConv0, TempBFEffectSat, PRCPArea, PRCPEvapArea, PRCPColFactor
!
! デフォルト値については初期化手続 "cloud_T1993base#CloudT1993baseInit"
! のソースコードを参照のこと.
!
! Refer to source codes in the initialization procedure
! "cloud_T1993base#CloudT1993baseInit" for the default values.
!
! 実行文 ; Executable statement
!
if ( cloud_T1993base_inited ) return
FlagSnow = ArgFlagSnow
! デフォルト値の設定
! Default values settings
!
RHThresholdCrtl = 0.5_DP
!!$ RHThresholdCrtl = 0.8_DP ! ECMWF IFS
RHThresholdSigmaMin = 1.0_DP
!!$ RHThresholdSigmaMin = 0.8_DP ! ECMWF IFS
RHThresholdOrd = 2.0_DP ! ECMWF IFS
CloudLifeTime0 = 1000.0_DP
CloudWatLifeTime0 = 1000.0_DP
CloudIceLifeTime0 = 1000.0_DP
QCloudWatEffConv0 = 1.0e-4_DP
QCloudIceEffConv0 = 1.0e-3_DP
TempBFEffectSat = 250.0_DP
! This value follows that by IFS CY38r1
! (p.82 of http://www.ecmwf.int/research/ifsdocs/CY38r1/IFSPart4.pdf)
! Actually, IFS CY38r1 uses 250.16 K.
PRCPArea = 1.0_DP
!!$ PRCPArea = 0.5_DP
PRCPEvapArea = 1.0_DP
PRCPColFactor = 300.0_DP
! This value comes from Sundqvist et al. (1989)
! IFS CY38r1 uses 100, but in addition ,
! precipitation area has to be considered.
! NAMELIST の読み込み
! NAMELIST is input
!
if ( trim(namelist_filename) /= '' ) then
call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)
rewind( unit_nml )
read( unit_nml, nml = cloud_T1993base_nml, iostat = iostat_nml ) ! (out)
close( unit_nml )
call NmlutilMsg( iostat_nml, module_name ) ! (in)
end if
! Initialization of modules used in this module
!
! 飽和比湿の算出
! Evaluate saturation specific humidity
!
call SaturateInit
! 大規模凝結 (非対流性凝結) (Manabe, 1965)
! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
!
call LScaleCondInit( FlagSnow )
! ヒストリデータ出力のためのへの変数登録
! Register of variables for history data output
!
!!$ call HistoryAutoAddVariable( 'EffCloudCover', &
!!$ & (/ 'lon ', 'lat ', 'time' /), &
!!$ & 'effective cloud cover', '1' )
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, 'RHThresholdCrtl = %f', d = (/ RHThresholdCrtl /) )
call MessageNotify( 'M', module_name, 'RHThresholdSigmaMin = %f', d = (/ RHThresholdSigmaMin /) )
call MessageNotify( 'M', module_name, 'RHThresholdOrd = %f', d = (/ RHThresholdOrd /) )
call MessageNotify( 'M', module_name, 'CloudLifeTime0 = %f', d = (/ CloudLifeTime0 /) )
call MessageNotify( 'M', module_name, 'CloudIceLifeTime0 = %f', d = (/ CloudIceLifeTime0 /) )
call MessageNotify( 'M', module_name, 'CloudWatLifeTime0 = %f', d = (/ CloudWatLifeTime0 /) )
call MessageNotify( 'M', module_name, 'QCloudWatEffConv0 = %f', d = (/ QCloudWatEffConv0 /) )
call MessageNotify( 'M', module_name, 'QCloudIceEffConv0 = %f', d = (/ QCloudIceEffConv0 /) )
call MessageNotify( 'M', module_name, 'TempBFEffectSat = %f', d = (/ TempBFEffectSat /) )
call MessageNotify( 'M', module_name, 'PRCPArea = %f', d = (/ PRCPArea /) )
call MessageNotify( 'M', module_name, 'PRCPEvapArea = %f', d = (/ PRCPEvapArea /) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
cloud_T1993base_inited = .true.
end subroutine CloudT1993baseInit
| Subroutine : | |
| 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_DQCloudWatDtCum(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_DQCloudIceDtCum(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(in ) |
| 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_QCloudWat(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout) |
| xyz_QCloudIce(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout) |
| xyz_CloudCover(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout) |
| xy_SurfRainFlux(0:imax-1, 1:jmax) : | real(DP), intent(out) |
| xy_SurfSnowFlux(0:imax-1, 1:jmax) : | real(DP), intent(out) |
subroutine CloudT1993baseWithIce( xyz_Press, xyr_Press, xyz_VirTemp, xyz_DQCloudWatDtCum, xyz_DQCloudIceDtCum, xyz_MoistConvDetTend, xyz_OMG, xyz_MoistConvSubsidMassFlux, xyz_DTempDtPhy, xyz_Temp, xyz_QH2OVap, xyz_QCloudWat, xyz_QCloudIce, 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, LatentHeatFusion, EpsV
! $ \epsilon_v $ .
! 水蒸気分子量比.
! Molecular weight of water vapor
! 飽和比湿の算出
! Evaluate saturation specific humidity
!
use saturate, only: xyz_CalcQVapSat , xyz_CalcDQVapSatDTemp
!!$ & xyz_CalcQVapSatOnLiq, xyz_CalcDQVapSatOnLiqDTemp, &
!!$ & xyz_CalcQVapSatOnSol, xyz_CalcDQVapSatOnSolDTemp
! 飽和比湿の算出
! Evaluate saturation specific humidity
!
use saturate, only : SaturateWatFraction
! 大規模凝結 (非対流性凝結)
! Large scale condensation
!
use lscond, only: LScaleCond1Grid
! 雲関系ルーチン
! Cloud-related routines
!
use cloud_utils, only : CloudUtilsPRCPStepPC1Grid, CloudUtilsPRCPEvap1Grid, CloudUtilConsChk
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_DQCloudWatDtCum (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_DQCloudIceDtCum (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_QCloudWat (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(inout) :: xyz_QCloudIce (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)
! Local variables
!
real(DP) :: xyz_TempB (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_QH2OVapB (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_QCloudWatB (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_QCloudIceB (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: LatentHeatSubl
real(DP) :: xyz_WatFrac (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_IceFrac (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_QH2OVapSat (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DQH2OVapSatDTemp (0:imax-1, 1:jmax, 1:kmax)
!!$ real(DP) :: xyz_QH2OVapSatOnLiq (0:imax-1, 1:jmax, 1:kmax)
!!$ real(DP) :: xyz_QH2OVapSatOnIce (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DQH2OVapSatDPressDenom(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DQH2OVapSatDPress(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DQH2OVapSatDt (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_ZOCloudProd (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_ZOCloudLoss (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DelCloudCoverStr(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactAl (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactAs (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactAl1 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactAs1 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactA20 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactAl2 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactAs2 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactBl (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_FactBs (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_Snow (0:imax-1, 1:jmax)
real(DP) :: xyz_Snow (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_FactCo (0:imax-1, 1:jmax)
real(DP) :: xy_FactBF (0:imax-1, 1:jmax)
real(DP) :: xy_QCloudWatConvThreshold(0:imax-1, 1:jmax)
real(DP) :: xy_QCloudIceConvThreshold(0:imax-1, 1:jmax)
real(DP) :: xy_FactIceTempDep (0:imax-1, 1:jmax)
real(DP) :: xy_FactConvThreshold (0:imax-1, 1:jmax)
real(DP) :: xy_RainConvFactor (0:imax-1, 1:jmax)
real(DP) :: xy_SnowConvFactor (0:imax-1, 1:jmax)
real(DP) :: xyz_DQH2OLiqDtLSC (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DQH2OSolDtLSC (0:imax-1, 1:jmax, 1:kmax)
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) :: RainConvFactor0
real(DP) :: SnowConvFactor0
real(DP), parameter :: C2 = 0.5_DP
! 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) :: xyz_DelMass( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: aaa_QH2OVapSat(1,1,1)
real(DP) :: QH2OVapSat
real(DP) :: QCloudWatTentative
real(DP) :: QCloudIceTentative
real(DP) :: xyz_Sigma (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudProdRHThreshold(0:imax-1, 1:jmax, 1:kmax)
integer :: i
integer :: j
integer :: k
integer :: l
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. cloud_T1993base_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! Save current values
!
xyz_TempB = xyz_Temp
xyz_QH2OVapB = xyz_QH2OVap
xyz_QCloudWatB = xyz_QCloudWat
xyz_QCloudIceB = xyz_QCloudIce
! Latent heat for sublimation (sum of evaporation and fusion)
LatentHeatSubl = LatentHeat + LatentHeatFusion
do k = 1, kmax
xyz_DelMass(:,:,k) = ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
end do
! 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 )
! Calculate liquid water fraction and ice fraction
!
call SaturateWatFraction( xyz_Temp, xyz_WatFrac )
xyz_IceFrac = 1.0_DP - xyz_WatFrac
xyz_QH2OVapSat = xyz_CalcQVapSat( xyz_Temp, xyz_Press )
xyz_DQH2OVapSatDTemp = xyz_CalcDQVapSatDTemp( xyz_Temp, xyz_QH2OVapSat )
!!$ xyz_QH2OVapSatOnLiq = xyz_CalcQVapSatOnLiq( xyz_Temp, xyz_Press )
!!$ xyz_DQH2OVapSatOnLiqDTemp = xyz_CalcDQVapSatOnLiqDTemp( xyz_Temp, xyz_QH2OVapSatOnLiq )
!!$ xyz_QH2OVapSatOnIce = xyz_CalcQVapSatOnIce( xyz_Temp, xyz_Press )
!!$ xyz_DQH2OVapSatOnIceDTemp = xyz_CalcDQVapSatOnIceDTemp( xyz_Temp, xyz_QH2OVapSatOnIce )
xyz_DQH2OVapSatDPressDenom = 1.0_DP + xyz_CloudCover / CpDry * ( xyz_WatFrac * LatentHeat + xyz_IceFrac * LatentHeatSubl ) * xyz_DQH2OVapSatDTemp
xyz_DQH2OVapSatDPress = ( - xyz_QH2OVapSat + GasRDry * xyz_VirTemp / CpDry * xyz_DQH2OVapSatDTemp ) / ( xyz_Press * xyz_DQH2OVapSatDPressDenom )
xyz_DQH2OVapSatDt = xyz_DQH2OVapSatDPress * ( xyz_OMG + Grav * xyz_MoistConvSubsidMassFlux ) + xyz_DQH2OVapSatDTemp / xyz_DQH2OVapSatDPressDenom * xyz_DTempDtPhy
! set zero-one flag
do k = 1, kmax
xyz_Sigma(:,:,k) = xyz_Press(:,:,k) / xyr_Press(:,:,0)
end do
xyz_CloudProdRHThreshold = RHThresholdCrtl + ( 1.0_DP - xyz_QH2OVap / xyz_QH2OVapSat ) * ( max( ( xyz_Sigma - RHThresholdSigmaMin ), 0.0_DP ) / ( xyz_Sigma - RHThresholdSigmaMin ) )**RHThresholdOrd
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_DQH2OVapSatDt(i,j,k) < 0.0_DP ) then
if ( xyz_QH2OVap(i,j,k) >= xyz_CloudProdRHThreshold(i,j,k) * xyz_QH2OVapSat(i,j,k) ) then
xyz_ZOCloudProd(i,j,k) = 1.0_DP
else
xyz_ZOCloudProd(i,j,k) = 0.0_DP
end if
xyz_ZOCloudLoss(i,j,k) = 1.0_DP
else
xyz_ZOCloudProd(i,j,k) = 0.0_DP
xyz_ZOCloudLoss(i,j,k) = 0.0_DP
end if
end do
end do
end do
! Rain and snow at the surface
xy_Rain = 0.0_DP
xy_Snow = 0.0_DP
! Rain area
xy_RainArea = 0.0_DP
k_loop : do k = kmax, 1, -1
! Freezing/melting and evaporation of precipitation
!
do j = 1, jmax
do i = 0, imax-1
call CloudUtilsPRCPStepPC1Grid( xyr_Press(i,j,k-1), xyr_Press(i,j,k), xyz_Temp(i,j,k), xy_Rain(i,j), xy_Snow(i,j) )
end do
end do
do j = 1, jmax
do i = 0, imax-1
call CloudUtilsPRCPEvap1Grid( xyz_Press(i,j,k), xyr_Press(i,j,k-1), xyr_Press(i,j,k), PRCPArea, PRCPEvapArea, xyz_Temp(i,j,k), xyz_QH2OVap(i,j,k), xy_Rain(i,j), xy_Snow(i,j) )
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_ZOCloudProd(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_QCloudWat(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_ZOCloudProd(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_QCloudWat(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_FactE1(:,:,k) * 2.0_DP * DelTime
!
!
!
!
xyz_FactAl1(:,:,k) = xyz_DQCloudWatDtCum(:,:,k)
xyz_FactAs1(:,:,k) = xyz_DQCloudIceDtCum(:,:,k)
!
!
!
xyz_FactA20(:,:,k) = - xyz_CloudCover(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * xyz_ZOCloudProd(:,:,k) - xyz_DelCloudCoverStr(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * xyz_ZOCloudProd(:,:,k) / 2.0_DP - xyz_CloudCover(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * ( 1.0_DP - xyz_ZOCloudLoss(:,:,k) ) - 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_FactA20(i,j,k) * 2.0_DP * DelTime > xyz_QH2OVap(i,j,k) ) then
xyz_FactA20(i,j,k) = xyz_QH2OVap(i,j,k) / ( 2.0_DP * DelTime ) * ( 1.0_DP - 1.0e-14_DP )
end if
end do
end do
xyz_FactAl2(:,:,k) = xyz_WatFrac(:,:,k) * xyz_FactA20(:,:,k)
xyz_FactAs2(:,:,k) = xyz_IceFrac(:,:,k) * xyz_FactA20(:,:,k)
xyz_FactAl(:,:,k) = xyz_FactAl1(:,:,k) + xyz_FactAl2(:,:,k)
xyz_FactAs(:,:,k) = xyz_FactAs1(:,:,k) + xyz_FactAs2(:,:,k)
RainConvFactor0 = 1.0_DP / ( CloudWatLifeTime0 + 1.0d-100 )
! Factor for collection
xy_FactCo = 1.0_DP + PRCPColFactor * sqrt( xy_Rain + xy_Snow )
! Factor for Bergeron-Findeisen effect
! Sundqvist et al. (1989)
!!$ xy_FactBF = 1.0_DP &
!!$ & + C2 * sqrt( max( 268.0_DP - xyz_Temp(:,:,k), 0.0_DP ) )
! An original equation following that by IFS CY38r1
! (p.82 of http://www.ecmwf.int/research/ifsdocs/CY38r1/IFSPart4.pdf)
!!$ xy_FactBF = 1.0_DP &
!!$ & + C2 * sqrt( min( &
!!$ & max( 268.0_DP - xyz_Temp(:,:,k), 0.0_DP ), &
!!$ & max( 268.0_DP - TempBFEffectSat, 0.0_DP ) &
!!$ & ) &
!!$ & )
! Constant (no effect)
xy_FactBF = 1.0_DP
!
!!$ xy_QCloudWatConvThreshold = &
!!$ & QCloudWatEffConv0 / ( xy_FactCo * xy_FactBF )
xy_QCloudWatConvThreshold = ( QCloudWatEffConv0 + 1.0e-100_DP ) / ( xy_FactCo * xy_FactBF )
!!$ xy_QCloudWatConvThreshold = &
!!$ & QCloudWatEffConv0 / xy_FactCo
xy_FactConvThreshold = 1.0_DP - exp( - ( xyz_QCloudWat(:,:,k) / ( xyz_CloudCover(:,:,k) + 1.0e-10_DP ) / xy_QCloudWatConvThreshold )**2 )
!
!!$ xy_RainConvFactor = RainConvFactor0 * xy_FactCo * xy_FactBF
xy_RainConvFactor = RainConvFactor0 * xy_FactCo * xy_FactBF * xy_FactConvThreshold
SnowConvFactor0 = 1.0_DP / ( CloudIceLifeTime0 + 1.0d-100 )
xy_FactIceTempDep = exp( 0.025_DP * ( xyz_Temp(:,:,k) - 273.15_DP ) )
!!$ xy_QCloudIceConvThreshold = QCloudIceEffConv0
xy_QCloudIceConvThreshold = QCloudIceEffConv0 + 1.0e-100_DP
xy_FactConvThreshold = 1.0_DP - exp( - ( xyz_QCloudIce(:,:,k) / ( xyz_CloudCover(:,:,k) + 1.0e-10_DP ) / xy_QCloudIceConvThreshold )**2 )
!!$ xy_SnowConvFactor = SnowConvFactor0 * xy_FactIceTempDep
xy_SnowConvFactor = SnowConvFactor0 * xy_FactIceTempDep * xy_FactConvThreshold
!
!
!
!!$ xyz_FactBl(:,:,k) = 1.0_DP / ( CloudWatLifeTime0 + 1.0d-100 )
xyz_FactBl(:,:,k) = xy_RainConvFactor
!
!!$ xyz_FactBs(:,:,k) = 1.0_DP / ( CloudIceLifeTime0 + 1.0d-100 )
xyz_FactBs(:,:,k) = xy_SnowConvFactor
!
!
!
! Values at next time step are calculated.
!
do j = 1, jmax
do i = 0, imax-1
if ( xyz_FactBl(i,j,k) >= FactBlsThreshold ) then
xyz_QCloudWat(i,j,k) = xyz_QCloudWat(i,j,k) * exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) + xyz_FactAl(i,j,k) / xyz_FactBl(i,j,k) * ( 1.0_DP - exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) )
else
xyz_QCloudWat(i,j,k) = xyz_QCloudWat(i,j,k) + xyz_FactAl(i,j,k) * ( 2.0_DP * DelTime )
end if
end do
end do
! The value of cloud water amount is checked, and xyz_FactA is updated.
do j = 1, jmax
do i = 0, imax-1
if ( xyz_QCloudWat(i,j,k) < 0.0_DP ) then
if ( xyz_FactBl(i,j,k) >= FactBlsThreshold ) then
xyz_FactAl2(i,j,k) = - xyz_FactAl1(i,j,k) - xyz_FactBl(i,j,k) * xyz_QCloudWatB(i,j,k) * exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) / ( 1.0_DP - exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) )
else
!!$ xyz_FactAl2(i,j,k) = - xyz_FactAl1(i,j,k)
xyz_FactAl2(i,j,k) = - ( xyz_QCloudWatB(i,j,k) + xyz_FactAl1(i,j,k) * ( 2.0_DP * DelTime ) ) / ( 2.0_DP * DelTime )
!!$ call MessageNotify( 'E', module_name, 'Unexpected negative QCloudWat.' )
end if
xyz_QCloudWat(i,j,k) = 0.0_DP
end if
end do
end do
!
do j = 1, jmax
do i = 0, imax-1
if ( xyz_FactBs(i,j,k) >= FactBlsThreshold ) then
xyz_QCloudIce(i,j,k) = xyz_QCloudIce(i,j,k) * exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) + xyz_FactAs(i,j,k) / xyz_FactBs(i,j,k) * ( 1.0_DP - exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) )
else
xyz_QCloudIce(i,j,k) = xyz_QCloudIce(i,j,k) + xyz_FactAs(i,j,k) * ( 2.0_DP * DelTime )
end if
end do
end do
! The value of cloud water amount is checked, and xyz_FactA is updated.
do j = 1, jmax
do i = 0, imax-1
if ( xyz_QCloudIce(i,j,k) < 0.0_DP ) then
if ( xyz_FactBs(i,j,k) >= FactBlsThreshold ) then
xyz_FactAs2(i,j,k) = - xyz_FactAs1(i,j,k) - xyz_FactBs(i,j,k) * xyz_QCloudIceB(i,j,k) * exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) / ( 1.0_DP - exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) )
else
!!$ xyz_FactAs2(i,j,k) = - xyz_FactAs1(i,j,k)
xyz_FactAs2(i,j,k) = - ( xyz_QCloudIceB(i,j,k) + xyz_FactAs1(i,j,k) * ( 2.0_DP * DelTime ) ) / ( 2.0_DP * DelTime )
!!$ call MessageNotify( 'E', module_name, 'Unexpected negative QCloudIce.' )
end if
xyz_QCloudIce(i,j,k) = 0.0_DP
end if
end do
end do
! Al and As are updated because Al2 and As2 were updated above
xyz_FactAl(:,:,k) = xyz_FactAl1(:,:,k) + xyz_FactAl2(:,:,k)
xyz_FactAs(:,:,k) = xyz_FactAs1(:,:,k) + xyz_FactAs2(:,:,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_QH2OVap(:,:,k) = xyz_QH2OVap(:,:,k) - ( xyz_FactAl2(:,:,k) + xyz_FactAs2(:,:,k) ) * 2.0_DP * DelTime
!
!!$ xyz_Temp(:,:,k) = xyz_Temp(:,:,k) &
!!$ & + xyz_FactA2(:,:,k) * 2.0_DP * DelTime * LatentHeat / CpDry
xyz_Temp(:,:,k) = xyz_Temp(:,:,k) + ( xyz_FactAl2(:,:,k) * LatentHeat + xyz_FactAs2(:,:,k) * LatentHeatSubl ) * 2.0_DP * DelTime / CpDry
! Rain
!
do j = 1, jmax
do i = 0, imax-1
!!$ if ( xyz_FactBl(i,j,k) >= FactBlsThreshold ) then
!!$ xyz_Rain(i,j,k) = &
!!$ & xyz_FactAl(i,j,k) * 2.0_DP * DelTime &
!!$ & + xyz_QCloudWatB(i,j,k) &
!!$ & * ( 1.0_DP - exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) ) &
!!$ & - xyz_FactAl(i,j,k) / xyz_FactBl(i,j,k) &
!!$ & * ( 1.0_DP - exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) )
!!$ xyz_Rain(i,j,k) = xyz_Rain(i,j,k) / ( 2.0_DP * DelTime )
!!$ else
!!$ xyz_Rain(i,j,k) = 0.0_DP
!!$ end if
xyz_Rain(i,j,k) = xyz_QCloudWatB(i,j,k) + xyz_FactAl(i,j,k) * 2.0_DP * DelTime - xyz_QCloudWat (i,j,k)
xyz_Rain(i,j,k) = xyz_Rain(i,j,k) / ( 2.0_DP * DelTime )
end do
end do
! Snow
!
do j = 1, jmax
do i = 0, imax-1
!!$ if ( xyz_FactBs(i,j,k) >= FactBlsThreshold ) then
!!$ xyz_Snow(i,j,k) = &
!!$ & xyz_FactAs(i,j,k) * 2.0_DP * DelTime &
!!$ & + xyz_QCloudIceB (i,j,k) &
!!$ & * ( 1.0_DP - exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) ) &
!!$ & - xyz_FactAs(i,j,k) / xyz_FactBs(i,j,k) &
!!$ & * ( 1.0_DP - exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) )
!!$ xyz_Snow(i,j,k) = xyz_Snow(i,j,k) / ( 2.0_DP * DelTime )
!!$ else
!!$ xyz_Snow(i,j,k) = 0.0_DP
!!$ end if
xyz_Snow(i,j,k) = xyz_QCloudIceB (i,j,k) + xyz_FactAs(i,j,k) * 2.0_DP * DelTime - xyz_QCloudIce (i,j,k)
xyz_Snow(i,j,k) = xyz_Snow(i,j,k) / ( 2.0_DP * DelTime )
end do
end do
! Rain at the surface
xy_Rain = xy_Rain + xyz_Rain(:,:,k) * xyz_DelMass(:,:,k)
! Snow at the surface
xy_Snow = xy_Snow + xyz_Snow(:,:,k) * xyz_DelMass(:,:,k)
! Treatment of supersaturation
do j = 1, jmax
do i = 0, imax-1
QCloudWatTentative = xyz_QCloudWat(i,j,k)
QCloudIceTentative = xyz_QCloudIce(i,j,k)
call LScaleCond1Grid( xyz_Temp(i,j,k), xyz_QH2OVap(i,j,k), QCloudWatTentative, QCloudIceTentative, xyz_Press(i,j,k), xyr_Press(i,j,k-1), xyr_Press(i,j,k), xyz_DQH2OLiqDtLSC(i,j,k), xyz_DQH2OSolDtLSC(i,j,k) )
!!$ xy_Rain(i,j) = xy_Rain(i,j) &
!!$ & + ( xyz_QCloudWat(i,j,k) - QCloudWatTentative ) &
!!$ & / ( 2.0_DP * DelTime ) &
!!$ & * ( 2.0_DP * DelTime ) &
!!$ & * xyz_DelMass(i,j,k)
!!$ xy_Snow(i,j) = xy_Snow(i,j) &
!!$ & + ( xyz_QCloudIce(i,j,k) - QCloudIceTentative ) &
!!$ & / ( 2.0_DP * DelTime ) &
!!$ & * ( 2.0_DP * DelTime ) &
!!$ & * xyz_DelMass(i,j,k)
!!$ xyz_QCloudWat(i,j,k) = xyz_QCloudWat(i,j,k) &
!!$ & - ( xyz_QCloudWat(i,j,k) - QCloudWatTentative ) &
!!$ & / ( 2.0_DP * DelTime ) &
!!$ & * ( 2.0_DP * DelTime )
!!$ xyz_QCloudIce(i,j,k) = xyz_QCloudIce(i,j,k) &
!!$ & - ( xyz_QCloudIce(i,j,k) - QCloudIceTentative ) &
!!$ & / ( 2.0_DP * DelTime ) &
!!$ & * ( 2.0_DP * DelTime )
end do
end do
xy_Rain = xy_Rain + xyz_DQH2OLiqDtLSC(:,:,k) * xyz_DelMass(:,:,k)
xy_Snow = xy_Snow + xyz_DQH2OSolDtLSC(:,:,k) * xyz_DelMass(:,:,k)
! Evaporation
!
do j = 1, jmax
do i = 0, imax-1
if ( ( xyz_QCloudWat(i,j,k) < QCloudWaterEvapThreshold ) .and. ( xyz_QCloudIce(i,j,k) < QCloudWaterEvapThreshold ) ) then
xyz_CloudCover(i,j,k) = 0.0_DP
end if
end do
end do
! Cloud cover is restricted.
xyz_CloudCover(:,:,k) = min( max( xyz_CloudCover(:,:,k), 0.0_DP ), 1.0_DP )
! Check values
do j = 1, jmax
do i = 0, imax-1
if ( xyz_QH2OVap(i,j,k) < 0.0_DP ) then
write( 6, * ) 'QH2OVap is negative', i, j, k, xyz_QH2OVap(i,j,k)
end if
if ( xyz_QCloudWat(i,j,k) < 0.0_DP ) then
write( 6, * ) 'QCloudWat is negative', i, j, k, xyz_QCloudWat(i,j,k)
end if
if ( xyz_QCloudIce (i,j,k) < 0.0_DP ) then
write( 6, * ) 'QCloudIce is negative', i, j, k, xyz_QCloudIce (i,j,k)
end if
end do
end do
end do k_loop
! 大規模凝結 (非対流性凝結) (Manabe, 1965)
! Large scale condensation (non-convective condensation)
! (Manabe, 1965)
!
! It should be noted that H2OLiq and H2OSol have updated in above
! subroutine.
!
! This routine is commented out, since this is done in an above k_loop.
!!$ call LScaleCond1D3DWrapper( &
!!$ & xyz_Temp, xyz_QH2OVap, & ! (inout)
!!$ & xyz_QCloudWat, & ! (inout)
!!$ & xyz_QCloudIce, & ! (inout)
!!$ & xyz_Press, xyr_Press, & ! (in)
!!$ & xyz_DQH2OLiqDtLSC, xyz_DQH2OSolDtLSC & ! (out)
!!$ & )
xy_SurfRainFlux = xy_Rain
xy_SurfSnowFlux = xy_Snow
call CloudUtilConsChk( "CloudT1993baseWithIce", xyr_Press, xyz_TempB, xyz_QH2OVapB, xyz_QCloudWatB, xyz_QCloudIceB, xyz_Temp , xyz_QH2OVap , xyz_QCloudWat , xyz_QCloudIce , xy_SurfRainFlux, xy_SurfSnowFlux )
end subroutine CloudT1993baseWithIce
| Variable : | |||
| TempBFEffectSat : | real(DP), save
|
| Variable : | |||
| cloud_T1993base_inited = .false. : | logical, save
|
| Constant : | |||
| module_name = ‘cloud_T1993base‘ : | character(*), parameter
|
| Constant : | |||
| version = ’$Name: $’ // ’$Id: cloud_T1993base.f90,v 1.6 2015/01/29 12:20:40 yot Exp $’ : | character(*), parameter
|