MODULE lmdz_lscp_condensation !---------------------------------------------------------------- ! Module for condensation of clouds routines ! that are called in LSCP IMPLICIT NONE CONTAINS !********************************************************************************** SUBROUTINE condensation_lognormal( & klon, temp, qtot, qsat, gamma_cond, ratqs, & keepgoing, cldfra, qincld, qvc) !---------------------------------------------------------------------- ! This subroutine calculates the formation of clouds, using a ! statistical cloud scheme. It uses a generalised lognormal distribution ! of total water in the gridbox ! See Bony and Emanuel, 2001 !---------------------------------------------------------------------- USE lmdz_lscp_ini, ONLY: eps IMPLICIT NONE ! ! Input ! INTEGER, INTENT(IN) :: klon ! number of horizontal grid points ! REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K] REAL, INTENT(IN) , DIMENSION(klon) :: qtot ! total specific humidity (without precip) [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-] REAL, INTENT(IN) , DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-] LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed ! ! Output ! NB. those are in INOUT because of the convergence loop on temperature ! (in some cases, the values are not re-computed) but the values ! are never used explicitely ! REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra ! cloud fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: qincld ! cloud-mean in-cloud total specific water [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: qvc ! gridbox-mean vapor in the cloud [kg/kg] ! ! Local ! INTEGER :: i REAL :: pdf_std, pdf_k, pdf_delta REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2 ! !--Loop on klon ! DO i = 1, klon !--If a new calculation of the condensation is needed, !--i.e., temperature has not yet converged (or the cloud is !--formed elsewhere) IF (keepgoing(i)) THEN pdf_std = ratqs(i) * qtot(i) pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) ) pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) ) pdf_a = pdf_delta / ( pdf_k * SQRT(2.) ) pdf_b = pdf_k / (2. * SQRT(2.)) pdf_e1 = pdf_a - pdf_b pdf_e1 = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 ) pdf_e1 = 1. - ERF(pdf_e1) pdf_e2 = pdf_a + pdf_b pdf_e2 = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 ) pdf_e2 = 1. - ERF(pdf_e2) IF ( pdf_e1 .LT. eps ) THEN cldfra(i) = 0. qincld(i) = qsat(i) !--AB grid-mean vapor in the cloud - we assume saturation adjustment qvc(i) = 0. ELSE cldfra(i) = 0.5 * pdf_e1 qincld(i) = qtot(i) * pdf_e2 / pdf_e1 !--AB grid-mean vapor in the cloud - we assume saturation adjustment qvc(i) = qsat(i) * cldfra(i) ENDIF ENDIF ! end keepgoing ENDDO ! end loop on i END SUBROUTINE condensation_lognormal !********************************************************************************** !********************************************************************************** SUBROUTINE condensation_ice_supersat( & klon, dtime, pplay, paprsdn, paprsup, totfra_in, & cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, & temp, qtot_in, qsat, gamma_cond, ratqs, keepgoing, pt_pron_clds, & dzsed_abv, flsed_abv, cfsed_abv, dzsed, flsed, cfsed, flauto, & cldfra, qincld, qvc, issrfra, qissr, & dcf_sub, dcf_con, dcf_mix, dcf_sed, dcf_auto, & dqi_adj, dqi_sub, dqi_con, dqi_mix, dqi_sed, dqi_auto, & dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, dqvc_sed, dqvc_auto, & contfra_in, qvcont_in, qicont_in, nic_in, flight_dist, flight_fuel, & contfra, qcont, qvcont, Ncont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, & AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails, & dzsed_cont_abv, flsed_cont_abv, Nflsed_cont_abv, cfsed_cont_abv, & dzsed_cont, flsed_cont, Nflsed_cont, cfsed_cont, & dcfc_ini, dqic_ini, dqtc_ini, dnic_ini, dcfc_sub, dqic_sub, dqtc_sub, dnic_sub, & dcfc_mix, dqic_mix, dqtc_mix, dnic_mix, dnic_agg, dqic_adj, dqtc_adj, & dcfc_sed, dqic_sed, dqtc_sed, dnic_sed, dcfc_auto, dqic_auto, dqtc_auto, dnic_auto) !---------------------------------------------------------------------- ! This subroutine calculates the formation, evolution and dissipation ! of clouds, using a process-oriented treatment of the cloud properties ! (cloud fraction, vapor in the cloud, condensed water in the cloud). ! It allows for ice supersaturation in cold regions, in clear sky. ! If ok_unadjusted_clouds, it also allows for sub- and supersaturated ! cloud water vapors. ! It also allows for the formation and evolution of condensation trails ! (contrails) from aviation. ! Authors: Audran Borella, Etienne Vignon, Olivier Boucher ! April 2024 !---------------------------------------------------------------------- USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC USE lmdz_lscp_ini, ONLY: RLSTT, RTT, RD, RG, RV, RPI, EPS_W USE lmdz_lscp_ini, ONLY: eps, temp_nowater, ok_unadjusted_clouds USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice USE lmdz_lscp_ini, ONLY: N_ice_volume, corr_incld_depsub, nu_iwc_pdf_lscp USE lmdz_lscp_ini, ONLY: beta_pdf_lscp, temp_thresh_pdf_lscp USE lmdz_lscp_ini, ONLY: std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp USE lmdz_lscp_ini, ONLY: coef_mixing_lscp, coef_shear_lscp, chi_mixing USE lmdz_lscp_ini, ONLY: aspect_ratio_cirrus, cooling_rate_ice_thresh USE lmdz_lscp_ini, ONLY: ok_ice_sedim, ffallv_sed, cice_velo, dice_velo USE lmdz_lscp_ini, ONLY: chi_sedim USE lmdz_lscp_ini, ONLY: ok_plane_contrail, ok_unadjusted_contrails USE lmdz_lscp_ini, ONLY: aspect_ratio_contrails, coef_mixing_contrails, coef_shear_contrails USE lmdz_lscp_ini, ONLY: chi_mixing_contrails, eff2vol_radius_contrails USE lmdz_lscp_tools, ONLY: ICECRYSTAL_VELO USE lmdz_aviation, ONLY: contrails_formation IMPLICIT NONE ! ! Input ! INTEGER, INTENT(IN) :: klon ! number of horizontal grid points REAL, INTENT(IN) :: dtime ! time step [s] ! REAL, INTENT(IN) , DIMENSION(klon) :: pplay ! layer pressure [Pa] REAL, INTENT(IN) , DIMENSION(klon) :: paprsdn ! pressure at the lower interface [Pa] REAL, INTENT(IN) , DIMENSION(klon) :: paprsup ! pressure at the upper interface [Pa] REAL, INTENT(IN) , DIMENSION(klon) :: totfra_in ! total available fraction for stratiform clouds [-] REAL, INTENT(IN) , DIMENSION(klon) :: cldfra_in ! cloud fraction [-] REAL, INTENT(IN) , DIMENSION(klon) :: qvc_in ! gridbox-mean water vapor in cloud [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qliq_in ! specific liquid water content [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qice_in ! specific ice water content [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: shear ! vertical shear [s-1] REAL, INTENT(IN) , DIMENSION(klon) :: pbl_eps ! TKE dissipation [m2/s3] REAL, INTENT(IN) , DIMENSION(klon) :: cell_area ! cell area [m2] REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K] REAL, INTENT(IN) , DIMENSION(klon) :: qtot_in ! total specific humidity (without precip) [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-] REAL, INTENT(IN) , DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-] LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed LOGICAL, INTENT(IN) , DIMENSION(klon) :: pt_pron_clds ! .TRUE. if clouds are prognostic in this mesh REAL, INTENT(IN) , DIMENSION(klon) :: dzsed_abv ! sedimentated cloud height IN THE LAYER ABOVE [m] REAL, INTENT(IN) , DIMENSION(klon) :: flsed_abv ! sedimentated ice flux FROM THE LAYER ABOVE [kg/s/m2] REAL, INTENT(IN) , DIMENSION(klon) :: cfsed_abv ! cloud fraction IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: dzsed ! sedimentated cloud height [m] REAL, INTENT(INOUT), DIMENSION(klon) :: flsed ! sedimentated ice flux [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: cfsed ! sedimentated cloud fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: flauto ! autoconverted ice flux [kg/s/m2] ! ! Input for aviation ! REAL, INTENT(IN) , DIMENSION(klon) :: contfra_in ! input contrails fraction [-] REAL, INTENT(IN) , DIMENSION(klon) :: qvcont_in ! input contrails specific humidity [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: qicont_in ! input contrails ice specific humidity [kg/kg] REAL, INTENT(IN) , DIMENSION(klon) :: nic_in ! input contrails ice crystals concentration [#/kg] REAL, INTENT(IN) , DIMENSION(klon) :: flight_dist ! aviation distance flown concentration [m/s/m3] REAL, INTENT(IN) , DIMENSION(klon) :: flight_fuel ! aviation fuel consumption concentration [kg/s/m3] REAL, INTENT(IN) , DIMENSION(klon) :: dzsed_cont_abv! sedimentated contrails height IN THE LAYER ABOVE [m] REAL, INTENT(IN) , DIMENSION(klon) :: flsed_cont_abv! sedimentated ice flux in contrails FROM THE LAYER ABOVE [kg/s/m2] REAL, INTENT(IN) , DIMENSION(klon) :: Nflsed_cont_abv! sedimentated ice crystals concentration in contrails FROM THE LAYER ABOVE [#/s/m2] REAL, INTENT(IN) , DIMENSION(klon) :: cfsed_cont_abv! contrails fraction IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: dzsed_cont ! sedimentated contrails height [m] REAL, INTENT(INOUT), DIMENSION(klon) :: flsed_cont ! sedimentated ice flux in contrails [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: Nflsed_cont ! sedimentated ice crystals concentration flux in contrails [#/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: cfsed_cont ! sedimentated contrails fraction [-] ! ! Output ! NB. cldfra and qincld should be outputed as cf_seri and qi_seri, ! or as tendencies (maybe in the future) ! NB. those are in INOUT because of the convergence loop on temperature ! (in some cases, the values are not re-computed) but the values ! are never used explicitely ! REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra ! cloud fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: qincld ! cloud-mean in-cloud total specific water [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: qvc ! gridbox-mean vapor in the cloud [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: issrfra ! ISSR fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: qissr ! gridbox-mean ISSR specific water [kg/kg] ! ! Diagnostics for condensation and ice supersaturation ! NB. idem for the INOUT ! REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_sub ! cloud fraction tendency because of sublimation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_con ! cloud fraction tendency because of condensation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_mix ! cloud fraction tendency because of cloud mixing [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_sed ! cloud fraction tendency because of sedimentation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_auto ! cloud fraction tendency because of autoconversion [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_adj ! specific ice content tendency because of temperature adjustment [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_sub ! specific ice content tendency because of sublimation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_con ! specific ice content tendency because of condensation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_mix ! specific ice content tendency because of cloud mixing [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_sed ! specific ice content tendency because of sedimentation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_auto ! specific ice content tendency because of autoconversion [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_sed ! specific cloud water vapor tendency because of sedimentation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_auto! specific cloud water vapor tendency because of autoconversion [kg/kg/s] ! ! Diagnostics for aviation ! NB. idem for the INOUT ! REAL, INTENT(INOUT), DIMENSION(klon) :: contfra ! contrail fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: qcont ! contrail total specific humidity [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: qvcont ! contrail specific humidity [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: Ncont ! contrail ice crystals concentration [#/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: Tcritcont ! critical temperature for contrail formation [K] REAL, INTENT(INOUT), DIMENSION(klon) :: qcritcont ! critical specific humidity for contrail formation [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: potcontfraP ! potential persistent contrail fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: AEI_contrails ! Apparent emission index contrails [#/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: AEI_surv_contrails ! Apparent emission index contrails after vortex loss [#/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: fsurv_contrails ! Survival fraction after vortex loss [-] REAL, INTENT(INOUT), DIMENSION(klon) :: section_contrails ! Cross section of newly formed contrails [m2] REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_ini ! contrails cloud fraction tendency because of initial formation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_ini ! contrails ice specific humidity tendency because of initial formation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_ini ! contrails total specific humidity tendency because of initial formation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_ini ! contrails ice crystals concentration tendency because of initial formation [#/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_sub ! contrails cloud fraction tendency because of sublimation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_sub ! contrails ice specific humidity tendency because of sublimation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_sub ! contrails total specific humidity tendency because of sublimation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_sub ! contrails ice crystals concentration tendency because of sublimation [#/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_adj ! contrails ice specific humidity tendency because of phase adjustment [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_adj ! contrails total specific humidity tendency because of phase adjustment [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_mix ! contrails cloud fraction tendency because of mixing [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_mix ! contrails ice specific humidity tendency because of mixing [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_mix ! contrails total specific humidity tendency because of mixing [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_mix ! contrails ice crystals concentration tendency because of mixing [#/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_agg ! contrails ice crystals concentration tendency because of aggregation [#/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_sed ! contrails cloud fraction tendency because of ice sedimentation [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_sed ! contrails ice specific humidity tendency because of ice sedimentation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_sed ! contrails total specific humidity tendency because of ice sedimentation [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_sed ! contrails ice crystals concentration tendency because of ice sedimentation [#/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_auto ! contrails cloud fraction tendency because of ice autoconversion [s-1] REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_auto ! contrails ice specific humidity tendency because of ice autoconversion [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_auto ! contrails total specific humidity tendency because of ice autoconversion [kg/kg/s] REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_auto ! contrails ice crystals concentration tendency because of ice autoconversion [#/kg/s] ! ! Local ! INTEGER :: i LOGICAL :: ok_warm_cloud REAL, DIMENSION(klon) :: qcld, qzero REAL, DIMENSION(klon) :: clrfra, qclr ! ! for lognormal REAL :: pdf_std, pdf_k, pdf_delta REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2 ! ! for unadjusted clouds REAL :: qiceincld, qvapincld, qvapincld_new REAL :: qice_ratio ! ! for deposition / sublimation REAL :: pres_sat, kappa_depsub, tauinv_depsub REAL :: air_thermal_conduct, water_vapor_diff ! ! for dissipation REAL, DIMENSION(klon) :: temp_diss, qsati_diss, qiceincld_min REAL :: pdf_shape ! ! for condensation REAL, DIMENSION(klon) :: qsatl, dqsat_tmp REAL, DIMENSION(klon) :: pdf_alpha, pdf_scale, pdf_gamma REAL :: rhl_clr, pdf_loc REAL :: pdf_e3, pdf_x, pdf_y REAL :: dqt_con ! ! for sedimentation and autoconversion REAL, DIMENSION(klon) :: qised_abv REAL :: iwc, icefall_velo, coef_sed, qice_sedim, Nice_sedim REAL :: sedfra_abv, sedfra1, sedfra2, sedfra3, sedfra_tmp REAL :: dcf_sed1, dcf_sed2, dqvc_sed1, dqvc_sed2, dqi_sed1, dqi_sed2 REAL :: dz_auto, coef_auto, qice_auto ! ! for aggregation REAL :: dei, sticking_coef REAL :: gamma_agg, dispersion_coef ! ! for mixing REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix REAL :: cldfra_mix, clrfra_mix, sigma_mix REAL :: L_shear, shear_fra REAL :: qvapinmix, qiceinmix, qvapinmix_lim, qvapinclr_lim REAL :: pdf_fra_above_nuc, pdf_q_above_nuc REAL :: pdf_fra_above_lim, pdf_q_above_lim REAL :: pdf_fra_below_lim REAL :: mixed_fraction ! ! for cell properties REAL, DIMENSION(klon) :: rho REAL :: rhodz, dz ! ! for contrails REAL :: mice, Niceincld, qvapincont REAL, DIMENSION(klon) :: qised_cont_abv, Nised_cont_abv REAL :: dcfc_sed1, dcfc_sed2, dqtc_sed1, dnic_sed1, dqtc_sed2, dqic_sed1, dqic_sed2, dnic_sed2 REAL :: dz_auto_cont qzero(:) = 0. !--Calculation of qsat w.r.t. liquid CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsat_tmp) !--Calculation of qsat max for dissipation temp_diss(:) = temp(:) + cooling_rate_ice_thresh * dtime CALL calc_qsat_ecmwf(klon, temp_diss, qzero, pplay, RTT, 2, .FALSE., qsati_diss, dqsat_tmp) !--Additionally to a minimum in cloud water vapor, we impose a minimum !--on the in-cloud ice water content. It is calculated following !--Marti and Mauersberger (1993), see also Schiller et al. (2008) qiceincld_min(:) = qsati_diss(:) - qsat(:) ! !--Loop on klon ! DO i = 1, klon !--If a new calculation of the condensation is needed, !--i.e., temperature has not yet converged (or the cloud is !--formed elsewhere) IF (keepgoing(i)) THEN !--Initialisation issrfra(i) = 0. qissr(i) = 0. IF ( ok_ice_sedim ) THEN flauto(i) = 0. flsed(i) = 0. dzsed(i) = 0. cfsed(i) = 0. ENDIF IF ( ok_plane_contrail ) THEN contfra(i) = 0. qcont(i) = 0. qvcont(i) = 0. Ncont(i) = 0. IF ( ok_ice_sedim ) THEN flsed_cont(i) = 0. Nflsed_cont(i) = 0. dzsed_cont(i) = 0. cfsed_cont(i) = 0. ENDIF ENDIF !--If the temperature is higher than the threshold below which !--there is no liquid in the gridbox, we activate the usual scheme !--(generalised lognormal from Bony and Emanuel 2001) !--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for !--all clouds, and the lognormal scheme is not activated IF ( .NOT. pt_pron_clds(i) ) THEN pdf_std = ratqs(i) * qtot_in(i) pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot_in(i))**2 ) ) pdf_delta = LOG( qtot_in(i) / ( gamma_cond(i) * qsat(i) ) ) pdf_a = pdf_delta / ( pdf_k * SQRT(2.) ) pdf_b = pdf_k / (2. * SQRT(2.)) pdf_e1 = pdf_a - pdf_b pdf_e1 = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 ) pdf_e1 = 1. - ERF(pdf_e1) pdf_e2 = pdf_a + pdf_b pdf_e2 = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 ) pdf_e2 = 1. - ERF(pdf_e2) IF ( pdf_e1 .LT. eps ) THEN cldfra(i) = 0. qincld(i) = qsat(i) qvc(i) = 0. ELSE cldfra(i) = 0.5 * pdf_e1 qincld(i) = qtot_in(i) * pdf_e2 / pdf_e1 qvc(i) = qsat(i) * cldfra(i) ENDIF !--If the temperature is lower than temp_nowater, we use the new !--condensation scheme that allows for ice supersaturation ELSE !--Initialisation !--If the air mass is warm (liquid water can exist), !--all the memory is lost and the scheme becomes statistical, !--i.e., the sublimation and mixing processes are deactivated, !--and the condensation process is slightly adapted !--This can happen only if ok_weibull_warm_clouds = .TRUE. ok_warm_cloud = ( temp(i) .GT. temp_nowater ) !--The following barriers ensure that the traced cloud properties !--are consistent. In some rare cases, i.e. the cloud water vapor !--can be greater than the total water in the gridbox cldfra(i) = MAX(0., MIN(totfra_in(i), cldfra_in(i))) qcld(i) = MAX(0., MIN(qtot_in(i), qliq_in(i) + qice_in(i) + qvc_in(i))) qvc(i) = MAX(0., MIN(qcld(i), qvc_in(i))) !--Initialise clear fraction properties clrfra(i) = totfra_in(i) - cldfra(i) qclr(i) = qtot_in(i) - qcld(i) dcf_sub(i) = 0. dqi_sub(i) = 0. dqvc_sub(i) = 0. dqi_adj(i) = 0. dqvc_adj(i) = 0. dcf_con(i) = 0. dqi_con(i) = 0. dqvc_con(i) = 0. dcf_mix(i) = 0. dqi_mix(i) = 0. dqvc_mix(i) = 0. dcf_sed(i) = 0. dqi_sed(i) = 0. dqvc_sed(i) = 0. qised_abv(i) = 0. IF ( dzsed_abv(i) .GT. eps ) THEN !--If ice sedimentation is activated, the quantity of sedimentated ice was added !--to the total water vapor in the precipitation routine. Here we remove it !--(it will be reincluded later) !--The barrier is needed because although the sedimented ice was vaporised, !--inacurracies in the advection scheme might lead to considering this water !--to be cloudy (hence not in the clear sky region) qised_abv(i) = MIN(qclr(i), flsed_abv(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime) qclr(i) = qclr(i) - qised_abv(i) ENDIF !--Initialisation of the cell properties !--Dry density [kg/m3] rho(i) = pplay(i) / temp(i) / RD !--Dry air mass [kg/m2] rhodz = ( paprsdn(i) - paprsup(i) ) / RG !--Cell thickness [m] dz = rhodz / rho(i) !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment !--hypothesis is lost, and the vapor in the cloud is purely prognostic. ! !--The deposition equation is !-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) !--from Lohmann et al. (2016), where !--alpha is the deposition coefficient [-] !--mi is the mass of one ice crystal [kg] !--C is the capacitance of an ice crystal [m] !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-] !--R_v is the specific gas constant for humid air [J/kg/K] !--T is the temperature [K] !--esi is the saturation pressure w.r.t. ice [Pa] !--Dv is the diffusivity of water vapor [m2/s] !--Ls is the specific latent heat of sublimation [J/kg/K] !--ka is the thermal conductivity of dry air [J/m/s/K] ! !--alpha is a coefficient to take into account the fact that during deposition, a water !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays !--coherent (with the same structure). It has no impact for sublimation. !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016)) !--during deposition, and alpha = 1. during sublimation. !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus !-- C = capa_cond_cirrus * rm_ice ! !--We have qice = Nice * mi, where Nice is the ice crystal !--number concentration per kg of moist air !--HYPOTHESIS 1: the ice crystals are spherical, therefore !-- mi = 4/3 * pi * rm_ice**3 * rho_ice !--HYPOTHESIS 2: the ice crystals concentration is constant in the cloud ! !--The equation in terms of q_ice is valide locally, and the local ice water content !--follows a Gamma distribution with a factor nu_iwc_pdf_lscp. Therefore, by !--integrating the local equation over the PDF (entire cloud), a correcting factor !--must be included, equal to !-- corr_incld_depsub = GAMMA(nu + 1/3) / GAMMA(nu) / nu**(1/3) !--NB. this is equal to about 0.9, hence the correction is not big !--NB. to lighten the calculated, corr_incld_depsub is calculated in lmdz_lscp_ini ! !--As the deposition process does not create new ice crystals, !--and because we assume a same rm_ice value for all crystals !--therefore the sublimation process does not destroy ice crystals !--(or, in a limit case, it destroys all ice crystals), then !--Nice is a constant during the sublimation/deposition process !--hence dmi = dqi ! !--The deposition equation then reads for qi the in-cloud ice water content: !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat * corr_incld_depsub & !-- / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice !-- dqi/dt = alpha*4pi*capa_cond_cirrus*Nice*corr_incld_depsub & !-- / ( 4pi/3 N_ice rho_ice )**(1/3) & !-- / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) & !-- qi**(1/3) * (qvc - qsat) / qsat !--and we have !-- dqvc/dt = - alpha * kappa(T) * qi**(1/3) * (qvc - qsat) !-- dqi/dt = alpha * kappa(T) * qi**(1/3) * (qvc - qsat) ! !--This system of equations can be resolved with an exact !--explicit numerical integration, having one variable resolved !--explicitly, the other exactly. qvc is always the variable solved exactly. ! !--kappa is computed as an initialisation constant, as it depends only !--on temperature and other pre-computed values pres_sat = qsat(i) / ( EPS_W + ( 1. - EPS_W ) * qsat(i) ) * pplay(i) !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971) air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184 !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976) water_vapor_diff = 0.211 * ( temp(i) / RTT )**1.94 * ( 101325. / pplay(i) ) * 1.e-4 !--NB. the greater kappa_depsub, the more efficient is the !--deposition/sublimation process kappa_depsub = 4. * RPI * capa_cond_cirrus * N_ice_volume / rho(i) * corr_incld_depsub & / qsat(i) / ( 4. / 3. * RPI * N_ice_volume / rho(i) * rho_ice )**(1./3.) & / ( RV * temp(i) / water_vapor_diff / pres_sat & + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) ) !--If contrails are activated IF ( ok_plane_contrail ) THEN !--The following barriers are needed since the advection scheme does not !--conserve order relations contfra(i) = MAX(0., MIN(cldfra(i), contfra_in(i))) Ncont(i) = MAX(0., nic_in(i)) IF ( ok_unadjusted_contrails ) THEN qcont(i) = MAX(0., MIN(qcld(i), qvcont_in(i) + qicont_in(i))) qvcont(i) = MAX(0., MIN(qcont(i), qvcont_in(i))) ELSE qcont(i) = MAX(0., MIN(qcld(i), qvcont_in(i))) qvcont(i) = qsat(i) * contfra(i) ENDIF qised_cont_abv(i) = 0. Nised_cont_abv(i) = 0. IF ( dzsed_cont_abv(i) .GT. eps ) THEN !--If ice sedimentation is activated, the quantity of sedimentated ice was added !--to the total water vapor in the precipitation routine. Here we remove it !--(it will be reincluded later) !--The barrier is needed because although the sedimented ice was vaporised, !--inacurracies in the advection scheme might lead to considering this water !--to be cloudy (hence not in the clear sky region) qised_cont_abv(i) = MIN(qclr(i), flsed_cont_abv(i) & / ( paprsdn(i) - paprsup(i) ) * RG * dtime) Nised_cont_abv(i) = Nflsed_cont_abv(i) & / ( paprsdn(i) - paprsup(i) ) * RG * dtime qclr(i) = qclr(i) - qised_cont_abv(i) ENDIF dcfc_ini(i) = 0. dqic_ini(i) = 0. dqtc_ini(i) = 0. dnic_ini(i) = 0. dcfc_sub(i) = 0. dqic_sub(i) = 0. dqtc_sub(i) = 0. dnic_sub(i) = 0. dqic_adj(i) = 0. dqtc_adj(i) = 0. dcfc_mix(i) = 0. dqic_mix(i) = 0. dqtc_mix(i) = 0. dnic_mix(i) = 0. dnic_agg(i) = 0. dcfc_sed(i) = 0. dqic_sed(i) = 0. dqtc_sed(i) = 0. dnic_sed(i) = 0. dcfc_auto(i) = 0. dqic_auto(i) = 0. dqtc_auto(i) = 0. dnic_auto(i) = 0. ELSE contfra(i) = 0. qcont(i) = 0. qvcont(i) = 0. Ncont(i) = 0. ENDIF !---------------------------------------------------------------------- !-- SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CONTRAIL -- !---------------------------------------------------------------------- !--If there is a contrail IF ( contfra(i) .GT. eps ) THEN !--We remove contrails from the main class cldfra(i) = MAX(0., cldfra(i) - contfra(i)) qcld(i) = MAX(0., qcld(i) - qcont(i)) qvc(i) = MAX(0., MIN(qcld(i), qvc(i) - qvcont(i))) qvapincld = qvcont(i) / contfra(i) qiceincld = ( qcont(i) / contfra(i) - qvapincld ) !--The in-cloud ice crystals concentration is conserved Niceincld = Ncont(i) / contfra(i) !--If the ice water content is too low, the cloud is purely sublimated !--Most probably, we advected a cloud with no ice water content (possible !--if the entire cloud precipited for example) !--Same if the number concentration is too low (less than 1 #/m3) IF ( ( qiceincld .LT. eps ) .OR. ( ( Niceincld * rho(i) ) .LT. 1. ) ) THEN dcfc_sub(i) = - contfra(i) dqtc_sub(i) = - qvcont(i) dqic_sub(i) = - ( qcont(i) - qvcont(i) ) dnic_sub(i) = - Ncont(i) contfra(i) = 0. qcont(i) = 0. qvcont(i) = 0. Ncont(i) = 0. clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcfc_sub(i)) qclr(i) = qclr(i) - dqtc_sub(i) - dqic_sub(i) !--Else, the cloud is adjusted and sublimated ELSE IF ( ok_unadjusted_contrails ) THEN IF ( qvapincld .GE. qsat(i) ) THEN !--If the cloud is initially supersaturated !--Exact explicit integration (qvcont exact, qice explicit) tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub & * (Niceincld / N_ice_volume)**(2./3.) qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub ) ELSE !--If the cloud is initially subsaturated !--Exact explicit integration (qice exact, qvcont explicit) !--The barrier is set so that the resulting vapor in cloud !--cannot be greater than qsat !--qice_ratio is the ratio between the new ice content and !--the old one, it is comprised between 0 and 1 tauinv_depsub = qiceincld**(1./3.) * kappa_depsub & * (Niceincld / N_ice_volume)**(2./3.) qice_ratio = tauinv_depsub * dtime / 1.5 / qiceincld * ( qsat(i) - qvapincld ) !--The new vapor in the cloud is increased with the !--sublimated ice qvapincld_new = qvapincld + qiceincld * ( 1. - MAX(0., 1. - qice_ratio)**1.5 ) !--The new vapor in the cloud cannot be greater than qsat qvapincld_new = MIN(qvapincld_new, qsat(i)) !--If all the ice is sublimated IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0. ENDIF ! qvapincld .GT. qsat ELSE !--We keep the saturation adjustment hypothesis, and the vapor in the !--cloud is set equal to the saturation vapor IF ( ( qvapincld + qiceincld ) .GT. qsat(i) ) THEN qvapincld_new = qsat(i) ELSE qvapincld_new = 0. ENDIF ENDIF ! ok_unadjusted_contrails !--------------------------------------- !-- DISSIPATION OF THE CONTRAIL -- !--------------------------------------- !--If the dissipation process must be activated IF ( ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) ) .GT. qvapincld ) THEN !--Gamma distribution starting at qvapincld pdf_shape = nu_iwc_pdf_lscp / qiceincld pdf_y = pdf_shape * ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) - qvapincld ) pdf_e1 = GAMMAINC ( nu_iwc_pdf_lscp , pdf_y ) pdf_e2 = GAMMAINC ( nu_iwc_pdf_lscp + 1. , pdf_y ) !--Tendencies and diagnostics dcfc_sub(i) = - contfra(i) * pdf_e1 dqic_sub(i) = - contfra(i) * pdf_e2 / pdf_shape dqtc_sub(i) = dcfc_sub(i) * qvapincld dnic_sub(i) = dcfc_sub(i) * Niceincld !--Add tendencies contfra(i) = MAX(0., contfra(i) + dcfc_sub(i)) qcont(i) = qcont(i) + dqtc_sub(i) + dqic_sub(i) qvcont(i) = qvcont(i) + dqtc_sub(i) Ncont(i) = Ncont(i) + dnic_sub(i) clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcfc_sub(i)) qclr(i) = qclr(i) - dqtc_sub(i) - dqic_sub(i) ELSEIF ( qvapincld_new .EQ. 0. ) THEN !--If all the ice has been sublimated, we sublimate !--completely the cloud and do not activate the dissipation !--process !--Tendencies and diagnostics dcfc_sub(i) = - contfra(i) dqtc_sub(i) = - qvcont(i) dqic_sub(i) = - ( qcont(i) - qvcont(i) ) dnic_sub(i) = - Ncont(i) !--Add tendencies contfra(i) = 0. qcont(i) = 0. qvcont(i) = 0. Ncont(i) = 0. clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcfc_sub(i)) qclr(i) = qclr(i) - dqtc_sub(i) - dqic_sub(i) ENDIF ! ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) ) .GT. qvapincld !------------------------------------ !-- PHASE ADJUSTMENT -- !------------------------------------ IF ( qvapincld_new .GT. 0. ) THEN !--Adjustment of the IWC to the new vapor in cloud !--(this can be either positive or negative) dqtc_adj(i) = ( qvapincld_new * contfra(i) - qvcont(i) ) dqic_adj(i) = - dqtc_adj(i) !--Add tendencies !--The vapor in the cloud is updated, but not qcont as it is constant !--through this process, as well as contfra which is unmodified qvcont(i) = MAX(0., MIN(qcont(i), qvcont(i) + dqtc_adj(i))) ENDIF ENDIF ! qiceincld .LT. eps ENDIF ! contfra(i) .GT. eps !------------------------------------------------------------------- !-- SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD -- !------------------------------------------------------------------- !--If there is a cloud IF ( cldfra(i) .GT. eps ) THEN qvapincld = qvc(i) / cldfra(i) IF ( qvapincld .GT. gamma_cond(i) * qsat(i) ) THEN qvapincld = gamma_cond(i) * qsat(i) qvc(i) = qvapincld * cldfra(i) ENDIF qiceincld = ( qcld(i) / cldfra(i) - qvapincld ) !--If the ice water content is too low, the cloud is purely sublimated !--Most probably, we advected a cloud with no ice water content (possible !--if the entire cloud precipited for example) IF ( qiceincld .LT. eps ) THEN dcf_sub(i) = - cldfra(i) dqvc_sub(i) = - qvc(i) dqi_sub(i) = - ( qcld(i) - qvc(i) ) cldfra(i) = 0. qcld(i) = 0. qvc(i) = 0. clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcf_sub(i)) qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i) !--Else, the cloud is adjusted and sublimated ELSE IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN IF ( qvapincld .GE. qsat(i) ) THEN !--If the cloud is initially supersaturated !--Exact explicit integration (qvc exact, qice explicit) tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub ) ELSE !--If the cloud is initially subsaturated !--Exact explicit integration (qice exact, qvc explicit) !--The barrier is set so that the resulting vapor in cloud !--cannot be greater than qsat !--qice_ratio is the ratio between the new ice content and !--the old one, it is comprised between 0 and 1 tauinv_depsub = qiceincld**(1./3.) * kappa_depsub qice_ratio = tauinv_depsub * dtime / 1.5 / qiceincld * ( qsat(i) - qvapincld ) !--The new vapor in the cloud is increased with the !--sublimated ice qvapincld_new = qvapincld + qiceincld * ( 1. - MAX(0., 1. - qice_ratio)**1.5 ) !--The new vapor in the cloud cannot be greater than qsat qvapincld_new = MIN(qvapincld_new, qsat(i)) !--If all the ice is sublimated IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0. ENDIF ! qvapincld .GT. qsat ELSE !--We keep the saturation adjustment hypothesis, and the vapor in the !--cloud is set equal to the saturation vapor IF ( ( qvapincld + qiceincld ) .GT. qsat(i) ) THEN qvapincld_new = qsat(i) ELSE qvapincld_new = 0. ENDIF ENDIF ! ok_unadjusted_clouds !------------------------------------ !-- DISSIPATION OF THE CLOUD -- !------------------------------------ !--If the dissipation process must be activated IF ( ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) ) .GT. qvapincld ) THEN !--Gamma distribution starting at qvapincld pdf_shape = nu_iwc_pdf_lscp / qiceincld pdf_y = pdf_shape * ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) - qvapincld ) pdf_e1 = GAMMAINC ( nu_iwc_pdf_lscp , pdf_y ) pdf_e2 = GAMMAINC ( nu_iwc_pdf_lscp + 1. , pdf_y ) !--Tendencies and diagnostics dcf_sub(i) = - cldfra(i) * pdf_e1 dqi_sub(i) = - cldfra(i) * pdf_e2 / pdf_shape dqvc_sub(i) = dcf_sub(i) * qvapincld !--Add tendencies cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i)) qcld(i) = qcld(i) + dqvc_sub(i) + dqi_sub(i) qvc(i) = qvc(i) + dqvc_sub(i) clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcf_sub(i)) qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i) ELSEIF ( qvapincld_new .EQ. 0. ) THEN !--If all the ice has been sublimated, we sublimate !--completely the cloud and do not activate the dissipation !--process !--Tendencies and diagnostics dcf_sub(i) = - cldfra(i) dqvc_sub(i) = - qvc(i) dqi_sub(i) = - ( qcld(i) - qvc(i) ) !--Add tendencies cldfra(i) = 0. qcld(i) = 0. qvc(i) = 0. clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcf_sub(i)) qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i) ENDIF ! ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) ) .GT. qvapincld !------------------------------------ !-- PHASE ADJUSTMENT -- !------------------------------------ IF ( qvapincld_new .GT. 0. ) THEN !--Adjustment of the IWC to the new vapor in cloud !--(this can be either positive or negative) dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) ) dqi_adj(i) = - dqvc_adj(i) !--Add tendencies !--The vapor in the cloud is updated, but not qcld as it is constant !--through this process, as well as cldfra which is unmodified qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i))) ENDIF ENDIF ! qiceincld .LT. eps ENDIF ! cldfra(i) .GT. eps !-------------------------------------------------------------------------- !-- CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS -- !-------------------------------------------------------------------------- !--This section relies on a distribution of water in the clear-sky region of !--the mesh. !--If there is a clear-sky region IF ( clrfra(i) .GT. eps ) THEN !--New PDF rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. !--rhl_clr is limited to 80%, because the fit is not valid above. !--This corresponds to rhi_clr ~= 140% at 210K, and ~= 105% at 245K rhl_clr = MAX(0., MIN(80., rhl_clr)) !--Calculation of the properties of the PDF !--Parameterization from IAGOS observations !--pdf_alpha, pdf_scale and pdf_gamma will be reused below !--Coefficient for standard deviation: !-- tuning coef * (clear sky area**0.25) * (function of temperature) !--MAX on the sqrt of the area because the parameterization was derived for !-- L > 50 km pdf_e1 = beta_pdf_lscp * SQRT(MAX(5e4, SQRT( clrfra(i) * cell_area(i) ) )) & * MAX( MAX(210., MIN(245., temp(i))) - temp_thresh_pdf_lscp, 0. ) IF ( rhl_clr .GT. 50. ) THEN pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp ELSE pdf_std = pdf_e1 * rhl_clr / 50. ENDIF pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * & MAX( temp_nowater - MAX(210., MIN(245., temp(i))), 0. ) pdf_alpha(i) = EXP( rhl_clr / 100. ) * pdf_e3 pdf_alpha(i) = MIN(10., pdf_alpha(i)) !--Avoid overflows !IF ( ok_warm_cloud ) THEN ! !--If the statistical scheme is activated, the standard deviation is adapted ! !--to depend on the pressure level. It is multiplied by ratqs, so that near the ! !--surface std is almost 0, and upper than about 450 hPa the std is left untouched ! pdf_std = pdf_std * ratqs(i) !ENDIF pdf_gamma(i) = GAMMA(1. + 1. / pdf_alpha(i)) !--Barrier to avoid overflows pdf_scale(i) = MAX(eps, MIN(rhl_clr / pdf_gamma(i), pdf_std / SQRT( & GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 ))) !--rhl_clr is not bounded for the location parameter rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) !--Calculation of the newly condensed water and fraction (pronostic) !--Integration of the clear sky PDF between gamma_cond*qsat and +inf !--NB. the calculated values are clear-sky averaged pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100. pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_nuc = 0. pdf_q_above_nuc = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_nuc = 1. pdf_q_above_nuc = qclr(i) / clrfra(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_nuc = EXP( - pdf_y ) pdf_q_above_nuc = ( pdf_e3 + pdf_loc * pdf_fra_above_nuc ) * qsatl(i) / 100. ENDIF IF ( pdf_fra_above_nuc .GT. eps ) THEN dcf_con(i) = clrfra(i) * pdf_fra_above_nuc dqt_con = clrfra(i) * pdf_q_above_nuc !--Barriers (should be useless dcf_con(i) = MIN(dcf_con(i), clrfra(i)) dqt_con = MIN(dqt_con, qclr(i)) IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute !--the new vapor qvapincld. The timestep is divided by two because we do not !--know when the condensation occurs qvapincld = gamma_cond(i) * qsat(i) qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i) tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) & * EXP( - dtime / 2. * tauinv_depsub ) ELSE !--We keep the saturation adjustment hypothesis, and the vapor in the !--newly formed cloud is set equal to the saturation vapor. qvapincld_new = qsat(i) ENDIF !--Tendency on cloud vapor and diagnostic dqvc_con(i) = qvapincld_new * dcf_con(i) dqi_con(i) = dqt_con - dqvc_con(i) !--Add tendencies cldfra(i) = cldfra(i) + dcf_con(i) qcld(i) = qcld(i) + dqt_con qvc(i) = qvc(i) + dqvc_con(i) clrfra(i) = clrfra(i) - dcf_con(i) qclr(i) = qclr(i) - dqt_con ENDIF ! pdf_fra_above_nuc .GT. eps ELSE !--Default value for the clear sky distribution: homogeneous distribution pdf_alpha(i) = 1. pdf_gamma(i) = 1. pdf_scale(i) = eps ENDIF ! clrfra(i) .GT. eps !-------------------------------------- !-- CLOUD MIXING -- !-------------------------------------- !--This process mixes the cloud with its surroundings: the subsaturated clear sky, !--and the supersaturated clear sky. It is activated if the cloud is big enough, !--but does not cover the entire mesh. ! IF ( ( cldfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN !-- PART 1 - TURBULENT DIFFUSION !--Clouds within the mesh are assumed to be ellipses. The length of the !--semi-major axis is a and the length of the semi-minor axis is b. !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer. !--In particular, it is 0 if cldfra = 1. !--clouds_perim is the total perimeter of the clouds within the mesh, !--not considering interfaces with other meshes (only the interfaces with clear !--sky are taken into account). !-- !--The area of each cloud is A = a * b * RPI, !--and the perimeter of each cloud is !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) ) !-- !--With cell_area the area of the cell, we have: !-- cldfra = A * N_cld_mix / cell_area !-- clouds_perim = P * N_cld_mix !-- bovera = aspect_ratio_cirrus !--P / a is a function of b / a only, that we can calculate !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) ) !--The clouds perimeter is imposed using the formula from Morcrette 2012, !--based on observations. !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra ) !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have: !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra ) !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) ) a_mix = Povera / coef_mixing_lscp / RPI / ( 1. - cldfra(i) ) / bovera !--and finally, !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a ) N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) ) * cell_area(i) & / Povera / a_mix !--The time required for turbulent diffusion to homogenize a region of size !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014) !--We compute L_mix and assume that the cloud is mixed over this length L_mix = SQRT( dtime**3 * pbl_eps(i) ) !--The mixing length cannot be greater than the semi-minor axis. In this case, !--the entire cloud is mixed. L_mix = MIN(L_mix, a_mix * bovera) !--The fraction of clear sky mixed is !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area clrfra_mix = N_cld_mix * RPI / cell_area(i) & * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 ) !--The fraction of clear sky mixed is !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area cldfra_mix = N_cld_mix * RPI / cell_area(i) & * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 ) !-- PART 2 - SHEARING !--The clouds are then sheared. We keep the shape and number !--assumptions from before. The clouds are sheared along their !--semi-major axis (a_mix), on the entire cell heigh dz. !--The increase in size is L_shear = coef_shear_lscp * shear(i) * dz * dtime !--therefore, the fraction of clear sky mixed is !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area !--and the fraction of cloud mixed is !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area !--(note that they are equal) shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i) !--and the environment and cloud mixed fractions are the same, !--which we add to the previous calculated mixed fractions. !--We therefore assume that the sheared clouds and the turbulent !--mixed clouds are different. clrfra_mix = clrfra_mix + shear_fra cldfra_mix = cldfra_mix + shear_fra !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES clrfra_mix = MIN(clrfra(i), clrfra_mix) cldfra_mix = MIN(cldfra(i), cldfra_mix) !--We compute the limit vapor in clear sky where the mixed cloud could not !--survive if all the ice crystals were sublimated. Note that here we assume, !--for growth or reduction of the cloud, that the mixed cloud is adjusted !--to saturation, ie the vapor in the mixed cloud is qsat. This is only a !--diagnostic, and if the cloud size is increased, we add the new vapor to the !--cloud's vapor without condensing or sublimating ice crystals IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN qiceinmix = ( qcld(i) - qvc(i) ) / cldfra(i) / ( 1. + clrfra_mix / cldfra_mix ) tauinv_depsub = qiceinmix**(1./3.) * kappa_depsub !qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime * tauinv_depsub ) ) qvapinmix_lim = qsat(i) - qiceinmix * MAX(1., 1.5 / ( dtime * tauinv_depsub )) qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) & - qvc(i) / cldfra(i) * cldfra_mix / clrfra_mix ELSE !--NB. if tau_depsub = 0 (ie tauinv_depsub = inf), we get the same result as above qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) & - qcld(i) / cldfra(i) * cldfra_mix / clrfra_mix ENDIF IF ( qvapinclr_lim .LT. 0. ) THEN !--Whatever we do, the cloud will increase in size dcf_mix(i) = clrfra_mix IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN dqvc_mix(i) = clrfra_mix * qclr(i) / clrfra(i) ELSE dqvc_mix(i) = clrfra_mix * qsat(i) dqi_mix(i) = clrfra_mix * ( qclr(i) / clrfra(i) - qsat(i) ) ENDIF ELSE !--We then calculate the clear sky part where the humidity is lower than !--qvapinclr_lim, and the part where it is higher than qvapinclr_lim !--This is the clear-sky PDF calculated in the condensation section. Note !--that if we are here, we necessarily went through the condensation part !--because the clear sky fraction can only be reduced by condensation. !--Thus the `pdf_xxx` variables are well defined. rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_x = qvapinclr_lim / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qsat(i) / qsatl(i) * 100. pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_lim = 0. pdf_q_above_lim = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_lim = clrfra(i) pdf_q_above_lim = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) pdf_q_above_lim = ( pdf_e3 * clrfra(i) & + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. ENDIF pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim !--sigma_mix is the ratio of the surroundings of the clouds where mixing !--increases the size of the cloud, to the total surroundings of the clouds. !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing !--decreases the size of the clouds !sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim ) !--The chance that the air surrounding cirrus clouds is moist is higher then !--the actual distribution of humidity in the gridbox. !--This is quantified with chi_mixing. For chi_mixing=1, the two are equal. !--If chi_mixing > 1, the chance that moist air is around cirrus clouds is increased sigma_mix = MAX(1. - pdf_fra_below_lim / clrfra_mix, & MIN(pdf_fra_above_lim / clrfra_mix, & chi_mixing * pdf_fra_above_lim & / ( pdf_fra_below_lim + chi_mixing * pdf_fra_above_lim ))) IF ( pdf_fra_above_lim .GT. eps ) THEN IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN dcf_mix(i) = clrfra_mix * sigma_mix dqvc_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim ELSE dcf_mix(i) = clrfra_mix * sigma_mix dqvc_mix(i) = clrfra_mix * sigma_mix * qsat(i) dqi_mix(i) = clrfra_mix * sigma_mix & * (pdf_q_above_lim / pdf_fra_above_lim - qsat(i)) ENDIF ENDIF IF ( pdf_fra_below_lim .GT. eps ) THEN dcf_mix(i) = dcf_mix(i) - cldfra_mix * ( 1. - sigma_mix ) dqvc_mix(i) = dqvc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & * qvc(i) / cldfra(i) dqi_mix(i) = dqi_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & * ( qcld(i) - qvc(i) ) / cldfra(i) ENDIF ENDIF ENDIF ! ( cldfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) !-------------------------------------- !-- CONTRAIL MIXING -- !-------------------------------------- IF ( ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN !-- PART 1 - TURBULENT DIFFUSION !--Clouds within the mesh are assumed to be ellipses. The length of the !--semi-major axis is a and the length of the semi-minor axis is b. !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer. !--In particular, it is 0 if cldfra = 1. !--clouds_perim is the total perimeter of the clouds within the mesh, !--not considering interfaces with other meshes (only the interfaces with clear !--sky are taken into account). !-- !--The area of each cloud is A = a * b * RPI, !--and the perimeter of each cloud is !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) ) !-- !--With cell_area the area of the cell, we have: !-- cldfra = A * N_cld_mix / cell_area !-- clouds_perim = P * N_cld_mix !-- bovera = aspect_ratio_contrails !--P / a is a function of b / a only, that we can calculate !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) ) !--The clouds perimeter is imposed using the formula from Morcrette 2012, !--based on observations. !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra ) !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have: !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra ) !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) ) a_mix = Povera / coef_mixing_contrails / RPI / ( 1. - contfra(i) ) / bovera !--and finally, !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a ) N_cld_mix = coef_mixing_contrails * contfra(i) * ( 1. - contfra(i) ) & * cell_area(i) / Povera / a_mix !--The time required for turbulent diffusion to homogenize a region of size !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014) !--We compute L_mix and assume that the cloud is mixed over this length L_mix = SQRT( dtime**3 * pbl_eps(i) ) !--The mixing length cannot be greater than the semi-minor axis. In this case, !--the entire cloud is mixed. L_mix = MIN(L_mix, a_mix * bovera) !--The fraction of clear sky mixed is !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area clrfra_mix = N_cld_mix * RPI / cell_area(i) & * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 ) !--The fraction of clear sky mixed is !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area cldfra_mix = N_cld_mix * RPI / cell_area(i) & * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 ) !-- PART 2 - SHEARING !--The clouds are then sheared. We keep the shape and number !--assumptions from before. The clouds are sheared with a random orientation !--of the wind, on average we assume that the wind and the semi-major axis !--make a 45 degrees angle. Moreover, the contrails only mix !--along their semi-minor axis (b), because it is easier to compute. !--With this, the clouds increase in size along b only, by a factor !--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind) L_shear = coef_shear_contrails * shear(i) * dz * dtime !--therefore, the fraction of clear sky mixed is !-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area !--and the fraction of cloud mixed is !-- N_cld_mix * ( a * b - a * (b - L_shear * SQRT(2.) / 2.) ) * RPI / 2. / cell_area !--(note that they are equal) shear_fra = RPI * L_shear * a_mix * SQRT(2.) / 2. / 2. * N_cld_mix / cell_area(i) !shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i) !--and the environment and cloud mixed fractions are the same, !--which we add to the previous calculated mixed fractions. !--We therefore assume that the sheared clouds and the turbulent !--mixed clouds are different. clrfra_mix = clrfra_mix + shear_fra cldfra_mix = cldfra_mix + shear_fra !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES clrfra_mix = MIN(clrfra(i), clrfra_mix) cldfra_mix = MIN(contfra(i), cldfra_mix) !--We compute the limit vapor in clear sky where the mixed cloud could not !--survive if all the ice crystals were sublimated. Note that here we assume, !--for growth or reduction of the cloud, that the mixed cloud is adjusted !--to saturation, ie the vapor in the mixed cloud is qsat. This is only a !--diagnostic, and if the cloud size is increased, we add the new vapor to the !--cloud's vapor without condensing or sublimating ice crystals IF ( ok_unadjusted_contrails ) THEN qiceinmix = ( qcont(i) - qvcont(i) ) / contfra(i) / ( 1. + clrfra_mix / cldfra_mix ) tauinv_depsub = qiceinmix**(1./3.) * kappa_depsub & * (Niceincld / N_ice_volume)**(2./3.) !qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime * tauinv_depsub ) ) qvapinmix_lim = qsat(i) - qiceinmix * MAX(1., 1.5 / ( dtime * tauinv_depsub )) qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) & - qvcont(i) / contfra(i) * cldfra_mix / clrfra_mix ELSE !--NB. if tau_depsub = 0 (ie tauinv_depsub = inf), we get the same result as above qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) & - qcont(i) / contfra(i) * cldfra_mix / clrfra_mix ENDIF IF ( qvapinclr_lim .LT. 0. ) THEN !--Whatever we do, contrails will increase in size dcfc_mix(i) = dcfc_mix(i) + clrfra_mix IF ( ok_unadjusted_contrails ) THEN dqtc_mix(i) = clrfra_mix * qclr(i) / clrfra(i) ELSE dqtc_mix(i) = clrfra_mix * qsat(i) dqic_mix(i) = clrfra_mix * ( qclr(i) / clrfra(i) - qsat(i) ) ENDIF ELSE !--We then calculate the clear sky part where the humidity is lower than !--qvapinclr_lim, and the part where it is higher than qvapinclr_lim !--This is the clear-sky PDF calculated in the condensation section. Note !--that if we are here, we necessarily went through the condensation part !--because the clear sky fraction can only be reduced by condensation. !--Thus the `pdf_xxx` variables are well defined. rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_x = qvapinclr_lim / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qsat(i) / qsatl(i) * 100. pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_lim = 0. pdf_q_above_lim = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_lim = clrfra(i) pdf_q_above_lim = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) pdf_q_above_lim = ( pdf_e3 * clrfra(i) & + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. ENDIF pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim !--sigma_mix is the ratio of the surroundings of the clouds where mixing !--increases the size of the cloud, to the total surroundings of the clouds. !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing !--decreases the size of the clouds !--For aviation, we increase the chance that the air surrounding contrails !--is moist. This is quantified with chi_mixing_contrails sigma_mix = MAX(1. - pdf_fra_below_lim / clrfra_mix, & MIN(pdf_fra_above_lim / clrfra_mix, & chi_mixing_contrails * pdf_fra_above_lim & / ( pdf_fra_below_lim + chi_mixing_contrails * pdf_fra_above_lim ))) IF ( pdf_fra_above_lim .GT. eps ) THEN IF ( ok_unadjusted_contrails ) THEN dcfc_mix(i) = clrfra_mix * sigma_mix dqtc_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim ELSE dcfc_mix(i) = clrfra_mix * sigma_mix dqtc_mix(i) = clrfra_mix * sigma_mix * qsat(i) dqic_mix(i) = clrfra_mix * sigma_mix & * ( pdf_q_above_lim / pdf_fra_above_lim - qsat(i) ) ENDIF ENDIF IF ( pdf_fra_below_lim .GT. eps ) THEN dcfc_mix(i) = dcfc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) dqtc_mix(i) = dqtc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & * qvcont(i) / contfra(i) dqic_mix(i) = dqic_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & * ( qcont(i) - qvcont(i) ) / contfra(i) dnic_mix(i) = dnic_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & * Ncont(i) / contfra(i) ENDIF ENDIF ENDIF ! ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) IF ( contfra(i) .GT. eps ) THEN !--We balance the mixing tendencies between the different cloud classes mixed_fraction = dcf_mix(i) + dcfc_mix(i) IF ( mixed_fraction .GT. clrfra(i) ) THEN mixed_fraction = clrfra(i) / mixed_fraction dcf_mix(i) = dcf_mix(i) * mixed_fraction dqvc_mix(i) = dqvc_mix(i) * mixed_fraction dqi_mix(i) = dqi_mix(i) * mixed_fraction dcfc_mix(i) = dcfc_mix(i) * mixed_fraction dqtc_mix(i) = dqtc_mix(i) * mixed_fraction dqic_mix(i) = dqic_mix(i) * mixed_fraction dnic_mix(i) = dnic_mix(i) * mixed_fraction ENDIF !--Add tendencies contfra(i) = contfra(i) + dcfc_mix(i) qcont(i) = qcont(i) + dqtc_mix(i) + dqic_mix(i) qvcont(i) = qvcont(i) + dqtc_mix(i) Ncont(i) = Ncont(i) + dnic_mix(i) clrfra(i) = clrfra(i) - dcfc_mix(i) qclr(i) = qclr(i) - dqtc_mix(i) ENDIF !--Add tendencies cldfra(i) = cldfra(i) + dcf_mix(i) qcld(i) = qcld(i) + dqvc_mix(i) + dqi_mix(i) qvc(i) = qvc(i) + dqvc_mix(i) clrfra(i) = clrfra(i) - dcf_mix(i) qclr(i) = qclr(i) - dqvc_mix(i) - dqi_mix(i) !---------------------------------------- !-- ICE CRYSTALS AGGREGATION -- !---------------------------------------- !--Aggregation of ice crystals only occur for 2-moment microphysical clouds, !--i.e., contrails for now IF ( ok_plane_contrail ) THEN IF ( contfra(i) .GT. eps ) THEN mice = ( qcont(i) - qvcont(i) ) / MAX(eps, Ncont(i)) icefall_velo = ICECRYSTAL_VELO(mice, temp(i), pplay(i)) !--Effective radius (in meters) dei = (mice / rho_ice * 3. / 4. / RPI)**(1./3.) / eff2vol_radius_contrails !--Effective radius to effective diameter dei = 8. / 3. / SQRT(3.) * dei !--Lin 1983's formula sticking_coef = EXP( 0.025 * MIN( ( temp(i) - RTT ), 0.) ) gamma_agg = 1. dispersion_coef = 0.1 dnic_agg(i) = - 0.5 * gamma_agg * RPI / 6. * sticking_coef * dispersion_coef & * dei**2 * icefall_velo * ( Ncont(i) / contfra(i) )**2 !--Gridbox average dnic_agg(i) = dnic_agg(i) * contfra(i) !--Add tendencies Ncont(i) = Ncont(i) + dnic_agg(i) ENDIF ENDIF !--------------------------------------- !-- ICE SEDIMENTATION -- !--------------------------------------- IF ( ok_ice_sedim ) THEN !--Initialisation dz_auto = 0. dz_auto_cont = 0. dcf_sed1 = 0. dqvc_sed1 = 0. dqi_sed1 = 0. dcf_sed2 = 0. dqvc_sed2 = 0. dqi_sed2 = 0. dcfc_sed1 = 0. dqtc_sed1 = 0. dqic_sed1 = 0. dnic_sed1 = 0. dcfc_sed2 = 0. dqtc_sed2 = 0. dqic_sed2 = 0. dnic_sed2 = 0. ! !--First, the current ice is sedimentated into the layer below. The ice fallspeed !--velocity is calculated and the quantity of sedimentated ice is computed !--accordingly. The decrease in cloud fraction is calculated such that !--in-cloud ice water content is constant. ! IF ( cldfra(i) .GT. eps ) THEN iwc = rho(i) * ( qcld(i) - qvc(i) ) / cldfra(i) icefall_velo = ffallv_sed * cice_velo * MAX(0., iwc)**dice_velo !--Sedimentation coef_sed = MAX(0., MIN(1., 2. - icefall_velo * dtime / dz)) dzsed(i) = MIN(dz, icefall_velo * dtime) * coef_sed cfsed(i) = cldfra(i) qice_sedim = ( qcld(i) - qvc(i) ) * dzsed(i) / dz !--Tendencies dcf_sed1 = - cldfra(i) * dzsed(i) / dz dqi_sed1 = - qice_sedim dqvc_sed1 = - qvc(i) * dzsed(i) / dz !--Convert to flux flsed(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime !--Autoconversion coef_auto = MAX(0., MIN(1., icefall_velo * dtime / dz - 1.)) dz_auto = MIN(dz, icefall_velo * dtime) * coef_auto qice_auto = ( qcld(i) - qvc(i) ) * dz_auto / dz !--Tendencies dcf_auto(i) = - cldfra(i) * dz_auto / dz dqi_auto(i) = - qice_auto dqvc_auto(i) = - qvc(i) * dz_auto / dz !--Convert to flux flauto(i) = flauto(i) + qice_auto * ( paprsdn(i) - paprsup(i) ) / RG / dtime !--Save the vapor in the cloud (will be needed later) qvapincld = qvc(i) / cldfra(i) !--Add tendencies cldfra(i) = cldfra(i) + dcf_sed1 + dcf_auto(i) qcld(i) = qcld(i) + dqvc_sed1 + dqi_sed1 + dqvc_auto(i) + dqi_auto(i) qvc(i) = qvc(i) + dqvc_sed1 + dqvc_auto(i) clrfra(i) = clrfra(i) - dcf_sed1 - dcf_auto(i) qclr(i) = qclr(i) - dqvc_sed1 - dqi_sed1 - dqvc_auto(i) - dqi_auto(i) ENDIF ! IF ( contfra(i) .GT. eps ) THEN mice = ( qcont(i) - qvcont(i) ) / MAX(eps, Ncont(i)) icefall_velo = ICECRYSTAL_VELO(mice, temp(i), pplay(i)) !--Sedimentation coef_sed = MAX(0., MIN(1., 2. - icefall_velo * dtime / dz)) dzsed_cont(i) = MIN(dz, icefall_velo * dtime) * coef_sed cfsed_cont(i) = contfra(i) qice_sedim = ( qcont(i) - qvcont(i) ) * dzsed_cont(i) / dz Nice_sedim = Ncont(i) * dzsed_cont(i) / dz !--Tendencies dcfc_sed1 = - contfra(i) * dzsed_cont(i) / dz dqic_sed1 = - qice_sedim dqtc_sed1 = - qvcont(i) * dzsed_cont(i) / dz dnic_sed1 = - Nice_sedim !--Convert to flux flsed_cont(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime Nflsed_cont(i) = Nice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime !--Autoconversion coef_auto = MAX(0., MIN(1., icefall_velo * dtime / dz - 1.)) dz_auto_cont = MIN(dz, icefall_velo * dtime) * coef_auto qice_auto = ( qcont(i) - qvcont(i) ) * dz_auto_cont / dz !--Tendencies dcfc_auto(i) = - contfra(i) * dz_auto_cont / dz dqic_auto(i) = - qice_auto dqtc_auto(i) = - qvcont(i) * dz_auto_cont / dz dnic_auto(i) = - Ncont(i) * dz_auto_cont / dz !--Convert to flux flauto(i) = flauto(i) + qice_auto * ( paprsdn(i) - paprsup(i) ) / RG / dtime !--Save the vapor in the cloud (will be needed later) qvapincont = qvcont(i) / contfra(i) !--Add tendencies contfra(i) = contfra(i) + dcfc_sed1 + dcfc_auto(i) qcont(i) = qcont(i) + dqtc_sed1 + dqtc_auto(i) + dqic_sed1 + dqic_auto(i) qvcont(i) = qvcont(i) + dqtc_sed1 + dqtc_auto(i) Ncont(i) = Ncont(i) + dnic_sed1 + dnic_auto(i) clrfra(i) = clrfra(i) - dcfc_sed1 - dcfc_auto(i) qclr(i) = qclr(i) - dqtc_sed1 - dqtc_auto(i) - dqic_sed1 - dqic_auto(i) ENDIF ! !--Then, the sedimentated ice from above is added to the gridbox !--The falling ice is partly considered a new cloud in the gridbox. !--BEWARE with this parameterisation, we can create a new cloud with a much !--different ice water content and water vapor content than the existing cloud !--(if it exists). This implies than unphysical fluxes of ice and vapor !--occur between the existing cloud and the newly formed cloud (from sedimentation). !--Note also that currently, the precipitation scheme assume a maximum !--random overlap, meaning that the new formed clouds will not be affected !--by vertical wind shear. ! sedfra_abv = MIN(dzsed_abv(i), dz) / dz * cfsed_abv(i) IF ( sedfra_abv .GT. eps ) THEN !--Three regions to be defined: (1) clear-sky, (2) cloudy, and !--(3) region where it was previously cloudy but now clear-sky because of !--ice sedimentation !--(1) we use the clear-sky PDF to determine the humidity and increase the !--cloud fraction / sublimate falling ice. NB it can also fall into !--another cloud class !--(2) we do not increase cloud fraction and add the falling ice to the cloud !--(3) we increase the cloud fraction but use in-cloud water vapor as !--background vapor sedfra2 = MIN(cfsed(i), cfsed_abv(i)) & * MAX(0., MIN(dz, dzsed_abv(i)) - dzsed(i) - dz_auto) / dz sedfra3 = MIN(cfsed(i), cfsed_abv(i)) & * MIN(MIN(dz, dzsed_abv(i)), dzsed(i) + dz_auto) / dz sedfra1 = sedfra_abv - sedfra3 - sedfra2 qiceincld = qised_abv(i) / sedfra_abv !--For region (1) IF ( sedfra1 .GT. eps ) THEN IF ( clrfra(i) .GT. eps ) THEN IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN tauinv_depsub = qiceincld**(1./3.) * kappa_depsub qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) ) ELSE qvapinclr_lim = qsat(i) - qiceincld ENDIF rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_x = qvapinclr_lim / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qsat(i) / qsatl(i) * 100. pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_lim = 0. pdf_q_above_lim = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_lim = clrfra(i) pdf_q_above_lim = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) pdf_q_above_lim = ( pdf_e3 * clrfra(i) & + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. ENDIF pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim !--The first three lines allow to favor ISSR rather than subsaturated sky for !--the growth. The fourth line indicates that the ice is falling only in !--the clear-sky area. The rest falls into other cloud classes sedfra_tmp = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, & sedfra1 * chi_sedim * pdf_fra_above_lim & / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) & * clrfra(i) / (totfra_in(i) - cldfra(i)) IF ( pdf_fra_above_lim .GT. eps ) THEN IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN dcf_sed2 = dcf_sed2 + sedfra_tmp dqvc_sed2 = dqvc_sed2 + sedfra_tmp * pdf_q_above_lim / pdf_fra_above_lim dqi_sed2 = dqi_sed2 + sedfra_tmp * qiceincld ELSE dcf_sed2 = dcf_sed2 + sedfra_tmp dqvc_sed2 = dqvc_sed2 + sedfra_tmp * qsat(i) dqi_sed2 = dqi_sed2 + sedfra_tmp & * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i)) ENDIF ENDIF ENDIF ! clrfra(i) .GT. eps !--If the sedimented ice falls into a contrail, it increases its content IF ( contfra(i) .GT. eps ) THEN sedfra_tmp = sedfra1 * contfra(i) / (totfra_in(i) - cldfra(i)) dqic_sed2 = dqic_sed2 + sedfra_tmp * qiceincld ENDIF ENDIF ! sedfra1 .GT. eps !--For region (2) dqi_sed2 = dqi_sed2 + sedfra2 * qiceincld !--For region (3) IF ( cldfra(i) .GT. eps ) THEN dcf_sed2 = dcf_sed2 + sedfra3 dqvc_sed2 = dqvc_sed2 + sedfra3 * qvapincld dqi_sed2 = dqi_sed2 + sedfra3 * qiceincld ENDIF ENDIF ! sedfra_abv .GT. eps IF ( ok_plane_contrail ) THEN sedfra_abv = MIN(dzsed_cont_abv(i), dz) / dz * cfsed_cont_abv(i) IF ( sedfra_abv .GT. eps ) THEN !--Three regions to be defined: (1) clear-sky, (2) cloudy, and !--(3) region where it was previously cloudy but now clear-sky because of !--ice sedimentation !--(1) we use the clear-sky PDF to determine the humidity and increase the !--cloud fraction / sublimate falling ice. NB it can also fall into !--another cloud class !--(2) we do not increase cloud fraction and add the falling ice to the cloud !--(3) we increase the cloud fraction but use in-cloud water vapor as !--background vapor sedfra2 = MIN(cfsed_cont(i), cfsed_cont_abv(i)) & * MAX(0., MIN(dz, dzsed_cont_abv(i)) & - dzsed_cont(i) - dz_auto_cont) / dz sedfra3 = MIN(cfsed_cont(i), cfsed_cont_abv(i)) & * MIN(MIN(dz, dzsed_cont_abv(i)), & dzsed_cont(i) + dz_auto_cont) / dz sedfra1 = sedfra_abv - sedfra3 - sedfra2 qiceincld = qised_cont_abv(i) / sedfra_abv Niceincld = Nised_cont_abv(i) / sedfra_abv !--For region (1) IF ( sedfra1 .GT. eps ) THEN IF ( clrfra(i) .GT. eps ) THEN IF ( ok_unadjusted_contrails ) THEN tauinv_depsub = qiceincld**(1./3.) * kappa_depsub & * (Niceincld / N_ice_volume)**(2./3.) qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) ) ELSE qvapinclr_lim = qsat(i) - qiceincld ENDIF rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_x = qvapinclr_lim / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qsat(i) / qsatl(i) * 100. pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows pdf_fra_above_lim = 0. pdf_q_above_lim = 0. ELSEIF ( pdf_y .LT. -10. ) THEN pdf_fra_above_lim = clrfra(i) pdf_q_above_lim = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) pdf_q_above_lim = ( pdf_e3 * clrfra(i) & + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. ENDIF pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim !--The first three lines allow to favor ISSR rather than subsaturated sky for !--the growth. The fourth line indicates that the ice is falling only in !--the clear-sky area. The rest falls into other cloud classes sedfra_tmp = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, & sedfra1 * chi_sedim * pdf_fra_above_lim & / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) & * clrfra(i) / (totfra_in(i) - contfra(i)) IF ( pdf_fra_above_lim .GT. eps ) THEN IF ( ok_unadjusted_contrails ) THEN dcfc_sed2 = dcfc_sed2 + sedfra_tmp dqtc_sed2 = dqtc_sed2 + sedfra_tmp * pdf_q_above_lim / pdf_fra_above_lim dqic_sed2 = dqic_sed2 + sedfra_tmp * qiceincld dnic_sed2 = dnic_sed2 + sedfra_tmp * Niceincld ELSE dcfc_sed2 = dcfc_sed2 + sedfra_tmp dqtc_sed2 = dqtc_sed2 + sedfra_tmp * qsat(i) dqic_sed2 = dqic_sed2 + sedfra_tmp & * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i)) dnic_sed2 = dnic_sed2 + sedfra_tmp * Niceincld ENDIF ENDIF ENDIF ! clrfra(i) .GT. eps !--If the sedimented ice falls into a natural cirrus, it increases its content IF ( cldfra(i) .GT. eps ) THEN sedfra_tmp = sedfra1 * cldfra(i) / (totfra_in(i) - contfra(i)) dqi_sed2 = dqi_sed2 + sedfra_tmp * qiceincld ENDIF ENDIF ! sedfra1 .GT. eps !--For region (2) dqic_sed2 = dqic_sed2 + sedfra2 * qiceincld dnic_sed2 = dnic_sed2 + sedfra2 * Niceincld !--For region (3) IF ( contfra(i) .GT. eps ) THEN dcfc_sed2 = dcfc_sed2 + sedfra3 dqtc_sed2 = dqtc_sed2 + sedfra3 * qvapincont dqic_sed2 = dqic_sed2 + sedfra3 * qiceincld dnic_sed2 = dnic_sed2 + sedfra3 * Niceincld ENDIF ENDIF ! sedfra_abv .GT. eps ENDIF ! ok_plane_contrails !--Add tendencies cldfra(i) = cldfra(i) + dcf_sed2 qcld(i) = qcld(i) + dqvc_sed2 + dqi_sed2 qvc(i) = qvc(i) + dqvc_sed2 clrfra(i) = clrfra(i) - dcf_sed2 qclr(i) = qclr(i) - dqvc_sed2 - dqi_sed2 !--We re-include sublimated sedimentated ice into clear sky water vapor qclr(i) = qclr(i) + qised_abv(i) !--Diagnostics dcf_sed(i) = dcf_sed1 + dcf_sed2 dqvc_sed(i) = dqvc_sed1 + dqvc_sed2 dqi_sed(i) = dqi_sed1 + dqi_sed2 if ( ok_plane_contrail ) THEN !--Add tendencies contfra(i) = contfra(i) + dcfc_sed2 qcont(i) = qcont(i) + dqtc_sed2 + dqic_sed2 qvcont(i) = qvcont(i) + dqtc_sed2 Ncont(i) = Ncont(i) + dnic_sed2 clrfra(i) = clrfra(i) - dcfc_sed2 qclr(i) = qclr(i) - dqtc_sed2 - dqic_sed2 !--We re-include sublimated sedimentated ice into clear sky water vapor qclr(i) = qclr(i) + qised_cont_abv(i) !--Diagnostics dcfc_sed(i) = dcfc_sed1 + dcfc_sed2 dqtc_sed(i) = dqtc_sed1 + dqtc_sed2 dqic_sed(i) = dqic_sed1 + dqic_sed2 dnic_sed(i) = dnic_sed1 + dnic_sed2 ENDIF ENDIF ! ok_ice_sedim !--We put back contrails in the clouds class IF ( contfra(i) .GT. 0. ) THEN cldfra(i) = cldfra(i) + contfra(i) qcld(i) = qcld(i) + qcont(i) qvc(i) = qvc(i) + qvcont(i) ENDIF !--Diagnose ISSRs IF ( clrfra(i) .GT. eps ) THEN rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qsat(i) / qsatl(i) * 100. pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows issrfra(i) = 0. qissr(i) = 0. ELSEIF ( pdf_y .LT. -10. ) THEN issrfra(i) = clrfra(i) qissr(i) = qclr(i) ELSE pdf_y = EXP( pdf_y ) pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) issrfra(i) = EXP( - pdf_y ) * clrfra(i) qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100. ENDIF IF ( issrfra(i) .LT. eps ) THEN issrfra(i) = 0. qissr(i) = 0. ENDIF ELSE issrfra(i) = 0. qissr(i) = 0. ENDIF !------------------------------------------- !-- FINAL BARRIERS AND OUTPUTS -- !------------------------------------------- cldfra(i) = MIN(cldfra(i), totfra_in(i)) qcld(i) = MIN(qcld(i), qtot_in(i)) qvc(i) = MIN(qvc(i), qcld(i)) IF ( cldfra(i) .LT. eps ) THEN !--If the cloud is too small, it is sublimated. cldfra(i) = 0. qcld(i) = 0. qvc(i) = 0. qincld(i) = qsat(i) ELSE qincld(i) = qcld(i) / cldfra(i) ENDIF ! cldfra .LT. eps !--Diagnostics dcf_sub(i) = dcf_sub(i) / dtime dcf_con(i) = dcf_con(i) / dtime dcf_mix(i) = dcf_mix(i) / dtime dcf_sed(i) = dcf_sed(i) / dtime dcf_auto(i) = dcf_auto(i) / dtime dqi_adj(i) = dqi_adj(i) / dtime dqi_sub(i) = dqi_sub(i) / dtime dqi_con(i) = dqi_con(i) / dtime dqi_mix(i) = dqi_mix(i) / dtime dqi_sed(i) = dqi_sed(i) / dtime dqi_auto(i) = dqi_auto(i) / dtime dqvc_adj(i) = dqvc_adj(i) / dtime dqvc_sub(i) = dqvc_sub(i) / dtime dqvc_con(i) = dqvc_con(i) / dtime dqvc_mix(i) = dqvc_mix(i) / dtime dqvc_sed(i) = dqvc_sed(i) / dtime dqvc_auto(i) = dqvc_auto(i) / dtime ENDIF ! pt_pron_clds(i) ENDIF ! end keepgoing ENDDO ! end loop on i !---------------------------------------- !-- CONTRAILS AND AVIATION -- !---------------------------------------- IF ( ok_plane_contrail ) THEN rho = pplay / temp / RD CALL contrails_formation( & klon, dtime, pplay, temp, rho, qsat, qsatl, gamma_cond, & flight_dist, flight_fuel, clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, & keepgoing, pt_pron_clds, & Tcritcont, qcritcont, potcontfraP, potcontfraNP, & AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails, & dcfc_ini, dqic_ini, dqtc_ini, dnic_ini) DO i = 1, klon IF ( keepgoing(i) .AND. pt_pron_clds(i) ) THEN dcfc_ini(i) = MIN(dcfc_ini(i), issrfra(i)) dqtc_ini(i) = MIN(dqtc_ini(i), qissr(i)) dqic_ini(i) = MIN(dqic_ini(i) + dqtc_ini(i), qissr(i)) - dqtc_ini(i) !--Add tendencies cldfra(i) = cldfra(i) + dcfc_ini(i) qcld(i) = qcld(i) + dqtc_ini(i) + dqic_ini(i) qvc(i) = qvc(i) + dqtc_ini(i) issrfra(i) = issrfra(i) - dcfc_ini(i) qissr(i) = qissr(i) - dqtc_ini(i) - dqic_ini(i) clrfra(i) = clrfra(i) - dcfc_ini(i) qclr(i) = qclr(i) - dqtc_ini(i) - dqic_ini(i) contfra(i) = contfra(i) + dcfc_ini(i) qcont(i) = qcont(i) + dqtc_ini(i) + dqic_ini(i) qvcont(i) = qvcont(i) + dqtc_ini(i) Ncont(i) = Ncont(i) + dnic_ini(i) !------------------------------------------- !-- FINAL BARRIERS AND OUTPUTS -- !------------------------------------------- IF ( (contfra(i) .LT. eps) .OR. (qcont(i) .LT. qvcont(i)) .OR. (Ncont(i) .LT. eps) ) THEN cldfra(i) = cldfra(i) - contfra(i) qcld(i) = qcld(i) - qcont(i) qvc(i) = qvc(i) - qvcont(i) contfra(i) = 0. qcont(i) = 0. qvcont(i) = 0. Ncont(i) = 0. ENDIF IF ( cldfra(i) .LT. eps ) THEN !--If the cloud is too small, it is sublimated. cldfra(i) = 0. qcld(i) = 0. qvc(i) = 0. contfra(i)= 0. qcont(i) = 0. qvcont(i) = 0. Ncont(i) = 0. qincld(i) = qsat(i) ELSE qincld(i) = qcld(i) / cldfra(i) ENDIF cldfra(i) = MAX(cldfra(i), contfra(i)) qcld(i) = MAX(qcld(i), qcont(i)) qvc(i) = MAX(qvc(i), qvcont(i)) !--Diagnostics dcfc_ini(i) = dcfc_ini(i) / dtime dqic_ini(i) = dqic_ini(i) / dtime dqtc_ini(i) = dqtc_ini(i) / dtime dnic_ini(i) = dnic_ini(i) / dtime dcfc_sub(i) = dcfc_sub(i) / dtime dqic_sub(i) = dqic_sub(i) / dtime dqtc_sub(i) = dqtc_sub(i) / dtime dnic_sub(i) = dnic_sub(i) / dtime dcfc_mix(i) = dcfc_mix(i) / dtime dqic_mix(i) = dqic_mix(i) / dtime dqtc_mix(i) = dqtc_mix(i) / dtime dnic_mix(i) = dnic_mix(i) / dtime dnic_agg(i) = dnic_agg(i) / dtime dcfc_sed(i) = dcfc_sed(i) / dtime dqic_sed(i) = dqic_sed(i) / dtime dqtc_sed(i) = dqtc_sed(i) / dtime dnic_sed(i) = dnic_sed(i) / dtime dcfc_auto(i) = dcfc_auto(i) / dtime dqic_auto(i) = dqic_auto(i) / dtime dqtc_auto(i) = dqtc_auto(i) / dtime dnic_auto(i) = dnic_auto(i) / dtime ENDIF ! keepgoing ENDDO ! loop on klon ENDIF ! ok_plane_contrail END SUBROUTINE condensation_ice_supersat !********************************************************************************** !********************************************************************************** SUBROUTINE condensation_cloudth(klon, & & temp,qt,qt_th,frac_th,zpspsk,play,thetal_th, & & ratqs,sigma_qtherm,qsth,qsenv,qcloud,ctot,ctotth,ctot_vol, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) ! This routine computes the condensation of clouds in convective boundary layers ! with thermals assuming two separate distribution of the saturation deficit in ! the thermal plumes and in the environment ! It is based on the work of Arnaud Jam (Jam et al. 2013, BLM) ! Author : Etienne Vignon (LMDZ/CNRS) ! Date: February 2025 ! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam ! IMPORTANT NOTE: we assume iflag_cloudth_vert=7 !----------------------------------------------------------------------------------- use lmdz_lscp_ini, only: iflag_cloudth_vert,iflag_ratqs,iflag_cloudth_vert_noratqs use lmdz_lscp_ini, only: vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin use lmdz_lscp_ini, only: RTT, RG, RPI, RD, RV, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th IMPLICIT NONE !------------------------------------------------------------------------------ ! Declarations !------------------------------------------------------------------------------ ! INPUT/OUTPUT INTEGER, INTENT(IN) :: klon REAL, DIMENSION(klon), INTENT(IN) :: temp ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip REAL, DIMENSION(klon), INTENT(IN) :: qt ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip REAL, DIMENSION(klon), INTENT(IN) :: qt_th ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip REAL, DIMENSION(klon), INTENT(IN) :: thetal_th ! Liquid potential temperature in thermals [K]: has not seen the evap of precip REAL, DIMENSION(klon), INTENT(IN) :: frac_th ! Fraction of the mesh covered by thermals [0-1] REAL, DIMENSION(klon), INTENT(IN) :: zpspsk ! Exner potential REAL, DIMENSION(klon), INTENT(IN) :: play ! Pressure of layers [Pa] REAL, DIMENSION(klon), INTENT(IN) :: ratqs ! Parameter that determines the width of the water distrib [-] REAL, DIMENSION(klon), INTENT(IN) :: sigma_qtherm ! Parameter determining the width of the distrib in thermals [-] REAL, DIMENSION(klon), INTENT(IN) :: qsth ! Saturation specific humidity in thermals REAL, DIMENSION(klon), INTENT(IN) :: qsenv ! Saturation specific humidity in environment REAL, DIMENSION(klon), INTENT(INOUT) :: ctot ! Cloud fraction [0-1] REAL, DIMENSION(klon), INTENT(INOUT) :: ctotth ! Cloud fraction [0-1] in thermals REAL, DIMENSION(klon), INTENT(INOUT) :: ctot_vol ! Volume cloud fraction [0-1] REAL, DIMENSION(klon), INTENT(INOUT) :: qcloud ! In cloud total water content [kg/kg] REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment ! LOCAL VARIABLES INTEGER itap,ind1,l,ig,iter,k INTEGER iflag_topthermals, niter REAL qcth(klon) REAL qcenv(klon) REAL qctot(klon) REAL cth(klon) REAL cenv(klon) REAL cth_vol(klon) REAL cenv_vol(klon) REAL qt_env(klon), thetal_env(klon) REAL sqrtpi,sqrt2,sqrt2pi REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth REAL zdelta,qsatbef,zcor REAL Tbefth(klon), Tbefenv(klon) REAL qlbef REAL dqsatenv(klon), dqsatth(klon) REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) REAL qincloud(klon) REAL alenvl, aenvl REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc !------------------------------------------------------------------------------ ! Initialisation !------------------------------------------------------------------------------ sqrt2pi=sqrt(2.*rpi) sqrt2=sqrt(2.) sqrtpi=sqrt(rpi) !------------------------------------------------------------------------------- ! Thermal fraction calculation and standard deviation of the distribution !------------------------------------------------------------------------------- ! initialisations and calculation of temperature, humidity and saturation specific humidity cloudth_senv(:) = 0. cloudth_sth(:) = 0. cloudth_sigmaenv(:) = 0. cloudth_sigmath(:) = 0. DO ind1=1,klon Tbefenv(ind1) = temp(ind1) thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1) Tbefth(ind1) = thetal_th(ind1)*zpspsk(ind1) qt_env(ind1) = (qt(ind1)-frac_th(ind1)*qt_th(ind1))/(1.-frac_th(ind1)) !qt = a*qtth + (1-a)*qtenv ENDDO DO ind1=1,klon IF (frac_th(ind1).GT.min_frac_th_cld) THEN !Thermal and environnement ! Environment: alenv=(RD/RV*RLVTT*qsenv(ind1))/(rd*thetal_env(ind1)**2) aenv=1./(1.+(alenv*RLVTT/rcpd)) senv=aenv*(qt_env(ind1)-qsenv(ind1)) ! Thermals: alth=(RD/RV*RLVTT*qsth(ind1))/(rd*thetal_th(ind1)**2) ath=1./(1.+(alth*RLVTT/rcpd)) sth=ath*(qt_th(ind1)-qsth(ind1)) ! Standard deviation of the distributions ! environment sigma1s_fraca = (sigma1s_factor**0.5)*(frac_th(ind1)**sigma1s_power) / & & (1-frac_th(ind1))*((sth-senv)**2)**0.5 IF (cloudth_ratqsmin>0.) THEN sigma1s_ratqs = cloudth_ratqsmin*qt(ind1) ELSE sigma1s_ratqs = ratqs(ind1)*qt(ind1) ENDIF sigma1s = sigma1s_fraca + sigma1s_ratqs IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then sigma1s = ratqs(ind1)*qt(ind1)*aenv ENDIF ! thermals sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((frac_th(ind1)+0.02)**sigma2s_power))+0.002*qt_th(ind1) IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1).ne.0) then sigma2s = sigma_qtherm(ind1)*ath ENDIF ! surface cloud fraction deltasenv=aenv*vert_alpha*sigma1s deltasth=ath*vert_alpha_th*sigma2s xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) exp_xenv1 = exp(-1.*xenv1**2) exp_xenv2 = exp(-1.*xenv2**2) xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) exp_xth1 = exp(-1.*xth1**2) exp_xth2 = exp(-1.*xth2**2) cth(ind1)=0.5*(1.-1.*erf(xth1)) cenv(ind1)=0.5*(1.-1.*erf(xenv1)) ctot(ind1)=frac_th(ind1)*cth(ind1)+(1.-1.*frac_th(ind1))*cenv(ind1) ctotth(ind1)=frac_th(ind1)*cth(ind1) !volume cloud fraction and condensed water !environnement IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 IntJ_CF=0.5*(1.-1.*erf(xenv2)) IF (deltasenv .LT. 1.e-10) THEN qcenv(ind1)=IntJ cenv_vol(ind1)=IntJ_CF ELSE IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) qcenv(ind1)=IntJ+IntI1+IntI2+IntI3 cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF ENDIF !thermals IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 IntJ_CF=0.5*(1.-1.*erf(xth2)) IF (deltasth .LT. 1.e-10) THEN qcth(ind1)=IntJ cth_vol(ind1)=IntJ_CF ELSE IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) qcth(ind1)=IntJ+IntI1+IntI2+IntI3 cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF ENDIF ! total qctot(ind1)=frac_th(ind1)*qcth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1) ctot_vol(ind1)=frac_th(ind1)*cth_vol(ind1)+(1.-1.*frac_th(ind1))*cenv_vol(ind1) IF (cenv(ind1).LT.min_neb_th.and.cth(ind1).LT.min_neb_th) THEN ctot(ind1)=0. ctot_vol(ind1)=0. qcloud(ind1)=qsenv(ind1) qincloud(ind1)=0. ELSE qincloud(ind1)=qctot(ind1)/ctot(ind1) !to prevent situations with cloud condensed water greater than available total water qincloud(ind1)=min(qincloud(ind1),qt(ind1)/ctot(ind1)) ! we assume that water vapor in cloud is qsenv qcloud(ind1)=qincloud(ind1)+qsenv(ind1) ENDIF ! Outputs used to check the PDFs cloudth_senv(ind1) = senv cloudth_sth(ind1) = sth cloudth_sigmaenv(ind1) = sigma1s cloudth_sigmath(ind1) = sigma2s ENDIF ! selection of grid points concerned by thermals ENDDO !loop on klon RETURN END SUBROUTINE condensation_cloudth !***************************************************************************************** !***************************************************************************************** ! pre-cmip7 routines are below and are becoming obsolete !***************************************************************************************** !***************************************************************************************** SUBROUTINE cloudth(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) use lmdz_lscp_ini, only: iflag_cloudth_vert,iflag_ratqs USE yomcst_mod_h USE yoethf_mod_h IMPLICIT NONE !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== INCLUDE "FCTTRE.h" INTEGER itap,ind1,ind2 INTEGER ngrid,klev,klon,l,ig real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv REAL ztv(ngrid,klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid,klev) REAL fraca(ngrid,klev+1) REAL zpspsk(ngrid,klev) REAL paprs(ngrid,klev+1) REAL pplay(ngrid,klev) REAL ztla(ngrid,klev) REAL zthl(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi REAL rdd,cppd,Lv REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,xth,xenv REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Gestion de deux versions de cloudth !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (iflag_cloudth_vert.GE.1) THEN CALL cloudth_vert(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t) RETURN ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------------------------- ! Initialisation des variables r?elles !------------------------------------------------------------------------------- sigma1(:,ind2)=0. sigma2(:,ind2)=0. qlth(:,ind2)=0. qlenv(:,ind2)=0. qltot(:,ind2)=0. rneb(:,ind2)=0. qcloud(:)=0. cth(:,ind2)=0. cenv(:,ind2)=0. ctot(:,ind2)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ?cart-types des distributions !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) ! zqenv(ind1)=po(ind1) Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) ath=1./(1.+(alth*Lv/cppd)) sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !------------------------------------------------------------------------------ ! Calcul des ?cart-types pour s !------------------------------------------------------------------------------ ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) then ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 !------------------------------------------------------------------------------ ! Calcul de l'eau condens?e et de la couverture nuageuse !------------------------------------------------------------------------------ sqrt2pi=sqrt(2.*pi) xth=sth/(sqrt(2.)*sigma2s) xenv=senv/(sqrt(2.)*sigma1s) cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (ctot(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) endif else ! gaussienne environnement seule zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sigma1s=ratqs(ind1,ind2)*zqenv(ind1) sqrt2pi=sqrt(2.*pi) xenv=senv/(sqrt(2.)*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif enddo return ! end END SUBROUTINE cloudth !=========================================================================== SUBROUTINE cloudth_vert(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t) !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== USE yoethf_mod_h use lmdz_lscp_ini, only: iflag_cloudth_vert, vert_alpha USE yomcst_mod_h IMPLICIT NONE INCLUDE "FCTTRE.h" INTEGER itap,ind1,ind2 INTEGER ngrid,klev,klon,l,ig REAL ztv(ngrid,klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid,klev) REAL fraca(ngrid,klev+1) REAL zpspsk(ngrid,klev) REAL paprs(ngrid,klev+1) REAL pplay(ngrid,klev) REAL ztla(ngrid,klev) REAL zthl(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi REAL rdd,cppd,Lv,sqrt2,sqrtpi REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,xth,xenv REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur ! Change the width of the PDF used for vertical subgrid scale heterogeneity ! (J Jouhaud, JL Dufresne, JB Madeleine) REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) !------------------------------------------------------------------------------ ! Initialisation des variables r?elles !------------------------------------------------------------------------------ sigma1(:,ind2)=0. sigma2(:,ind2)=0. qlth(:,ind2)=0. qlenv(:,ind2)=0. qltot(:,ind2)=0. rneb(:,ind2)=0. qcloud(:)=0. cth(:,ind2)=0. cenv(:,ind2)=0. ctot(:,ind2)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) sqrt2=sqrt(2.) sqrtpi=sqrt(pi) !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ?cart-types des distributions !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) ! zqenv(ind1)=po(ind1) Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) ath=1./(1.+(alth*Lv/cppd)) sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !------------------------------------------------------------------------------ ! Calcul des ?cart-types pour s !------------------------------------------------------------------------------ sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) then ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 !------------------------------------------------------------------------------ ! Calcul de l'eau condens?e et de la couverture nuageuse !------------------------------------------------------------------------------ sqrt2pi=sqrt(2.*pi) xth=sth/(sqrt(2.)*sigma2s) xenv=senv/(sqrt(2.)*sigma1s) cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) IF (iflag_cloudth_vert == 1) THEN !------------------------------------------------------------------------------- ! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) ! deltasenv=aenv*0.01*po(ind1) ! deltasth=ath*0.01*zqta(ind1,ind2) xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s) xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s) xth1=(sth-deltasth)/(sqrt(2.)*sigma2s) xth2=(sth+deltasth)/(sqrt(2.)*sigma2s) coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) cth(ind1,ind2)=0.5*(1.+1.*erf(xth2)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1)) IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlenv(ind1,ind2)=IntJ ! print*, qlenv(ind1,ind2),'VERIF EAU' IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1)) ! IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2)) ! IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1)) IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2)) IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1)) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlth(ind1,ind2)=IntJ ! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ELSE IF (iflag_cloudth_vert == 2) THEN !------------------------------------------------------------------------------- ! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) ! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) ! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) deltasenv=aenv*vert_alpha*sigma1s deltasth=ath*vert_alpha*sigma2s xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) ! coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) ! coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2) IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2)) ! IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) ! IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) ! IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlenv(ind1,ind2)=IntJ ! print*, qlenv(ind1,ind2),'VERIF EAU' IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2) IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2)) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlth(ind1,ind2)=IntJ ! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ENDIF ! of if (iflag_cloudth_vert==1 or 2) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & ! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) endif ! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif' else ! gaussienne environnement seule zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sigma1s=ratqs(ind1,ind2)*zqenv(ind1) sqrt2pi=sqrt(2.*pi) xenv=senv/(sqrt(2.)*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif enddo return ! end END SUBROUTINE cloudth_vert SUBROUTINE cloudth_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,sigma_qtherm,zqs,t, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) use lmdz_lscp_ini, only: iflag_cloudth_vert USE yomcst_mod_h USE yoethf_mod_h IMPLICIT NONE !=========================================================================== ! Author : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== INCLUDE "FCTTRE.h" integer, intent(in) :: ind2 integer, intent(in) :: ngrid,klev real, dimension(ngrid,klev), intent(in) :: ztv real, dimension(ngrid), intent(in) :: po real, dimension(ngrid,klev), intent(in) :: zqta real, dimension(ngrid,klev+1), intent(in) :: fraca real, dimension(ngrid), intent(out) :: qcloud real, dimension(ngrid,klev), intent(out) :: ctot real, dimension(ngrid,klev), intent(out) :: ctot_vol real, dimension(ngrid,klev), intent(in) :: zpspsk real, dimension(ngrid,klev+1), intent(in) :: paprs real, dimension(ngrid,klev), intent(in) :: pplay real, dimension(ngrid,klev), intent(in) :: ztla real, dimension(ngrid,klev), intent(inout) :: zthl real, dimension(ngrid,klev), intent(in) :: ratqs,sigma_qtherm real, dimension(ngrid), intent(in) :: zqs real, dimension(ngrid,klev), intent(in) :: t real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv REAL zqenv(ngrid) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL cth_vol(ngrid,klev) REAL cenv_vol(ngrid,klev) REAL rneb(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi REAL rdd,cppd,Lv REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2 REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) INTEGER :: ind1,l, ig IF (iflag_cloudth_vert.GE.1) THEN CALL cloudth_vert_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,sigma_qtherm,zqs,t, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) RETURN ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------------------------- ! Initialisation des variables r?elles !------------------------------------------------------------------------------- sigma1(:,ind2)=0. sigma2(:,ind2)=0. qlth(:,ind2)=0. qlenv(:,ind2)=0. qltot(:,ind2)=0. rneb(:,ind2)=0. qcloud(:)=0. cth(:,ind2)=0. cenv(:,ind2)=0. ctot(:,ind2)=0. cth_vol(:,ind2)=0. cenv_vol(:,ind2)=0. ctot_vol(:,ind2)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) sqrt2=sqrt(2.) sqrtpi=sqrt(pi) !------------------------------------------------------------------------------- ! Cloud fraction in the thermals and standard deviation of the PDFs !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 !po = qt de l'environnement ET des thermique !zqenv = qt environnement !zqsatenv = qsat environnement !zthl = Tl environnement Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 ath=1./(1.+(alth*Lv/cppd)) !al, p84 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 !zqta = qt thermals !zqsatth = qsat thermals !ztla = Tl thermals !------------------------------------------------------------------------------ ! s standard deviations !------------------------------------------------------------------------------ ! tests ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1) ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2) ! final option sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) !------------------------------------------------------------------------------ ! Condensed water and cloud cover !------------------------------------------------------------------------------ xth=sth/(sqrt2*sigma2s) xenv=senv/(sqrt2*sigma1s) cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) ctot_vol(ind1,ind2)=ctot(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) if (ctot(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) endif else ! Environnement only, follow the if l.110 zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sigma1s=ratqs(ind1,ind2)*zqenv(ind1) xenv=senv/(sqrt2*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=ctot(ind1,ind2) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183 enddo ! from the loop on ngrid l.108 return ! end END SUBROUTINE cloudth_v3 !=========================================================================== SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,sigma_qtherm,zqs,t, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== use yoethf_mod_h use lmdz_lscp_ini, only : iflag_cloudth_vert,iflag_ratqs use lmdz_lscp_ini, only : vert_alpha,vert_alpha_th, sigma1s_factor, sigma1s_power , sigma2s_factor , sigma2s_power , cloudth_ratqsmin , iflag_cloudth_vert_noratqs USE yomcst_mod_h IMPLICIT NONE INCLUDE "FCTTRE.h" INTEGER itap,ind1,ind2 INTEGER ngrid,klev,klon,l,ig real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv REAL ztv(ngrid,klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid,klev) REAL fraca(ngrid,klev+1) REAL zpspsk(ngrid,klev) REAL paprs(ngrid,klev+1) REAL pplay(ngrid,klev) REAL ztla(ngrid,klev) REAL zthl(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) REAL cth_vol(ngrid,klev) REAL cenv_vol(ngrid,klev) REAL ctot_vol(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi REAL rdd,cppd,Lv REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL ratqs(ngrid,klev),sigma_qtherm(ngrid,klev) ! determine la largeur de distribution de vapeur ! Change the width of the PDF used for vertical subgrid scale heterogeneity ! (J Jouhaud, JL Dufresne, JB Madeleine) REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) REAL rhodz(ngrid,klev) REAL zrho(ngrid,klev) REAL dz(ngrid,klev) DO ind1 = 1, ngrid !Layer calculation rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2 zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3 dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre END DO !------------------------------------------------------------------------------ ! Initialize !------------------------------------------------------------------------------ sigma1(:,ind2)=0. sigma2(:,ind2)=0. qlth(:,ind2)=0. qlenv(:,ind2)=0. qltot(:,ind2)=0. rneb(:,ind2)=0. qcloud(:)=0. cth(:,ind2)=0. cenv(:,ind2)=0. ctot(:,ind2)=0. cth_vol(:,ind2)=0. cenv_vol(:,ind2)=0. ctot_vol(:,ind2)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) sqrt2=sqrt(2.) sqrtpi=sqrt(pi) !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ecart-types des distributions !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 !zqenv = qt environnement !zqsatenv = qsat environnement !zthl = Tl environnement Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 ath=1./(1.+(alth*Lv/cppd)) !al, p84 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 !zqta = qt thermals !zqsatth = qsat thermals !ztla = Tl thermals !------------------------------------------------------------------------------ ! s standard deviation !------------------------------------------------------------------------------ sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / & & (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 ! sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 IF (cloudth_ratqsmin>0.) THEN sigma1s_ratqs = cloudth_ratqsmin*po(ind1) ELSE sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1) ENDIF sigma1s = sigma1s_fraca + sigma1s_ratqs sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1,ind2).ne.0) then sigma2s = sigma_qtherm(ind1,ind2)*ath ENDIF ENDIF ! tests ! sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1) ! sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) then ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 IF (iflag_cloudth_vert == 1) THEN !------------------------------------------------------------------------------- ! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs !------------------------------------------------------------------------------- deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s) xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s) xth1=(sth-deltasth)/(sqrt(2.)*sigma2s) xth2=(sth+deltasth)/(sqrt(2.)*sigma2s) coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) cth(ind1,ind2)=0.5*(1.+1.*erf(xth2)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) ! Environment IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1)) IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! Thermal IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1)) IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2)) IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1)) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ELSE IF (iflag_cloudth_vert >= 3) THEN IF (iflag_cloudth_vert < 5) THEN !------------------------------------------------------------------------------- ! Version 3: Changes by J. Jouhaud; condensation for q > -delta s !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) ! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) ! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) IF (iflag_cloudth_vert == 3) THEN deltasenv=aenv*vert_alpha*sigma1s deltasth=ath*vert_alpha_th*sigma2s ELSE IF (iflag_cloudth_vert == 4) THEN IF (iflag_cloudth_vert_noratqs == 1) THEN deltasenv=vert_alpha*max(sigma1s_fraca,1e-10) deltasth=vert_alpha_th*sigma2s ELSE deltasenv=vert_alpha*sigma1s deltasth=vert_alpha_th*sigma2s ENDIF ENDIF xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) exp_xenv1 = exp(-1.*xenv1**2) exp_xenv2 = exp(-1.*xenv2**2) xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) exp_xth1 = exp(-1.*xth1**2) exp_xth2 = exp(-1.*xth2**2) !CF_surfacique cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) !CF_volumique & eau condense !environnement IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 IntJ_CF=0.5*(1.-1.*erf(xenv2)) if (deltasenv .lt. 1.e-10) then qlenv(ind1,ind2)=IntJ cenv_vol(ind1,ind2)=IntJ_CF else IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF endif !thermique IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 IntJ_CF=0.5*(1.-1.*erf(xth2)) if (deltasth .lt. 1.e-10) then qlth(ind1,ind2)=IntJ cth_vol(ind1,ind2)=IntJ_CF else IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF endif qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) ELSE IF (iflag_cloudth_vert == 5) THEN sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) & /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) & +ratqs(ind1,ind2)*po(ind1) !Environment sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) !Thermals !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) xth=sth/(sqrt(2.)*sigma2s) xenv=senv/(sqrt(2.)*sigma1s) !Volumique cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) !print *,'jeanjean_CV=',ctot_vol(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) !Surfacique !Neggers !beta=0.0044 !inverse_rho=1.+beta*dz(ind1,ind2) !print *,'jeanjean : beta=',beta !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) !Brooks a_Brooks=0.6694 b_Brooks=0.1882 A_Maj_Brooks=0.1635 !-- sans shear !A_Maj_Brooks=0.17 !-- ARM LES !A_Maj_Brooks=0.18 !-- RICO LES !A_Maj_Brooks=0.19 !-- BOMEX LES Dx_Brooks=200000. f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) !print *,'jeanjean_f=',f_Brooks cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.)) cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.)) ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) !print *,'JJ_ctot_1',ctot(ind1,ind2) ENDIF ! of if (iflag_cloudth_vert<5) ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4) ! if (ctot(ind1,ind2).lt.1.e-10) then if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. ctot_vol(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & ! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) endif else ! gaussienne environnement seule zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sth=0. sigma1s=ratqs(ind1,ind2)*zqenv(ind1) sigma2s=0. sqrt2pi=sqrt(2.*pi) xenv=senv/(sqrt(2.)*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=ctot(ind1,ind2) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ! ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492 ! Outputs used to check the PDFs cloudth_senv(ind1,ind2) = senv cloudth_sth(ind1,ind2) = sth cloudth_sigmaenv(ind1,ind2) = sigma1s cloudth_sigmath(ind1,ind2) = sigma2s enddo ! from the loop on ngrid l.333 return ! end END SUBROUTINE cloudth_vert_v3 ! SUBROUTINE cloudth_v6(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,T, & & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) USE yoethf_mod_h USE lmdz_lscp_ini, only: iflag_cloudth_vert USE yomcst_mod_h IMPLICIT NONE INCLUDE "FCTTRE.h" !Domain variables INTEGER ngrid !indice Max lat-lon INTEGER klev !indice Max alt real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv INTEGER ind1 !indice in [1:ngrid] INTEGER ind2 !indice in [1:klev] !thermal plume fraction REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox !temperatures REAL T(ngrid,klev) !temperature REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables) REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90) REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th) REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env) !pressure REAL paprs(ngrid,klev+1) !pressure at the interface of levels REAL pplay(ngrid,klev) !pressure at the middle of the level !humidity REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution REAL po(ngrid) !total water (qt) REAL zqenv(ngrid) !total water in the environment (qt_env) REAL zqta(ngrid,klev) !total water in the thermals (qt_th) REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th) REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env) REAL qlth(ngrid,klev) !condensed water in the thermals REAL qlenv(ngrid,klev) !condensed water in the environment REAL qltot(ngrid,klev) !condensed water in the gridbox !cloud fractions REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox !PDF of saturation deficit variables REAL rdd,cppd,Lv REAL Tbef,zdelta,qsatbef,zcor REAL alth,alenv,ath,aenv REAL sth,senv !saturation deficits in the thermals and environment REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF !cloud fraction variables REAL xth,xenv REAL inverse_rho,beta !Neggers et al. (2011) method REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method !Incloud total water variables REAL zqs(ngrid) !q_sat REAL qcloud(ngrid) !eau totale dans le nuage !Some arithmetic variables REAL pi,sqrt2,sqrt2pi !Depth of the layer REAL dz(ngrid,klev) !epaisseur de la couche en metre REAL rhodz(ngrid,klev) REAL zrho(ngrid,klev) DO ind1 = 1, ngrid rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2] zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3] dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m] END DO !------------------------------------------------------------------------------ ! Initialization !------------------------------------------------------------------------------ qlth(:,ind2)=0. qlenv(:,ind2)=0. qltot(:,ind2)=0. cth_vol(:,ind2)=0. cenv_vol(:,ind2)=0. ctot_vol(:,ind2)=0. cth_surf(:,ind2)=0. cenv_surf(:,ind2)=0. ctot_surf(:,ind2)=0. qcloud(:)=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2=sqrt(2.) sqrt2pi=sqrt(2.*pi) DO ind1=1,ngrid !------------------------------------------------------------------------------- !Both thermal and environment in the gridbox !------------------------------------------------------------------------------- IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN !-------------------------------------------- !calcul de qsat_env !-------------------------------------------- Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef !-------------------------------------------- !calcul de s_env !-------------------------------------------- alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 these Arnaud Jam aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul de qsat_th !-------------------------------------------- Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef !-------------------------------------------- !calcul de s_th !-------------------------------------------- alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 these Arnaud Jam ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul standard deviations bi-Gaussian PDF !-------------------------------------------- sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) & /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) & +ratqs(ind1,ind2)*po(ind1) xth=sth/(sqrt2*sigma_th) xenv=senv/(sqrt2*sigma_env) !-------------------------------------------- !Cloud fraction by volume CF_vol !-------------------------------------------- cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) !-------------------------------------------- !Condensed water qc !-------------------------------------------- qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2)) qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) !-------------------------------------------- !Cloud fraction by surface CF_surf !-------------------------------------------- !Method Neggers et al. (2011) : ok for cumulus clouds only !beta=0.0044 (Jouhaud et al.2018) !inverse_rho=1.+beta*dz(ind1,ind2) !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho !Method Brooks et al. (2005) : ok for all types of clouds a_Brooks=0.6694 b_Brooks=0.1882 A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) !-------------------------------------------- !Incloud Condensed water qcloud !-------------------------------------------- if (ctot_surf(ind1,ind2) .lt. 1.e-10) then ctot_vol(ind1,ind2)=0. ctot_surf(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1) endif !------------------------------------------------------------------------------- !Environment only in the gridbox !------------------------------------------------------------------------------- ELSE !-------------------------------------------- !calcul de qsat_env !-------------------------------------------- Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef !-------------------------------------------- !calcul de s_env !-------------------------------------------- alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 these Arnaud Jam aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul standard deviations Gaussian PDF !-------------------------------------------- zqenv(ind1)=po(ind1) sigma_env=ratqs(ind1,ind2)*zqenv(ind1) xenv=senv/(sqrt2*sigma_env) !-------------------------------------------- !Cloud fraction by volume CF_vol !-------------------------------------------- ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) !-------------------------------------------- !Condensed water qc !-------------------------------------------- qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2)) !-------------------------------------------- !Cloud fraction by surface CF_surf !-------------------------------------------- !Method Neggers et al. (2011) : ok for cumulus clouds only !beta=0.0044 (Jouhaud et al.2018) !inverse_rho=1.+beta*dz(ind1,ind2) !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho !Method Brooks et al. (2005) : ok for all types of clouds a_Brooks=0.6694 b_Brooks=0.1882 A_Maj_Brooks=0.1635 !-- sans dependence au shear Dx_Brooks=200000. f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) !-------------------------------------------- !Incloud Condensed water qcloud !-------------------------------------------- if (ctot_surf(ind1,ind2) .lt. 1.e-8) then ctot_vol(ind1,ind2)=0. ctot_surf(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2) endif END IF ! From the separation (thermal/envrionnement) et (environnement only) ! Outputs used to check the PDFs cloudth_senv(ind1,ind2) = senv cloudth_sth(ind1,ind2) = sth cloudth_sigmaenv(ind1,ind2) = sigma_env cloudth_sigmath(ind1,ind2) = sigma_th END DO ! From the loop on ngrid return END SUBROUTINE cloudth_v6 END MODULE lmdz_lscp_condensation