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
|