!$gpum horizontal klon ngrid MODULE lmdz_lscp_main IMPLICIT NONE CONTAINS !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE lscp(klon, klev, dtime, missing_val, & paprs, pplay, omega, temp, qt, ql_seri, qi_seri, & ratqs, sigma_qtherm, ptconv, cfcon_old, qvcon_old, & qccon_old, cfcon, qvcon, qccon, & d_t, d_q, d_ql, d_qi, rneb, rneblsvol, & pfraclr, pfracld, & cldfraliq, cldfraliqth, & sigma2_icefracturb,sigma2_icefracturbth, & mean_icefracturb,mean_icefracturbth, & radocond, radicefrac, rain, snow, & frac_impa, frac_nucl, beta, & prfl, psfl, rhcl, qta, fraca, & tv, pspsk, tla, thl, wth, iflag_cld_th, & iflag_ice_thermo, distcltop, temp_cltop, & tke, tke_dissip, & entr_therm, detr_therm, & cell_area, stratomask, & cf_seri, qvc_seri, u_seri, v_seri, & qsub, qissr, qcld, subfra, issrfra, gamma_cond, & dcf_sub, dcf_con, dcf_mix, & dqi_sed, dcf_sed, dqvc_sed, & dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, & dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati, & cfc_seri, qtc_seri, qic_seri, nic_seri, & qice_cont, flight_dist, flight_fuel, & contfrarad, qradice_cont, & Tcritcont, qcritcont, potcontfraP, potcontfraNP, & cloudth_sth, & cloudth_senv, cloudth_sigmath, cloudth_sigmaenv, & qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, & dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,& dqsmelt, dqsfreez) !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD) ! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD) !-------------------------------------------------------------------------------- ! Date: 01/2021 !-------------------------------------------------------------------------------- ! Aim: Large Scale Clouds and Precipitation (LSCP) ! ! This code is a new version of the fisrtilp.F90 routine, which is itself a ! merge of 'first' (superrsaturation physics, P. LeVan K. Laval) ! and 'ilp' (il pleut, L. Li) ! ! Compared to the original fisrtilp code, lscp ! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.) ! -> consider always precipitation thermalisation (fl_cor_ebil>0) ! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T) ! -> option "oldbug" by JYG has been removed ! -> iflag_t_glace >0 always ! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always) ! -> rectangular distribution from L. Li no longer available ! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt) !-------------------------------------------------------------------------------- ! References: ! ! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2 ! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y ! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3 ! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379 ! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046 ! - Touzze-Peifert Ludo, PhD thesis, p117-124 ! ------------------------------------------------------------------------------- ! Code structure: ! ! P0> Thermalization of the precipitation coming from the overlying layer ! P1> Evaporation of the precipitation (falling from the k+1 level) ! P2> Cloud formation (at the k level) ! P2.A.1> With the PDFs, calculation of cloud properties using the inital ! values of T and Q ! P2.A.2> Coupling between condensed water and temperature ! P2.A.3> Calculation of final quantities associated with cloud formation ! P2.B> Release of Latent heat after cloud formation ! P3> Autoconversion to precipitation (k-level) ! P4> Wet scavenging !------------------------------------------------------------------------------ ! Some preliminary comments (JBM) : ! ! The cloud water that the radiation scheme sees is not the same that the cloud ! water used in the physics and the dynamics ! ! During the autoconversion to precipitation (P3 step), radocond (cloud water used ! by the radiation scheme) is calculated as an average of the water that remains ! in the cloud during the precipitation and not the water remaining at the end ! of the time step. The latter is used in the rest of the physics and advected ! by the dynamics. ! ! In summary: ! ! Radiation: ! xflwc(newmicro)+xfiwc(newmicro) = ! radocond=lwcon(nc)+iwcon(nc) ! ! Notetheless, be aware of: ! ! radocond .NE. ocond(nc) ! i.e.: ! lwcon(nc)+iwcon(nc) .NE. ocond(nc) ! but oliq+(ocond-oliq) .EQ. ocond ! (which is not trivial) !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! USE de modules contenant des fonctions. USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat USE lmdz_lscp_condensation, ONLY : cloudth, cloudth_v3, cloudth_v6, condensation_cloudth USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld ! Use du module lmdz_lscp_ini contenant les constantes USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0 USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat USE lmdz_lscp_ini, ONLY : min_frac_th_cld, temp_nowater, rho_ice USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG, RPI USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac USE lmdz_lscp_ini, ONLY : ok_weibull_warm_clouds, ok_no_issr_strato, ok_ice_sedim USE lmdz_lscp_ini, ONLY : ok_plane_contrail, ok_unadjusted_contrails USE lmdz_lscp_ini, ONLY : ok_nodeep_lscp, ok_higher_cirrus_cover USE lmdz_lscp_ini, ONLY : ok_lscp_mergecond, gamma_mixth USE lmdz_lscp_ini, ONLY : eff2vol_radius_contrails ! Temporary call for Lamquin et al (2012) diagnostics USE phys_local_var_mod, ONLY : issrfra100to150, issrfra150to200, issrfra200to250 USE phys_local_var_mod, ONLY : issrfra250to300, issrfra300to400, issrfra400to500 USE phys_local_var_mod, ONLY : AEI_contrails, AEI_surv_contrails USE phys_local_var_mod, ONLY : fsurv_contrails, section_contrails USE phys_local_var_mod, ONLY : dcfc_ini, dqic_ini, dqtc_ini, dnic_ini USE phys_local_var_mod, ONLY : dcfc_sub, dqic_sub, dqtc_sub, dnic_sub USE phys_local_var_mod, ONLY : dcfc_mix, dqic_mix, dqtc_mix, dnic_mix USE phys_local_var_mod, ONLY : dnic_agg, dqic_adj, dqtc_adj USE phys_local_var_mod, ONLY : dcfc_sed, dqic_sed, dqtc_sed, dnic_sed USE phys_local_var_mod, ONLY : dcfc_auto, dqic_auto, dqtc_auto, dnic_auto USE phys_local_var_mod, ONLY : dcf_auto, dqi_auto, dqvc_auto USE phys_local_var_mod, ONLY : nice_ygcont, iwc_ygcont, rvol_ygcont, tau_ygcont USE phys_local_var_mod, ONLY : nice_vscont, iwc_vscont, rvol_vscont, tau_vscont USE phys_local_var_mod, ONLY : nice_cont, iwc_cont, rvol_cont, tau_cont IMPLICIT NONE !=============================================================================== ! VARIABLES DECLARATION !=============================================================================== ! INPUT VARIABLES: !----------------- INTEGER, INTENT(IN) :: klon,klev ! number of horizontal grid points and vertical levels INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics REAL, INTENT(IN) :: dtime ! time step [s] REAL, INTENT(IN) :: missing_val ! missing value for output REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa] REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa] REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K) REAL, DIMENSION(klon,klev), INTENT(IN) :: omega ! vertical velocity [Pa/s] REAL, DIMENSION(klon,klev), INTENT(IN) :: qt ! total specific humidity (in vapor phase in input) [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: ql_seri ! liquid specific content [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: qi_seri ! ice specific content [kg/kg] REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke ! turbulent kinetic energy [m2/s2] REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke_dissip ! TKE dissipation [m2/s3] REAL, DIMENSION(klon,klev), INTENT(IN) :: entr_therm ! thermal plume entrainment rate * dz [kg/s/m2] REAL, DIMENSION(klon,klev), INTENT(IN) :: detr_therm ! thermal plume detrainment rate * dz [kg/s/m2] LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active REAL, DIMENSION(klon,klev), INTENT(IN) :: cfcon_old ! cloud fraction from deep convection from previous timestep [-] REAL, DIMENSION(klon,klev), INTENT(INOUT):: qvcon_old ! in-cloud vapor specific humidity from deep convection from previous timestep [kg/kg] REAL, DIMENSION(klon,klev), INTENT(INOUT):: qccon_old ! in-cloud condensed specific humidity from deep convection from previous timestep [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: cfcon ! cloud fraction from deep convection [-] REAL, DIMENSION(klon,klev), INTENT(IN) :: qvcon ! in-cloud vapor specific humidity from deep convection [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: qccon ! in-cloud condensed specific humidity from deep convection [kg/kg] !Inputs associated with thermal plumes INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds REAL, DIMENSION(klon,klev), INTENT(IN) :: tv ! virtual potential temperature [K] REAL, DIMENSION(klon,klev), INTENT(IN) :: qta ! specific humidity within thermals [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-] REAL, DIMENSION(klon,klev), INTENT(IN) :: pspsk ! exner potential (p/100000)**(R/cp) REAL, DIMENSION(klon,klev), INTENT(IN) :: tla ! liquid potential temperature within thermals [K] REAL, DIMENSION(klon,klev), INTENT(IN) :: wth ! vertical velocity within thermals [m/s] REAL, DIMENSION(klon,klev), INTENT(IN) :: sigma_qtherm ! controls the saturation deficit distribution width in thermals [-] ! INPUT/OUTPUT variables !------------------------ REAL, DIMENSION(klon,klev), INTENT(INOUT) :: thl ! liquid potential temperature [K] REAL, DIMENSION(klon,klev), INTENT(INOUT) :: ratqs ! function that sets the large-scale water distribution width ! INPUT/OUTPUT condensation and ice supersaturation !-------------------------------------------------- REAL, DIMENSION(klon,klev), INTENT(INOUT):: cf_seri ! cloud fraction [-] REAL, DIMENSION(klon,klev), INTENT(INOUT):: qvc_seri ! cloudy water vapor [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: u_seri ! eastward wind [m/s] REAL, DIMENSION(klon,klev), INTENT(IN) :: v_seri ! northward wind [m/s] REAL, DIMENSION(klon), INTENT(IN) :: cell_area ! area of each cell [m2] REAL, DIMENSION(klon,klev), INTENT(IN) :: stratomask ! fraction of stratosphere (0 or 1) ! INPUT/OUTPUT aviation !-------------------------------------------------- REAL, DIMENSION(klon,klev), INTENT(INOUT):: cfc_seri ! contrail fraction [-] REAL, DIMENSION(klon,klev), INTENT(INOUT):: qtc_seri ! contrail total specific humidity [kg/kg] REAL, DIMENSION(klon,klev), INTENT(INOUT):: qic_seri ! contrail ice specific humidity [kg/kg] REAL, DIMENSION(klon,klev), INTENT(INOUT):: nic_seri ! contrail ice crystals concentration [#/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: flight_dist ! aviation distance flown concentration [m/s/m3] REAL, DIMENSION(klon,klev), INTENT(IN) :: flight_fuel ! aviation fuel consumption concentration [kg/s/m3] ! OUTPUT variables !----------------- REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K] REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneblsvol ! cloud fraction per unit volume [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: cldfraliq ! liquid fraction of cloud fraction [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: cldfraliqth ! liquid fraction of cloud fraction in thermals [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: sigma2_icefracturb ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: mean_icefracturb ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: sigma2_icefracturbth ! Variance of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: mean_icefracturbth ! Mean of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: radocond ! condensed water used in the radiation scheme [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-] REAL, DIMENSION(klon), INTENT(OUT) :: rain ! surface large-scale rainfall [kg/s/m2] REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2] REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2] REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2] REAL, DIMENSION(klon,klev), INTENT(OUT) :: distcltop ! distance to cloud top [m] REAL, DIMENSION(klon,klev), INTENT(OUT) :: temp_cltop ! temperature of cloud top [K] REAL, DIMENSION(klon,klev), INTENT(OUT) :: beta ! conversion rate of condensed water ! fraction of aerosol scavenging through impaction and nucleation (for on-line) REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-] ! for condensation and ice supersaturation REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsub !--specific total water content in sub-saturated clear sky region [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: qissr !--specific total water content in supersat region [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld !--specific total water content in cloudy region [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: subfra !--mesh fraction of subsaturated clear sky [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: issrfra !--mesh fraction of ISSR [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_cond !--coefficient governing the ice nucleation RHi threshold [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_sub !--cloud fraction tendency because of sublimation [s-1] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_con !--cloud fraction tendency because of condensation [s-1] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_mix !--cloud fraction tendency because of cloud mixing [s-1] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_adj !--specific ice content tendency because of temperature adjustment [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_sub !--specific ice content tendency because of sublimation [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_con !--specific ice content tendency because of condensation [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_mix !--specific ice content tendency because of cloud mixing [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_adj !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_sub !--specific cloud water vapor tendency because of sublimation [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_con !--specific cloud water vapor tendency because of condensation [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_mix !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_sed !--ice water content tendency due to sedmentation of ice crystals [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_sed !--cloud fraction tendency due to sedimentation of ice crystals [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_sed !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl !--saturation specific humidity wrt liquid [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsati !--saturation specific humidity wrt ice [kg/kg] ! for contrails and aviation REAL, DIMENSION(klon,klev), INTENT(OUT) :: qice_cont !--condensed water in contrails [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: contfrarad !--contrail fraction for radiation [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: qradice_cont !--condensed water in contrails used in the radiation scheme [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcritcont !--critical temperature for contrail formation [K] REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcritcont !--critical specific humidity for contrail formation [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: potcontfraP !--potential persistent contrail fraction [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: potcontfraNP !--potential non-persistent contrail fraction [-] ! for POPRECIP REAL, DIMENSION(klon,klev), INTENT(OUT) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqreva !--rain tendendy due to evaporation [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqssub !--snow tendency due to sublimation [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsrim !--snow tendency due to riming [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsmelt !--snow tendency due to melting [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrmelt !--rain tendency due to melting [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsfreez !--snow tendency due to freezing [kg/kg/s] REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrfreez !--rain tendency due to freezing [kg/kg/s] ! for thermals REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth !--mean saturation deficit in thermals REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv !--mean saturation deficit in environment REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath !--std of saturation deficit in thermals REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv !--std of saturation deficit in environment ! LOCAL VARIABLES: !---------------- REAL, DIMENSION(klon) :: qliq_in, qice_in, qvc_in, cldfra_in REAL, DIMENSION(klon,klev) :: ctot, rnebth, ctot_vol REAL, DIMENSION(klon,klev) :: wls !-- large scalce vertical velocity [m/s] REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi REAL, DIMENSION(klon) :: zqsth, zqslth, zqsith, zdqsth, zdqslth, zdqsith REAL :: zdelta REAL, DIMENSION(klon) :: zdqsdT_raw REAL, DIMENSION(klon) :: gammasat,dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T REAL, DIMENSION(klon) :: Tbef,Tbefth,Tbefthm1,qlibef,DT ! temperature, humidity and temp. variation during condensation iteration REAL :: num,denom REAL :: cste REAL, DIMENSION(klon) :: qincloud REAL, DIMENSION(klon) :: zrfl, zifl REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqth, zqn, zqnl REAL, DIMENSION(klon) :: zoliql, zoliqi REAL, DIMENSION(klon) :: zt, zp REAL, DIMENSION(klon) :: zfice, zficeth, zficeenv, zneb, zcf, zsnow REAL, DIMENSION(klon) :: dzfice, dzficeth, dzficeenv REAL, DIMENSION(klon) :: qtot, zeroklon ! Variables precipitation energy conservation REAL, DIMENSION(klon) :: zmqc REAL :: zalpha_tr REAL :: zfrac_lessi REAL, DIMENSION(klon) :: zprec_cond REAL, DIMENSION(klon) :: zlh_solid REAL, DIMENSION(klon) :: ztupnew REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl REAL, DIMENSION(klon) :: zrflclr, zrflcld REAL, DIMENSION(klon) :: ziflclr, ziflcld REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld REAL, DIMENSION(klon) :: tot_zneb REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop, zdeltaz REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl, zqliqth, zqiceth, zqvapclth, sursat_e, invtau_e ! for icefrac_lscp_turb ! for ice sedimentation REAL, DIMENSION(klon) :: dzsed, flsed, cfsed REAL, DIMENSION(klon) :: dzsed_abv, flsed_abv, cfsed_abv REAL, DIMENSION(klon) :: flauto REAL :: qice_sedim ! for quantity of condensates seen by radiation REAL, DIMENSION(klon) :: zradocond, zradoice REAL, DIMENSION(klon) :: zrho_up, zvelo_up REAL :: dz, coef_cover ! for condensation and ice supersaturation REAL, DIMENSION(klon) :: qvc, qvcl, shear REAL, DIMENSION(klon) :: zrneb_deep, zcond_deep REAL :: delta_z, deepconv_coef ! for contrails REAL, DIMENSION(klon) :: contfra, qcont, qvcont, Ncont REAL, DIMENSION(klon) :: totfra_in, qtot_in LOGICAL, DIMENSION(klon) :: pt_pron_clds REAL, DIMENSION(klon) :: dzsed_cont, flsed_cont, Nflsed_cont, cfsed_cont REAL, DIMENSION(klon) :: dzsed_cont_abv, flsed_cont_abv, Nflsed_cont_abv, cfsed_cont_abv REAL :: rho, rhodz, iwp_cont, rei_cont !--for Lamquin et al 2012 diagnostics REAL, DIMENSION(klon) :: issrfra100to150UP, issrfra150to200UP, issrfra200to250UP REAL, DIMENSION(klon) :: issrfra250to300UP, issrfra300to400UP, issrfra400to500UP INTEGER i, k, kk, iter INTEGER, DIMENSION(klon) :: n_i INTEGER ncoreczq LOGICAL iftop LOGICAL, DIMENSION(klon) :: stratiform_or_all_distrib,pticefracturb LOGICAL, DIMENSION(klon) :: keepgoing CHARACTER (len = 20) :: modname = 'lscp' CHARACTER (len = 80) :: abort_message !=============================================================================== ! INITIALISATION !=============================================================================== ! Few initial checks IF (iflag_fisrtilp_qsat .LT. 0) THEN abort_message = 'lscp cannot be used with iflag_fisrtilp<0' CALL abort_physic(modname,abort_message,1) ENDIF ! AA for 'safety' reasons zalpha_tr = 0. zfrac_lessi = 0. beta(:,:)= 0. ! Initialisation of variables: prfl(:,:) = 0.0 psfl(:,:) = 0.0 d_t(:,:) = 0.0 d_q(:,:) = 0.0 d_ql(:,:) = 0.0 d_qi(:,:) = 0.0 rneb(:,:) = 0.0 rnebth(:,:)=0.0 pfraclr(:,:)=0.0 pfracld(:,:)=0.0 cldfraliq(:,:)=0. sigma2_icefracturb(:,:)=0. mean_icefracturb(:,:)=0. cldfraliqth(:,:)=0. sigma2_icefracturbth(:,:)=0. mean_icefracturbth(:,:)=0. radocond(:,:) = 0.0 radicefrac(:,:) = 0.0 frac_nucl(:,:) = 1.0 frac_impa(:,:) = 1.0 rain(:) = 0.0 snow(:) = 0.0 zrfl(:) = 0.0 zifl(:) = 0.0 zneb(:) = seuil_neb zrflclr(:) = 0.0 ziflclr(:) = 0.0 zrflcld(:) = 0.0 ziflcld(:) = 0.0 tot_zneb(:) = 0.0 zeroklon(:) = 0.0 zdistcltop(:)=0.0 ztemp_cltop(:) = 0.0 ztupnew(:)=0.0 ctot_vol(:,:)=0.0 rneblsvol(:,:)=0.0 znebprecip(:)=0.0 znebprecipclr(:)=0.0 znebprecipcld(:)=0.0 distcltop(:,:)=0. temp_cltop(:,:)=0. !--Ice supersaturation gamma_cond(:,:) = 1. qissr(:,:) = 0. issrfra(:,:) = 0. dcf_sub(:,:) = 0. dcf_con(:,:) = 0. dcf_mix(:,:) = 0. dcf_sed(:,:) = 0. dcf_auto(:,:) = 0. dqi_adj(:,:) = 0. dqi_sub(:,:) = 0. dqi_con(:,:) = 0. dqi_mix(:,:) = 0. dqi_sed(:,:) = 0. dqi_auto(:,:) = 0. dqvc_adj(:,:) = 0. dqvc_sub(:,:) = 0. dqvc_con(:,:) = 0. dqvc_mix(:,:) = 0. dqvc_sed(:,:) = 0. dqvc_auto(:,:) = 0. qvc(:) = 0. shear(:) = 0. flsed(:) = 0. flauto(:) = 0. flsed_cont(:) = 0. Nflsed_cont(:) = 0. pt_pron_clds(:) = .FALSE. !--Contrails dcfc_ini(:,:) = 0. dqic_ini(:,:) = 0. dqtc_ini(:,:) = 0. dnic_ini(:,:) = 0. dcfc_sub(:,:) = 0. dqic_sub(:,:) = 0. dqtc_sub(:,:) = 0. dnic_sub(:,:) = 0. dqic_adj(:,:) = 0. dqtc_adj(:,:) = 0. dcfc_mix(:,:) = 0. dqic_mix(:,:) = 0. dqtc_mix(:,:) = 0. dnic_mix(:,:) = 0. dnic_agg(:,:) = 0. dcfc_sed(:,:) = 0. dqic_sed(:,:) = 0. dqtc_sed(:,:) = 0. dnic_sed(:,:) = 0. dcfc_auto(:,:) = 0. dqic_auto(:,:) = 0. dqtc_auto(:,:) = 0. dnic_auto(:,:) = 0. Tcritcont(:,:) = missing_val qcritcont(:,:) = missing_val potcontfraP(:,:) = missing_val potcontfraNP(:,:) = missing_val AEI_contrails(:,:) = missing_val AEI_surv_contrails(:,:) = missing_val fsurv_contrails(:,:) = missing_val section_contrails(:,:) = missing_val !--for Lamquin et al (2012) diagnostics issrfra100to150(:) = 0. issrfra100to150UP(:) = 0. issrfra150to200(:) = 0. issrfra150to200UP(:) = 0. issrfra200to250(:) = 0. issrfra200to250UP(:) = 0. issrfra250to300(:) = 0. issrfra250to300UP(:) = 0. issrfra300to400(:) = 0. issrfra300to400UP(:) = 0. issrfra400to500(:) = 0. issrfra400to500UP(:) = 0. !-- poprecip qraindiag(:,:)= 0. qsnowdiag(:,:)= 0. dqreva(:,:) = 0. dqrauto(:,:) = 0. dqrmelt(:,:) = 0. dqrfreez(:,:) = 0. dqrcol(:,:) = 0. dqssub(:,:) = 0. dqsauto(:,:) = 0. dqsrim(:,:) = 0. dqsagg(:,:) = 0. dqsfreez(:,:) = 0. dqsmelt(:,:) = 0. zqupnew(:) = 0. zqvapclr(:) = 0. !-- cloud phase useful variables wls(:,:)=-omega(:,:) / RG / (pplay(:,:)/RD/temp(:,:)) zqliq(:)=0. zqice(:)=0. zqvapcl(:)=0. zqliqth(:)=0. zqiceth(:)=0. zqvapclth(:)=0. sursat_e(:)=0. invtau_e(:)=0. pticefracturb(:)=.FALSE. !=============================================================================== ! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM !=============================================================================== ncoreczq=0 DO k = klev, 1, -1 IF (k.LE.klev-1) THEN iftop=.false. ELSE iftop=.true. ENDIF ! Initialisation temperature and specific humidity ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt ! at the end of the klon loop, a temperature incremtent d_t due to all processes ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated ! d_t = temperature tendency due to lscp ! The temperature of the overlying layer is updated here because needed for thermalization DO i = 1, klon zt(i)=temp(i,k) zq(i)=qt(i,k) zp(i)=pplay(i,k) qliq_in(i) = ql_seri(i,k) qice_in(i) = qi_seri(i,k) zcf(i) = 0. zfice(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation dzfice(i) = 0.0 zficeth(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation dzficeth(i) = 0.0 zficeenv(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation dzficeenv(i) = 0.0 IF (.not. iftop) THEN ztupnew(i) = temp(i,k+1) + d_t(i,k+1) zqupnew(i) = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1) !--zqs(i) is the saturation specific humidity in the layer above zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i)) ENDIF !c_iso init of iso ENDDO IF ( ok_ice_supersat ) THEN cldfra_in(:) = cf_seri(:,k) qvc_in(:) = qvc_seri(:,k) ENDIF ! -------------------------------------------------------------------- ! P1> Precipitation processes, before cloud formation: ! Thermalization of precipitation falling from the overlying layer AND ! Precipitation evaporation/sublimation/melting !--------------------------------------------------------------------- !================================================================ ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R) IF ( ok_poprecip ) THEN CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, & zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, & zqvapclr, zqupnew, flsed, flsed_cont, & cldfra_in, qvc_in, qliq_in, qice_in, & zrfl, zrflclr, zrflcld, & zifl, ziflclr, ziflcld, & dqreva(:,k), dqssub(:,k) & ) !================================================================ ELSE CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, & zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, & flsed, flsed_cont, & zrfl, zrflclr, zrflcld, & zifl, ziflclr, ziflcld, & dqreva(:,k), dqssub(:,k) & ) ENDIF ! (ok_poprecip) ! Calculation of qsat,L/cp*dqsat/dT and ncoreczq counter !------------------------------------------------------- qtot(:)=zq(:)+zmqc(:) CALL calc_qsat_ecmwf(klon,zt,qtot,zp,RTT,0,.false.,zqs,zdqs) DO i = 1, klon zdelta = MAX(0.,SIGN(1.,RTT-zt(i))) zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta) IF (zq(i) .LT. 1.e-15) THEN ncoreczq=ncoreczq+1 zq(i)=1.e-15 ENDIF ! c_iso: do something similar for isotopes ENDDO ! ------------------------------------------------------------------------- ! P2> Cloud formation including condensation and cloud phase determination !-------------------------------------------------------------------------- ! ! We always assume a 'fractional cloud' approach ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution ! is prescribed and depends on large scale variables and boundary layer ! properties) ! The decrease in condensed part due tu latent heating is taken into ! account ! ------------------------------------------------------------------- ! P2.1> With the PDFs (log-normal, bigaussian) ! cloud properties calculation with the initial values of t and q ! ---------------------------------------------------------------- ! initialise gammasat and stratiform_or_all_distrib stratiform_or_all_distrib(:)=.TRUE. gammasat(:)=1. ! part of the code that is supposed to become obsolete soon !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF (.NOT. ok_lscp_mergecond) THEN IF (iflag_cld_th.GE.5) THEN ! Cloud cover and content in meshes affected by shallow convection, ! are retrieved from a bi-gaussian distribution of the saturation deficit ! following Jam et al. 2013 IF (iflag_cloudth_vert.LE.2) THEN ! Old version of Arnaud Jam CALL cloudth(klon,klev,k,tv, & zq,qta,fraca, & qincloud,ctot,pspsk,paprs,pplay,tla,thl, & ratqs,zqs,temp, & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN ! Default version of Arnaud Jam CALL cloudth_v3(klon,klev,k,tv, & zq,qta,fraca, & qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, & ratqs,sigma_qtherm,zqs,temp, & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) ELSEIF (iflag_cloudth_vert.EQ.6) THEN ! Jean Jouhaud's version, with specific separation between surface and volume ! cloud fraction Decembre 2018 CALL cloudth_v6(klon,klev,k,tv, & zq,qta,fraca, & qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, & ratqs,zqs,temp, & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) ENDIF DO i=1,klon rneb(i,k)=ctot(i,k) ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.) rneblsvol(i,k)=ctot_vol(i,k) zqn(i)=qincloud(i) !--AB grid-mean vapor in the cloud - we assume saturation adjustment qvc(i) = rneb(i,k) * zqs(i) ENDDO ! Cloud phase final determination for clouds after "old" cloudth calls ! for those clouds, only temperature dependent phase partitioning (eventually modulated by ! distance to cloud top) is available ! compute distance to cloud top when cloud phase is computed depending on temperature ! and distance to cloud top IF (iflag_t_glace .GE. 4) THEN CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop) ENDIF CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice) ENDIF IF (iflag_cld_th .EQ. 5) THEN ! When iflag_cld_th=5, we always assume ! bi-gaussian distribution stratiform_or_all_distrib(:) = .FALSE. ELSEIF (iflag_cld_th .GE. 6) THEN ! stratiform distribution only when no thermals stratiform_or_all_distrib(:) = fraca(:,k) < min_frac_th_cld ENDIF ENDIF ! .not. ok_lscp_mergecond !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ DO i = 1, klon !--Calculate the shear value (input for condensation and ice supersat) !--Cell thickness [m] delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * zt(i) * RD IF ( iftop ) THEN ! top shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. & + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. ) ELSEIF ( k .EQ. 1 ) THEN ! surface shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. & + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. ) ELSE ! other layers shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. & - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. & + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. & - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. ) ENDIF ENDDO IF ( ok_ice_supersat ) THEN !--Initialisation IF ( ok_plane_contrail ) THEN IF ( iftop ) THEN dzsed_cont_abv(:) = 0. flsed_cont_abv(:) = 0. Nflsed_cont_abv(:)= 0. cfsed_cont_abv(:) = 0. ELSE dzsed_cont_abv(:) = dzsed_cont(:) flsed_cont_abv(:) = flsed_cont(:) Nflsed_cont_abv(:)= Nflsed_cont(:) cfsed_cont_abv(:) = cfsed_cont(:) ENDIF contfra(:) = 0. qcont(:) = 0. qvcont(:) = 0. Ncont(:) = 0. dzsed_cont(:) = 0. flsed_cont(:) = 0. Nflsed_cont(:)= 0. cfsed_cont(:) = 0. ENDIF IF ( iftop ) THEN dzsed_abv(:) = 0. flsed_abv(:) = 0. cfsed_abv(:) = 0. ELSE dzsed_abv(:) = dzsed(:) flsed_abv(:) = flsed(:) cfsed_abv(:) = cfsed(:) ENDIF dzsed(:) = 0. flsed(:) = 0. cfsed(:) = 0. flauto(:) = 0. DO i = 1, klon pt_pron_clds(i) = .TRUE. ENDDO IF ( .NOT. ok_weibull_warm_clouds ) THEN DO i = 1, klon pt_pron_clds(i) = pt_pron_clds(i) .AND. ( zt(i) .LE. temp_nowater ) ENDDO ENDIF IF ( ok_no_issr_strato ) THEN DO i = 1, klon pt_pron_clds(i) = pt_pron_clds(i) .AND. ( stratomask(i,k) .EQ. 0. ) ENDDO ENDIF totfra_in(:) = 1. qtot_in(:) = zq(:) IF ( ok_nodeep_lscp ) THEN DO i = 1, klon !--If deep convection is activated, the condensation scheme activates !--only in the environment. NB. the clear sky fraction will the be !--maximised by 1. - cfcon(i,k) IF ( pt_pron_clds(i) .AND. ptconv(i,k) ) THEN totfra_in(i) = 1. - cfcon(i,k) qtot_in(i) = zq(i) - ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k) ENDIF ENDDO ENDIF DO i = 1, klon IF ( pt_pron_clds(i) ) THEN IF ( cfcon(i,k) .LT. cfcon_old(i,k) ) THEN !--If deep convection is weakening, we add the clouds that are not anymore !--'in' deep convection to the advected clouds cldfra_in(i) = cldfra_in(i) + ( cfcon_old(i,k) - cfcon(i,k) ) qvc_in(i) = qvc_in(i) + qvcon_old(i,k) * ( cfcon_old(i,k) - cfcon(i,k) ) qice_in(i) = qice_in(i) + qccon_old(i,k) * ( cfcon_old(i,k) - cfcon(i,k) ) ELSEIF ( cldfra_in(i) .GT. eps ) THEN !--Else if deep convection is strengthening, it consumes the existing cloud !--fraction (which does not at this moment represent deep convection) !deepconv_coef = 1. - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) deepconv_coef = MAX(0., 1. - ( cfcon(i,k) - cfcon_old(i,k) ) / cldfra_in(i) ) cldfra_in(i) = cldfra_in(i) * deepconv_coef qvc_in(i) = qvc_in(i) * deepconv_coef qice_in(i) = qice_in(i) * deepconv_coef IF ( ok_plane_contrail ) THEN !--If contrails are activated, their fraction is also reduced when deep !--convection is active cfc_seri(i,k) = cfc_seri(i,k) * deepconv_coef qtc_seri(i,k) = qtc_seri(i,k) * deepconv_coef qic_seri(i,k) = qic_seri(i,k) * deepconv_coef nic_seri(i,k) = nic_seri(i,k) * deepconv_coef ENDIF ENDIF ENDIF ENDDO ENDIF DT(:) = 0. n_i(:)=0 Tbef(:)=zt(:) qlibef(:)=0. Tbefth(:)=tla(:,k)*pspsk(:,k) IF (k .GT. 1) THEN Tbefthm1(:)=tla(:,k-1)*pspsk(:,k-1) ELSE Tbefthm1(:)=Tbefth(:) ENDIF zqth(:)=qta(:,k) zdeltaz(:)=(paprs(:,k)-paprs(:,k+1))/RG/zp(:)*RD*zt(:) ! Treatment of stratiform clouds (lognormale or ice-sursat) or all clouds (including cloudth ! in case of ok_lscp_mergecond) ! We iterate here to take into account the change in qsat(T) and ice fraction ! during the condensation process ! the increment in temperature is calculated using the first principle of ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction ! and a clear sky fraction) ! note that we assume that the vapor in the cloud is at saturation for this calculation DO iter=1,iflag_fisrtilp_qsat+1 keepgoing(:) = .FALSE. DO i=1,klon ! keepgoing = .true. while convergence is not satisfied IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. stratiform_or_all_distrib(i)) THEN ! if not convergence: ! we calculate a new iteration keepgoing(i) = .TRUE. ! P2.2.1> cloud fraction and condensed water mass calculation ! Calculated variables: ! rneb : cloud fraction ! zqn : total water within the cloud ! zcond: mean condensed water within the mesh ! rhcl: clear-sky relative humidity !--------------------------------------------------------------- ! new temperature that only serves in the iteration process: Tbef(i)=Tbef(i)+DT(i) ! total water including mass of precip qtot(i)=zq(i)+zmqc(i) ENDIF ENDDO ! Calculation of saturation specific humidity CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,0,.false.,zqs,zdqs) ! also in thermals CALL calc_qsat_ecmwf(klon,Tbefth,zqth,zp,RTT,0,.false.,zqsth,zdqsth) IF (iflag_icefrac .GE. 1) THEN ! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction ! and cloud condensed water content. Principle from Dietlitcher et al. 2018, GMD ! we make this option works only for the physically-based and tke-dependent param from Raillard et al. ! (iflag_icefrac>=1) CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,1,.false.,zqsl,zdqsl) CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,2,.false.,zqsi,zdqsi) DO i=1,klon zqs(i)=zfice(i)*zqsi(i)+(1.-zfice(i))*zqsl(i) zdqs(i)=zfice(i)*zdqsi(i)+zqsi(i)*dzfice(i)+(1.-zfice(i))*zdqsl(i)-zqsl(i)*dzfice(i) ENDDO ENDIF IF (iflag_icefrac .GE. 2) THEN ! same adjustment for thermals CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,1,.false.,zqslth,zdqslth) CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,2,.false.,zqsith,zdqsith) DO i=1,klon zqsth(i)=zficeth(i)*zqsith(i)+(1.-zficeth(i))*zqslth(i) zdqsth(i)=zficeth(i)*zdqsith(i)+zqsith(i)*dzficeth(i)+(1.-zficeth(i))*zdqslth(i)-zqslth(i)*dzficeth(i) ENDDO ENDIF CALL calc_gammasat(klon,Tbef,qtot,zp,gammasat,dgammasatdt) ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:) ! Cloud condensation based on subgrid pdf !--AB Activates a condensation scheme that allows for !--ice supersaturation and contrails evolution from aviation IF (ok_ice_supersat) THEN !--------------------------------------------- !-- CONDENSATION AND ICE SUPERSATURATION -- !--------------------------------------------- CALL condensation_ice_supersat( & klon, dtime, pplay(:,k), paprs(:,k), paprs(:,k+1), & totfra_in, cldfra_in, qvc_in, qliq_in, qice_in, & shear, tke_dissip(:,k), cell_area, Tbef, qtot_in, zqs, & gammasat, ratqs(:,k), keepgoing, pt_pron_clds, & dzsed_abv, flsed_abv, cfsed_abv, dzsed, flsed, cfsed, flauto, & rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), & dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), & dcf_sed(:,k), dcf_auto(:,k), & dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), & dqi_sed(:,k), dqi_auto(:,k), & dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), & dqvc_sed(:,k), dqvc_auto(:,k), & cfc_seri(:,k), qtc_seri(:,k), qic_seri(:,k), nic_seri(:,k), & flight_dist(:,k), flight_fuel(:,k), & contfra, qcont, qvcont, Ncont, & Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), & AEI_contrails(:,k), AEI_surv_contrails(:,k), & fsurv_contrails(:,k), section_contrails(:,k), & dzsed_cont_abv, flsed_cont_abv, Nflsed_cont_abv, cfsed_cont_abv, & dzsed_cont, flsed_cont, Nflsed_cont, cfsed_cont, & dcfc_ini(:,k), dqic_ini(:,k), dqtc_ini(:,k), dnic_ini(:,k), & dcfc_sub(:,k), dqic_sub(:,k), dqtc_sub(:,k), dnic_sub(:,k), & dcfc_mix(:,k), dqic_mix(:,k), dqtc_mix(:,k), dnic_mix(:,k), & dnic_agg(:,k), dqic_adj(:,k), dqtc_adj(:,k), & dcfc_sed(:,k), dqic_sed(:,k), dqtc_sed(:,k), dnic_sed(:,k), & dcfc_auto(:,k), dqic_auto(:,k), dqtc_auto(:,k), dnic_auto(:,k)) ELSE !--generalised lognormal condensation scheme (Bony and Emanuel 2001) CALL condensation_lognormal( & klon, Tbef, zq, zqs, gammasat, ratqs(:,k), & keepgoing, rneb(:,k), zqn, qvc) ENDIF ! .NOT. ok_ice_supersat ! Volume cloud fraction rneblsvol(:,k)=rneb(:,k) IF (ok_lscp_mergecond) THEN ! in that case we overwrite rneb, rneblsvol and zqn for shallow convective clouds only CALL condensation_cloudth(klon,Tbef,zq,zqth,fraca(:,k), & pspsk(:,k),zp,tla(:,k), & ratqs(:,k),sigma_qtherm(:,k),zqsth,zqs,zqn,rneb(:,k),rnebth(:,k),rneblsvol(:,k), & cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k)) ENDIF ! Cloud phase determination IF (iflag_icefrac .LE. 1) THEN ! "old" phase partitioning depending on temperature and eventually distance to cloud top everywhere IF (iflag_t_glace .GE. 4) THEN CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop) ENDIF CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice) ENDIF IF (iflag_icefrac .EQ. 1) THEN ! phase partitioning depending on turbulence, vertical velocity and ice crystal microphysics ! only in stratiform clouds in the mixed phase regime (Raillard et al. 2025) ! it overwrites temperature-dependent phase partitioning except for convective boundary layer clouds pticefracturb(:) = (fraca(:,k) < min_frac_th_cld) .AND. (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT) DO i=1,klon qincloud(i)=zqn(i) zcf(i)=MIN(MAX(rneb(i,k), 0.),1.) sursat_e(i) = 0. invtau_e(i) = 0. ENDDO CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), qice_in, & ziflcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, & zficeenv, dzficeenv, cldfraliq(:,k),sigma2_icefracturb(:,k),mean_icefracturb(:,k)) DO i=1,klon IF (pticefracturb(i)) THEN zfice(i)=zficeenv(i) dzfice(i)=dzficeenv(i) ENDIF ENDDO ELSEIF (iflag_icefrac .GE. 2) THEN ! compute also phase partitioning also in thermal clouds (neglecting tke in thermals as first assumption) ! moreover, given the upward velocity of thermals, we assume that precipitation falls in the environment ! hence we assume that no snow falls in thermals. ! we activate only in the mixed phase regime not to distrub the cirrus param at very cold T pticefracturb(:) = (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT) !Thermals DO i=1,klon IF (fraca(i,k) .GT. min_frac_th_cld) THEN zcf(i)=MIN(MAX(rnebth(i,k),0.), 1.)/fraca(i,k) qincloud(i)=zqn(i)*fraca(i,k) ELSE zcf(i) = 0. qincloud(i) = 0. ENDIF sursat_e(i)=cloudth_senv(i,k)/zqsi(i) invtau_e(i)=gamma_mixth*MAX(entr_therm(i,k)-detr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i) ENDDO CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbefth, zp, paprs(:,k), paprs(:,k+1), wth(:,k), qice_in, & zeroklon, qincloud, zcf, zeroklon, zeroklon, sursat_e, invtau_e, zqliqth, zqvapclth, zqiceth, & zficeth, dzficeth,cldfraliqth(:,k), sigma2_icefracturbth(:,k), mean_icefracturbth(:,k)) !Environment DO i=1,klon qincloud(i)=zqn(i)*(1.-fraca(i,k)) zcf(i)=MIN(MAX(rneb(i,k)-rnebth(i,k), 0.),1.)/(1.-fraca(i,k)) IF (k .GT. 1) THEN ! evaluate the mixing sursaturation using saturation deficit at level below ! as air pacels detraining into clouds have not (less) seen yet entrainement from above sursat_e(i)=cloudth_sth(i,k-1)/(zqsith(i)+zdqsith(i)*RCPD/RLSTT*(Tbefthm1(i)-Tbefth(i))) ! mixing is assumed to scales with intensity of net detrainment/entrainment rate (D/dz-E/dz) / rho invtau_e(i)=gamma_mixth*MAX(detr_therm(i,k)-entr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i) ELSE sursat_e(i)=0. invtau_e(i)=0. ENDIF ENDDO CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), qice_in, & ziflcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, & zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k)) ! adjust zfice to account for condensates in thermals'fraction DO i=1,klon IF ( zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i) .GT. 0.) THEN zfice(i)=MIN(1., (zqiceth(i)+zqice(i))/(zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i))) dzfice(i)=0. ! dxice/dT=0. always when using icefrac_lscp_turb ENDIF ENDDO ENDIF ! iflag_icefrac DO i=1,klon IF (keepgoing(i)) THEN ! P2.2.2> Approximative calculation of temperature variation DT ! due to condensation. ! Calculated variables: ! dT : temperature change due to condensation !--------------------------------------------------------------- IF (zfice(i).LT.1) THEN cste=RLVTT ELSE cste=RLSTT ENDIF IF ( ok_unadjusted_clouds .OR. ok_unadjusted_contrails ) THEN !--AB We relax the saturation adjustment assumption !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme IF ( rneb(i,k) .GT. eps ) THEN qlibef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k)) ELSE qlibef(i) = 0. ENDIF ELSE qlibef(i)=max(0.,zqn(i)-zqs(i)) ENDIF IF ( ok_ice_sedim ) THEN qice_sedim = (flauto(i) + flsed(i) + flsed_cont(i)) & / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime ! Add the ice that was sedimented, as it is not included in zqn qlibef(i) = qlibef(i) + qice_sedim ENDIF num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT & +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlibef(i) denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) & -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k) & *qlibef(i)*dzfice(i) ! here we update a provisory temperature variable that only serves in the iteration ! process DT(i)=num/denom n_i(i)=n_i(i)+1 ENDIF ! end keepgoing ENDDO ! end loop on i ENDDO ! iter=1,iflag_fisrtilp_qsat+1 ! P2.2> Final quantities calculation ! Calculated variables: ! rneb : cloud fraction ! zcond: mean condensed water in the mesh ! zqn : mean water vapor in the mesh ! zfice: ice fraction in clouds ! zt : temperature ! rhcl : clear-sky relative humidity ! ---------------------------------------------------------------- ! Water vapor and condensed water update, subsequent latent heat exchange for each cloud type !--------------------------------------------------------------------------------------------- DO i=1, klon ! Checks on rneb, rhcl and zqn IF (rneb(i,k) .LE. 0.0) THEN zqn(i) = 0.0 rneb(i,k) = 0.0 zcond(i) = 0.0 rhcl(i,k)=zq(i)/zqs(i) ELSE IF (rneb(i,k) .GE. 1.0) THEN zqn(i) = zq(i) rneb(i,k) = 1.0 IF ( ok_unadjusted_clouds .OR. ok_unadjusted_contrails ) THEN !--AB We relax the saturation adjustment assumption !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme zcond(i) = MAX(0., zqn(i) - qvc(i)) ELSE zcond(i) = MAX(0.0,zqn(i)-zqs(i)) ENDIF rhcl(i,k)=1.0 ELSE IF ( ok_unadjusted_clouds .OR. ok_unadjusted_contrails ) THEN !--AB We relax the saturation adjustment assumption !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i)) ELSE zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k) ENDIF ! following line is very strange and probably wrong rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i) ! Correct calculation of clear-sky relative humidity. To activate once ! issues related to zqn>zq in bi-gaussian clouds are corrected !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat !--This is slighly approximated, the actual formula is !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q) !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content ! rhcl(i,k) = ( zq(i) - zqn(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i) ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param) ENDIF ! water vapor update zq(i) = zq(i) - zcond(i) IF ( ok_ice_sedim ) THEN qice_sedim = (flauto(i) + flsed(i) + flsed_cont(i)) & / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime ! Remove the ice that was sedimented. As it is not included in zqn, ! we only remove it from the total water zq(i) = zq(i) - qice_sedim ! Temperature update due to phase change (sedimented ice was condensed) zt(i) = zt(i) + qice_sedim & * RLSTT / RCPD / ( 1. + RVTMP2 * ( zq(i) + zmqc(i) + zcond(i) ) ) ENDIF ! temperature update due to phase change zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & & * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) & +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) ENDDO ! ---------------------------------------------------------------- ! P3> Precipitation processes, after cloud formation ! - precipitation formation, melting/freezing ! ---------------------------------------------------------------- ! Initiate the quantity of liquid and solid condensates ! Note that in the following, zcond is the total amount of condensates ! including newly formed precipitations (i.e., condensates formed by the ! condensation process), while zoliq is the total amount of condensates in ! the cloud (i.e., on which precipitation processes have applied) DO i = 1, klon zoliq(i) = zcond(i) zoliql(i) = zcond(i) * ( 1. - zfice(i) ) zoliqi(i) = zcond(i) * zfice(i) ENDDO !================================================================ ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R) IF (ok_poprecip) THEN CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), zp, & ctot_vol(:,k), ptconv(:,k), & zt, zq, zoliql, zoliqi, zfice, & rneb(:,k), znebprecipclr, znebprecipcld, & zrfl, zrflclr, zrflcld, & zifl, ziflclr, ziflcld, & qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), & dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), & dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), & dqsmelt(:,k), dqsfreez(:,k) & ) DO i = 1, klon zoliq(i) = zoliql(i) + zoliqi(i) ENDDO !================================================================ ELSE CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, & ctot_vol(:,k), ptconv(:,k), pt_pron_clds, zdqsdT_raw, & zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, & rneb(:,k), znebprecipclr, znebprecipcld, & zneb, tot_zneb, zrho_up, zvelo_up, & zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, & zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) & ) ENDIF ! ok_poprecip IF ( ok_ice_sedim ) THEN zifl(:) = zifl(:) + flauto(:) ziflcld(:) = ziflcld(:) + flauto(:) ENDIF ! End of precipitation processes after cloud formation ! ---------------------------------------------------- !---------------------------------------------------------------------- ! P4> Calculation of cloud condensates amount seen by radiative scheme !---------------------------------------------------------------------- DO i=1,klon IF (ok_poprecip) THEN IF (ok_radocond_snow) THEN radocond(i,k) = zoliq(i) zradoice(i) = zoliqi(i) + qsnowdiag(i,k) ELSE radocond(i,k) = zoliq(i) zradoice(i) = zoliqi(i) ENDIF ELSE radocond(i,k) = zradocond(i) ENDIF ! calculate the percentage of ice in "radocond" so cloud+precip seen ! by the radiation scheme IF (radocond(i,k) .GT. 0.) THEN radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.) ENDIF ENDDO ! ---------------------------------------------------------------- ! P5> Wet scavenging ! ---------------------------------------------------------------- !Scavenging through nucleation in the layer DO i = 1,klon IF(zcond(i).GT.zoliq(i)+1.e-10) THEN beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime ELSE beta(i,k) = 0. ENDIF zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN IF (temp(i,k) .GE. temp_nowater) THEN zalpha_tr = a_tr_sca(3) ELSE zalpha_tr = a_tr_sca(4) ENDIF zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb)) frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi ! Nucleation with a factor of -1 instead of -0.5 zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb)) ENDIF ENDDO ! Scavenging through impaction in the underlying layer DO kk = k-1, 1, -1 DO i = 1, klon IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN IF (temp(i,kk) .GE. temp_nowater) THEN zalpha_tr = a_tr_sca(1) ELSE zalpha_tr = a_tr_sca(2) ENDIF zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb)) frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi ENDIF ENDDO ENDDO !------------------------------------------------------------ ! P6 > write diagnostics and outputs !------------------------------------------------------------ CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,1,.false.,qsatl(:,k),zdqs) CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,2,.false.,qsati(:,k),zdqs) !--AB Write diagnostics and tracers for ice supersaturation IF ( ok_plane_contrail ) THEN DO i = 1, klon IF ( zoliq(i) .LE. 0. ) THEN contfra(i) = 0. qcont(i) = 0. qvcont(i) = 0. Ncont(i) = 0. ENDIF ENDDO cfc_seri(:,k) = contfra(:) nic_seri(:,k) = Ncont(:) !--Ice water content of contrails qice_cont(:,k) = qcont(:) - qvcont(:) IF ( ok_unadjusted_contrails ) THEN qtc_seri(:,k) = qvcont(:) qic_seri(:,k) = qice_cont(:,k) ELSE qtc_seri(:,k) = qcont(:) ENDIF !--Radiative properties contfrarad(:,k) = contfra(:) qradice_cont(:,k) = qice_cont(:,k) ENDIF IF ( ok_ice_supersat ) THEN DO i = 1, klon !--We save the cloud properties that will be advected cf_seri(i,k) = rneb(i,k) qvc_seri(i,k) = qvc(i) !--We keep convective clouds properties in memory, and account for !--the sink of condensed water from precipitation IF ( ptconv(i,k) ) THEN qvcon_old(i,k) = qvcon(i,k) qccon_old(i,k) = qccon(i,k) ELSE qvcon_old(i,k) = 0. qccon_old(i,k) = 0. ENDIF !--If everything was precipitated, the remaining empty cloud is dissipated !--and everything is transfered to the subsaturated clear sky region !--NB. we do not change rneb, as it is a diagnostic only IF ( zoliq(i) .LE. 0. ) THEN cf_seri(i,k) = 0. qvc_seri(i,k) = 0. ENDIF !--Diagnostics gamma_cond(i,k) = gammasat(i) subfra(i,k) = totfra_in(i) - cf_seri(i,k) - issrfra(i,k) qsub(i,k) = qtot_in(i) - qvc_seri(i,k) - qissr(i,k) qcld(i,k) = qvc_seri(i,k) + zoliq(i) IF ( ok_higher_cirrus_cover .AND. pt_pron_clds(i) .AND. .NOT. ptconv(i,k) ) THEN !--Following Brooks et al. (2005) !--This is only valid for cirrus clouds !--We do not apply it do convective clouds IF ( ( rneb(i,k) .GT. eps ) .AND. ( rneb(i,k) .LT. (1. - eps) ) ) THEN dz = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * zt(i) * RD coef_cover = (0.0706 + 0.1274 * shear(i)**0.3015) & * dz**0.7679 / SQRT(cell_area(i))**0.2254 rneb(i,k) = 1. / (1. + EXP(-coef_cover) * (1. / rneb(i,k) - 1.)) IF ( contfrarad(i,k) .GT. eps ) THEN contfrarad(i,k) = cfc_seri(i,k) / cf_seri(i,k) * rneb(i,k) ENDIF ENDIF ENDIF !--Calculation of the ice supersaturated fraction following Lamquin et al (2012) !--methodology: in each layer, we make a maximum random overlap assumption for !--ice supersaturation IF ( ( paprs(i,k) .GT. 10000. ) .AND. ( paprs(i,k) .LE. 15000. ) ) THEN IF ( issrfra100to150UP(i) .GT. ( 1. - eps ) ) THEN issrfra100to150(i) = 1. ELSE issrfra100to150(i) = 1. - ( 1. - issrfra100to150(i) ) * & ( 1. - MAX( issrfra(i,k), issrfra100to150UP(i) ) ) & / ( 1. - issrfra100to150UP(i) ) issrfra100to150UP(i) = issrfra(i,k) ENDIF ELSEIF ( ( paprs(i,k) .GT. 15000. ) .AND. ( paprs(i,k) .LE. 20000. ) ) THEN IF ( issrfra150to200UP(i) .GT. ( 1. - eps ) ) THEN issrfra150to200(i) = 1. ELSE issrfra150to200(i) = 1. - ( 1. - issrfra150to200(i) ) * & ( 1. - MAX( issrfra(i,k), issrfra150to200UP(i) ) ) & / ( 1. - issrfra150to200UP(i) ) issrfra150to200UP(i) = issrfra(i,k) ENDIF ELSEIF ( ( paprs(i,k) .GT. 20000. ) .AND. ( paprs(i,k) .LE. 25000. ) ) THEN IF ( issrfra200to250UP(i) .GT. ( 1. - eps ) ) THEN issrfra200to250(i) = 1. ELSE issrfra200to250(i) = 1. - ( 1. - issrfra200to250(i) ) * & ( 1. - MAX( issrfra(i,k), issrfra200to250UP(i) ) ) & / ( 1. - issrfra200to250UP(i) ) issrfra200to250UP(i) = issrfra(i,k) ENDIF ELSEIF ( ( paprs(i,k) .GT. 25000. ) .AND. ( paprs(i,k) .LE. 30000. ) ) THEN IF ( issrfra250to300UP(i) .GT. ( 1. - eps ) ) THEN issrfra250to300(i) = 1. ELSE issrfra250to300(i) = 1. - ( 1. - issrfra250to300(i) ) * & ( 1. - MAX( issrfra(i,k), issrfra250to300UP(i) ) ) & / ( 1. - issrfra250to300UP(i) ) issrfra250to300UP(i) = issrfra(i,k) ENDIF ELSEIF ( ( paprs(i,k) .GT. 30000. ) .AND. ( paprs(i,k) .LE. 40000. ) ) THEN IF ( issrfra300to400UP(i) .GT. ( 1. - eps ) ) THEN issrfra300to400(i) = 1. ELSE issrfra300to400(i) = 1. - ( 1. - issrfra300to400(i) ) * & ( 1. - MAX( issrfra(i,k), issrfra300to400UP(i) ) ) & / ( 1. - issrfra300to400UP(i) ) issrfra300to400UP(i) = issrfra(i,k) ENDIF ELSEIF ( ( paprs(i,k) .GT. 40000. ) .AND. ( paprs(i,k) .LE. 50000. ) ) THEN IF ( issrfra400to500UP(i) .GT. ( 1. - eps ) ) THEN issrfra400to500(i) = 1. ELSE issrfra400to500(i) = 1. - ( 1. - issrfra400to500(i) ) * & ( 1. - MAX( issrfra(i,k), issrfra400to500UP(i) ) ) & / ( 1. - issrfra400to500UP(i) ) issrfra400to500UP(i) = issrfra(i,k) ENDIF ENDIF ENDDO ENDIF IF ( ok_plane_contrail ) THEN !--Other useful diagnostics DO i = 1, klon !--Very young countrails IF ( dcfc_ini(i,k) .GT. eps ) THEN rho = pplay(i,k) / zt(i) / RD nice_ygcont(i,k) = dnic_ini(i,k) / dcfc_ini(i,k) / 1e6 * rho iwc_ygcont(i,k) = dqic_ini(i,k) / dcfc_ini(i,k) * 1e3 * rho rvol_ygcont(i,k) = (dqic_ini(i,k) / MAX(eps, dnic_ini(i,k)) / rho_ice * 3. / 4. / RPI)**(1./3.) * 1e6 rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG iwp_cont = 1e3 * dqic_ini(i,k) / dcfc_ini(i,k) * rhodz rei_cont = MIN(200., MAX(10., rvol_ygcont(i,k) / eff2vol_radius_contrails)) tau_ygcont(i,k) = iwp_cont*(3.448e-3+2.431/rei_cont) ELSE nice_ygcont(i,k) = missing_val iwc_ygcont(i,k) = missing_val rvol_ygcont(i,k) = missing_val tau_ygcont(i,k) = missing_val ENDIF !--All contrails IF ( cfc_seri(i,k) .GT. 1e-8 ) THEN rho = pplay(i,k) / zt(i) / RD nice_cont(i,k) = nic_seri(i,k) / cfc_seri(i,k) / 1e6 * rho iwc_cont(i,k) = qice_cont(i,k) / cfc_seri(i,k) * 1e3 * rho rvol_cont(i,k) = (qice_cont(i,k) / MAX(eps, nic_seri(i,k)) / rho_ice * 3. / 4. / RPI)**(1./3.) * 1e6 rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG iwp_cont = 1e3 * qice_cont(i,k) / contfrarad(i,k) * rhodz rei_cont = MIN(200., MAX(10., rvol_cont(i,k) / eff2vol_radius_contrails)) tau_cont(i,k) = iwp_cont*(3.448e-3+2.431/rei_cont) IF ( tau_cont(i,k) .GT. 0.05 ) THEN nice_vscont(i,k) = nice_cont(i,k) iwc_vscont(i,k) = iwc_cont(i,k) rvol_vscont(i,k) = rvol_cont(i,k) tau_vscont(i,k) = tau_cont(i,k) ELSE nice_vscont(i,k) = missing_val iwc_vscont(i,k) = missing_val rvol_vscont(i,k) = missing_val tau_vscont(i,k) = missing_val ENDIF ELSE nice_cont(i,k) = missing_val iwc_cont(i,k) = missing_val rvol_cont(i,k) = missing_val tau_cont(i,k) = missing_val nice_vscont(i,k) = missing_val iwc_vscont(i,k) = missing_val rvol_vscont(i,k) = missing_val tau_vscont(i,k) = missing_val ENDIF ENDDO ENDIF ! Outputs: !------------------------------- ! Precipitation fluxes at layer interfaces ! + precipitation fractions + ! temperature and water species tendencies DO i = 1, klon psfl(i,k)=zifl(i) prfl(i,k)=zrfl(i) pfraclr(i,k)=znebprecipclr(i) pfracld(i,k)=znebprecipcld(i) distcltop(i,k)=zdistcltop(i) temp_cltop(i,k)=ztemp_cltop(i) d_q(i,k) = zq(i) - qt(i,k) d_t(i,k) = zt(i) - temp(i,k) IF (ok_bug_phase_lscp) THEN d_ql(i,k) = (1-zfice(i))*zoliq(i) d_qi(i,k) = zfice(i)*zoliq(i) ELSE d_ql(i,k) = zoliql(i) d_qi(i,k) = zoliqi(i) ENDIF ENDDO ENDDO ! loop on k from top to bottom ! Rain or snow at the surface (depending on the first layer temperature) DO i = 1, klon snow(i) = zifl(i) rain(i) = zrfl(i) ! c_iso final output ENDDO IF ( ok_ice_sedim ) THEN DO i = 1, klon snow(i) = snow(i) + flsed(i) + flsed_cont(i) ENDDO ENDIF IF (ncoreczq>0) THEN WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.' ENDIF RETURN END SUBROUTINE lscp !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ END MODULE lmdz_lscp_main