MODULE lmdz_lscp IMPLICIT NONE CONTAINS !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE lscp(klon, klev, dtime, missing_val, & paprs, pplay, temp, qt, qice_save, ptconv, ratqs, & d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri, & pfraclr, pfracld, & cldfraliq, sigma2_icefracturb, mean_icefracturb, & radocond, radicefrac, rain, snow, & frac_impa, frac_nucl, beta, & prfl, psfl, rhcl, qta, fraca, & tv, pspsk, tla, thl, iflag_cld_th, & iflag_ice_thermo, ok_ice_sursat, qsatl, qsats, & distcltop, temp_cltop, tke, tke_dissip, & qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & Tcontr, qcontr, qcontr2, fcontrN, fcontrP, & 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_cloudth, ONLY: cloudth, cloudth_v3, cloudth_v6, cloudth_mpc USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat USE lmdz_lscp_tools, ONLY: icefrac_lscp, icefrac_lscp_turb USE lmdz_lscp_tools, ONLY: fallice_velocity, distance_to_cloud_top USE ice_sursat_mod, ONLY: ice_sursat USE lmdz_lscp_poprecip, ONLY: poprecip_precld, poprecip_postcld ! Use du module lmdz_lscp_ini contenant les constantes USE lmdz_lscp_ini, ONLY: prt_level, lunout USE lmdz_lscp_ini, ONLY: seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min USE lmdz_lscp_ini, ONLY: ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc USE lmdz_lscp_ini, ONLY: iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min USE lmdz_lscp_ini, ONLY: coef_eva, coef_sub, cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con USE lmdz_lscp_ini, ONLY: iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo USE lmdz_lscp_ini, ONLY: iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG USE lmdz_lscp_ini, ONLY: ok_poprecip USE lmdz_lscp_ini, ONLY: iflag_icefrac USE lmdz_abort_physic, ONLY: abort_physic IMPLICIT NONE !=============================================================================== ! VARIABLES DECLARATION !=============================================================================== ! INPUT VARIABLES: !----------------- INTEGER, INTENT(IN) :: klon, klev ! number of horizontal grid points and vertical levels 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) :: qt ! total specific humidity (in vapor phase in input) [kg/kg] REAL, DIMENSION(klon, klev), INTENT(IN) :: qice_save ! ice specific from previous time step [kg/kg] INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics ! CR: if iflag_ice_thermo=2, ONLY convection is active LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated LOGICAL, DIMENSION(klon, klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active !Inputs associated with thermal plumes 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 temperature within thermals [K] REAL, DIMENSION(klon, klev), INTENT(IN) :: tke !--turbulent kinetic energy [m2/s2] REAL, DIMENSION(klon, klev), INTENT(IN) :: tke_dissip !--TKE dissipation [m2/s3] ! INPUT/OUTPUT variables !------------------------ REAL, DIMENSION(klon, klev), INTENT(INOUT) :: thl ! liquid potential temperature [K] REAL, DIMENSION(klon, klev), INTENT(INOUT) :: ratqs ! function of pressure that sets the large-scale ! Input sursaturation en glace REAL, DIMENSION(klon, klev), INTENT(INOUT) :: rneb_seri ! fraction nuageuse en memoire ! 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 [-] 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) :: 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), INTENT(OUT) :: qsatl ! saturation specific humidity wrt liquid [kg/kg] REAL, DIMENSION(klon, klev), INTENT(OUT) :: qsats ! saturation specific humidity wrt ice [kg/kg] 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 supersaturation and contrails parameterisation REAL, DIMENSION(klon, klev), INTENT(OUT) :: qclr ! specific total water content in clear sky 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) :: qss ! specific total water content in supersat region [kg/kg] REAL, DIMENSION(klon, klev), INTENT(OUT) :: qvc ! specific vapor content in clouds [kg/kg] REAL, DIMENSION(klon, klev), INTENT(OUT) :: rnebclr ! mesh fraction of clear sky [-] REAL, DIMENSION(klon, klev), INTENT(OUT) :: rnebss ! mesh fraction of ISSR [-] REAL, DIMENSION(klon, klev), INTENT(OUT) :: gamma_ss ! coefficient governing the ice nucleation RHi threshold [-] REAL, DIMENSION(klon, klev), INTENT(OUT) :: Tcontr ! threshold temperature for contrail formation [K] REAL, DIMENSION(klon, klev), INTENT(OUT) :: qcontr ! threshold humidity for contrail formation [kg/kg] REAL, DIMENSION(klon, klev), INTENT(OUT) :: qcontr2 ! // (2nd expression more consistent with LMDZ expression of q) REAL, DIMENSION(klon, klev), INTENT(OUT) :: fcontrN ! fraction of grid favourable to non-persistent contrails REAL, DIMENSION(klon, klev), INTENT(OUT) :: fcontrP ! fraction of grid favourable to persistent contrails 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 ! 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] ! LOCAL VARIABLES: !---------------- REAL, DIMENSION(klon) :: qsl, qsi ! saturation threshold at current vertical level REAL :: zct, zcl, zexpo REAL, DIMENSION(klon, klev) :: ctot REAL, DIMENSION(klon, klev) :: ctot_vol REAL, DIMENSION(klon) :: zqs, zdqs REAL :: zdelta, zcor, zcvm5 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, qlbef, DT ! temperature, humidity and temp. variation during lognormal iteration REAL :: num, denom REAL :: cste REAL, DIMENSION(klon) :: zpdf_sig, zpdf_k, zpdf_delta ! lognormal parameters REAL, DIMENSION(klon) :: Zpdf_a, zpdf_b, zpdf_e1, zpdf_e2 ! lognormal intermediate variables REAL :: erf REAL, DIMENSION(klon) :: zfice_th REAL, DIMENSION(klon) :: qcloud, qincloud_mpc REAL, DIMENSION(klon) :: zrfl, zrfln REAL :: zqev, zqevt REAL, DIMENSION(klon) :: zifl, zifln, ziflprev REAL :: zqev0, zqevi, zqevti REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn REAL, DIMENSION(klon) :: zoliql, zoliqi REAL, DIMENSION(klon) :: zt REAL, DIMENSION(klon, klev) :: zrho REAL, DIMENSION(klon) :: zdz, iwc REAL :: zchau, zfroi REAL, DIMENSION(klon) :: zfice, zneb, znebprecip REAL :: zmelt, zrain, zsnow, zprecip REAL, DIMENSION(klon) :: dzfice REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb REAL :: zsolid REAL, DIMENSION(klon) :: qtot, qzero REAL, DIMENSION(klon) :: dqsl, dqsi REAL :: smallestreal ! Variables for Bergeron process REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl REAL, DIMENSION(klon) :: zqpreci, zqprecl ! Variables precipitation energy conservation REAL, DIMENSION(klon) :: zmqc REAL :: zalpha_tr REAL :: zfrac_lessi REAL, DIMENSION(klon) :: zprec_cond REAL :: zmair REAL :: zcpair, zcpeau REAL, DIMENSION(klon) :: zlh_solid REAL, DIMENSION(klon) :: ztupnew REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl REAL :: zm_solid ! for liquid -> solid conversion REAL, DIMENSION(klon) :: zrflclr, zrflcld REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr REAL, DIMENSION(klon) :: ziflclr, ziflcld REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr REAL, DIMENSION(klon, klev) :: velo REAL :: vr, ffallv REAL :: qlmpc, qimpc, rnebmpc REAL, DIMENSION(klon, klev) :: radocondi, radocondl REAL :: effective_zneb REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl ! for icefrac_lscp_turb INTEGER i, k, n, kk, iter INTEGER, DIMENSION(klon) :: n_i INTEGER ncoreczq INTEGER, DIMENSION(klon, klev) :: mpc_bl_points LOGICAL iftop LOGICAL, DIMENSION(klon) :: lognormale LOGICAL, DIMENSION(klon) :: keepgoing CHARACTER (len = 20) :: modname = 'lscp' CHARACTER (len = 80) :: abort_message !=============================================================================== ! INITIALISATION !=============================================================================== ! Few initial checks IF (iflag_fisrtilp_qsat < 0) THEN abort_message = 'lscp cannot be used with iflag_fisrtilp<0' CALL abort_physic(modname, abort_message, 1) ENDIF ! Few initialisations znebprecip(:) = 0.0 ctot_vol(1:klon, 1:klev) = 0.0 rneblsvol(1:klon, 1:klev) = 0.0 smallestreal = 1.e-9 znebprecipclr(:) = 0.0 znebprecipcld(:) = 0.0 mpc_bl_points(:, :) = 0 IF (prt_level>9) WRITE(lunout, *) 'NUAGES4 A. JAM' ! 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 pfraclr(:, :) = 0.0 pfracld(:, :) = 0.0 cldfraliq(:, :) = 0. sigma2_icefracturb(:, :) = 0. mean_icefracturb(:, :) = 0. radocond(:, :) = 0.0 radicefrac(:, :) = 0.0 frac_nucl(:, :) = 1.0 frac_impa(:, :) = 1.0 rain(:) = 0.0 snow(:) = 0.0 zoliq(:) = 0.0 zfice(:) = 0.0 dzfice(:) = 0.0 zfice_turb(:) = 0.0 dzfice_turb(:) = 0.0 zqprecl(:) = 0.0 zqpreci(:) = 0.0 zrfl(:) = 0.0 zifl(:) = 0.0 ziflprev(:) = 0.0 zneb(:) = seuil_neb zrflclr(:) = 0.0 ziflclr(:) = 0.0 zrflcld(:) = 0.0 ziflcld(:) = 0.0 tot_zneb(:) = 0.0 tot_znebn(:) = 0.0 d_tot_zneb(:) = 0.0 qzero(:) = 0.0 zdistcltop(:) = 0.0 ztemp_cltop(:) = 0.0 ztupnew(:) = 0.0 !--ice supersaturation gamma_ss(:, :) = 1.0 qss(:, :) = 0.0 rnebss(:, :) = 0.0 Tcontr(:, :) = missing_val qcontr(:, :) = missing_val qcontr2(:, :) = missing_val fcontrN(:, :) = 0.0 fcontrP(:, :) = 0.0 distcltop(:, :) = 0. temp_cltop(:, :) = 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. !c_iso: variable initialisation for iso !=============================================================================== ! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM !=============================================================================== ncoreczq = 0 DO k = klev, 1, -1 IF (k<=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) 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 !================================================================ ! 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), pplay(:, k), & zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, & zqvapclr, zqupnew, & zrfl, zrflclr, zrflcld, & zifl, ziflclr, ziflcld, & dqreva(:, k), dqssub(:, k) & ) !================================================================ ELSE ! -------------------------------------------------------------------- ! P1> Thermalization of precipitation falling from the overlying layer ! -------------------------------------------------------------------- ! Computes air temperature variation due to enthalpy transported by ! precipitation. Precipitation is then thermalized with the air in the ! layer. ! The precipitation should remain thermalized throughout the different ! thermodynamical transformations. ! The corresponding water mass should ! be added when calculating the layer's enthalpy change with ! temperature ! See lmdzpedia page todoan ! todoan: check consistency with ice phase ! todoan: understand why several steps ! --------------------------------------------------------------------- IF (iftop) THEN DO i = 1, klon zmqc(i) = 0. ENDDO ELSE DO i = 1, klon zmair = (paprs(i, k) - paprs(i, k + 1)) / RG ! no condensed water so cp=cp(vapor+dry air) ! RVTMP2=rcpv/rcpd-1 zcpair = RCPD * (1.0 + RVTMP2 * zq(i)) zcpeau = RCPD * RVTMP2 ! zmqc: precipitation mass that has to be thermalized with ! layer's air so that precipitation at the ground has the ! same temperature as the lowermost layer zmqc(i) = (zrfl(i) + zifl(i)) * dtime / zmair ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer zt(i) = (ztupnew(i) * zmqc(i) * zcpeau + zcpair * zt(i)) & / (zcpair + zmqc(i) * zcpeau) ENDDO ENDIF ! -------------------------------------------------------------------- ! P2> Precipitation evaporation/sublimation/melting ! -------------------------------------------------------------------- ! A part of the precipitation coming from above is evaporated/sublimated/melted. ! -------------------------------------------------------------------- IF (iflag_evap_prec>=1) THEN ! Calculation of saturation specific humidity ! depending on temperature: CALL calc_qsat_ecmwf(klon, zt(:), qzero(:), pplay(:, k), RTT, 0, .FALSE., zqs(:), zdqs(:)) ! wrt liquid water CALL calc_qsat_ecmwf(klon, zt(:), qzero(:), pplay(:, k), RTT, 1, .FALSE., qsl(:), dqsl(:)) ! wrt ice CALL calc_qsat_ecmwf(klon, zt(:), qzero(:), pplay(:, k), RTT, 2, .FALSE., qsi(:), dqsi(:)) DO i = 1, klon ! if precipitation IF (zrfl(i) + zifl(i)>0.) THEN ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4). ! c_iso: likely important to distinguish cs from neb isotope precipitation IF (iflag_evap_prec>=4) THEN zrfl(i) = zrflclr(i) zifl(i) = ziflclr(i) ENDIF IF (iflag_evap_prec==1) THEN znebprecip(i) = zneb(i) ELSE znebprecip(i) = MAX(zneb(i), znebprecip(i)) ENDIF IF (iflag_evap_prec>4) THEN ! Max evaporation not to saturate the clear sky precip fraction ! i.e. the fraction where evaporation occurs zqev0 = MAX(0.0, (zqs(i) - zq(i)) * znebprecipclr(i)) ELSEIF (iflag_evap_prec == 4) THEN ! Max evaporation not to saturate the whole mesh ! Pay attention -> lead to unrealistic and excessive evaporation zqev0 = MAX(0.0, zqs(i) - zq(i)) ELSE ! Max evap not to saturate the fraction below the cloud zqev0 = MAX(0.0, (zqs(i) - zq(i)) * znebprecip(i)) ENDIF ! Evaporation of liquid precipitation coming from above ! dP/dz=beta*(1-q/qsat)*sqrt(P) ! formula from Sundquist 1988, Klemp & Wilhemson 1978 ! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4) IF (iflag_evap_prec==3) THEN zqevt = znebprecip(i) * coef_eva * (1.0 - zq(i) / qsl(i)) & * SQRT(zrfl(i) / max(1.e-4, znebprecip(i))) & * (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG ELSE IF (iflag_evap_prec>=4) THEN zqevt = znebprecipclr(i) * coef_eva * (1.0 - zq(i) / qsl(i)) & * SQRT(zrfl(i) / max(1.e-8, znebprecipclr(i))) & * (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG ELSE zqevt = 1. * coef_eva * (1.0 - zq(i) / qsl(i)) * SQRT(zrfl(i)) & * (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG ENDIF zqevt = MAX(0.0, MIN(zqevt, zrfl(i))) & * RG * dtime / (paprs(i, k) - paprs(i, k + 1)) ! sublimation of the solid precipitation coming from above IF (iflag_evap_prec==3) THEN zqevti = znebprecip(i) * coef_sub * (1.0 - zq(i) / qsi(i)) & * SQRT(zifl(i) / max(1.e-4, znebprecip(i))) & * (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG ELSE IF (iflag_evap_prec>=4) THEN zqevti = znebprecipclr(i) * coef_sub * (1.0 - zq(i) / qsi(i)) & * SQRT(zifl(i) / max(1.e-8, znebprecipclr(i))) & * (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG ELSE zqevti = 1. * coef_sub * (1.0 - zq(i) / qsi(i)) * SQRT(zifl(i)) & * (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG ENDIF zqevti = MAX(0.0, MIN(zqevti, zifl(i))) & * RG * dtime / (paprs(i, k) - paprs(i, k + 1)) ! A. JAM ! Evaporation limit: we ensure that the layer's fraction below ! the cloud or the whole mesh (depending on iflag_evap_prec) ! does not reach saturation. In this case, we ! redistribute zqev0 conserving the ratio liquid/ice IF (zqevt + zqevti>zqev0) THEN zqev = zqev0 * zqevt / (zqevt + zqevti) zqevi = zqev0 * zqevti / (zqevt + zqevti) ELSE zqev = zqevt zqevi = zqevti ENDIF ! New solid and liquid precipitation fluxes after evap and sublimation zrfln(i) = Max(0., zrfl(i) - zqev * (paprs(i, k) - paprs(i, k + 1)) & / RG / dtime) zifln(i) = Max(0., zifl(i) - zqevi * (paprs(i, k) - paprs(i, k + 1)) & / RG / dtime) ! vapor, temperature, precip fluxes update ! vapor is updated after evaporation/sublimation (it is increased) zq(i) = zq(i) - (zrfln(i) + zifln(i) - zrfl(i) - zifl(i)) & * (RG / (paprs(i, k) - paprs(i, k + 1))) * dtime ! zmqc is the total condensed water in the precip flux (it is decreased) zmqc(i) = zmqc(i) + (zrfln(i) + zifln(i) - zrfl(i) - zifl(i)) & * (RG / (paprs(i, k) - paprs(i, k + 1))) * dtime ! air and precip temperature (i.e., gridbox temperature) ! is updated due to latent heat cooling zt(i) = zt(i) + (zrfln(i) - zrfl(i)) & * (RG / (paprs(i, k) - paprs(i, k + 1))) * dtime & * RLVTT / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i))) & + (zifln(i) - zifl(i)) & * (RG / (paprs(i, k) - paprs(i, k + 1))) * dtime & * RLSTT / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i))) ! New values of liquid and solid precipitation zrfl(i) = zrfln(i) zifl(i) = zifln(i) ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout) ! due to evap + sublim IF (iflag_evap_prec>=4) THEN zrflclr(i) = zrfl(i) ziflclr(i) = zifl(i) IF(zrflclr(i) + ziflclr(i)<=0) THEN znebprecipclr(i) = 0.0 ENDIF zrfl(i) = zrflclr(i) + zrflcld(i) zifl(i) = ziflclr(i) + ziflcld(i) ENDIF ! c_iso duplicate for isotopes or loop on isotopes ! Melting: zmelt = ((zt(i) - RTT) / (ztfondue - RTT)) ! JYG ! precip fraction that is melted zmelt = MIN(MAX(zmelt, 0.), 1.) ! update of rainfall and snowfall due to melting IF (iflag_evap_prec>=4) THEN zrflclr(i) = zrflclr(i) + zmelt * ziflclr(i) zrflcld(i) = zrflcld(i) + zmelt * ziflcld(i) zrfl(i) = zrflclr(i) + zrflcld(i) ELSE zrfl(i) = zrfl(i) + zmelt * zifl(i) ENDIF ! c_iso: melting of isotopic precipi with zmelt (no fractionation) ! Latent heat of melting because of precipitation melting ! NB: the air + precip temperature is simultaneously updated zt(i) = zt(i) - zifl(i) * zmelt * (RG * dtime) / (paprs(i, k) - paprs(i, k + 1)) & * RLMLT / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i))) IF (iflag_evap_prec>=4) THEN ziflclr(i) = ziflclr(i) * (1. - zmelt) ziflcld(i) = ziflcld(i) * (1. - zmelt) zifl(i) = ziflclr(i) + ziflcld(i) ELSE zifl(i) = zifl(i) * (1. - zmelt) ENDIF ELSE ! if no precip, we reinitialize the cloud fraction used for the precip to 0 znebprecip(i) = 0. ENDIF ! (zrfl(i)+zifl(i).GT.0.) ENDDO ! loop on klon ENDIF ! (iflag_evap_prec>=1) ENDIF ! (ok_poprecip) ! -------------------------------------------------------------------- ! End precip evaporation ! -------------------------------------------------------------------- ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter !------------------------------------------------------- qtot(:) = zq(:) + zmqc(:) CALL calc_qsat_ecmwf(klon, zt(:), qtot(:), pplay(:, k), 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) < 1.e-15) THEN ncoreczq = ncoreczq + 1 zq(i) = 1.e-15 ENDIF ! c_iso: do something similar for isotopes ENDDO ! -------------------------------------------------------------------- ! P2> Cloud formation !--------------------------------------------------------------------- ! Unlike fisrtilp, 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 qincloud_mpc gammasat(:) = 1. qincloud_mpc(:) = 0. IF (iflag_cld_th>=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<=2) THEN ! Old version of Arnaud Jam CALL cloudth(klon, klev, k, tv, & zq, qta, fraca, & qcloud, ctot, pspsk, paprs, pplay, tla, thl, & ratqs, zqs, temp, & cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv) ELSEIF (iflag_cloudth_vert>=3 .AND. iflag_cloudth_vert<=5) THEN ! Default version of Arnaud Jam CALL cloudth_v3(klon, klev, k, tv, & zq, qta, fraca, & qcloud, ctot, ctot_vol, pspsk, paprs, pplay, tla, thl, & ratqs, zqs, temp, & cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv) ELSEIF (iflag_cloudth_vert==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, & qcloud, ctot, ctot_vol, pspsk, paprs, pplay, tla, thl, & ratqs, zqs, temp, & cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv) ELSEIF (iflag_cloudth_vert == 7) THEN ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment ! for boundary-layer mixed phase clouds CALL cloudth_mpc(klon, klev, k, mpc_bl_points, zt, zq, qta(:, k), fraca(:, k), & pspsk(:, k), paprs(:, k + 1), paprs(:, k), pplay(:, k), tla(:, k), & ratqs(:, k), qcloud, qincloud_mpc, zfice_th, ctot(:, k), ctot_vol(:, k), & cloudth_sth(:, k), cloudth_senv(:, k), cloudth_sigmath(:, k), cloudth_sigmaenv(:, k)) ENDIF DO i = 1, klon rneb(i, k) = ctot(i, k) rneblsvol(i, k) = ctot_vol(i, k) zqn(i) = qcloud(i) ENDDO ENDIF IF (iflag_cld_th <= 4) THEN ! lognormal lognormale(:) = .TRUE. ELSEIF (iflag_cld_th >= 6) THEN ! lognormal distribution when no thermals lognormale(:) = fraca(:, k) < min_frac_th_cld ELSE ! When iflag_cld_th=5, we always assume ! bi-gaussian distribution lognormale(:) = .FALSE. ENDIF DT(:) = 0. n_i(:) = 0 Tbef(:) = zt(:) qlbef(:) = 0. ! Treatment of non-boundary layer clouds (lognormale) ! condensation with qsat(T) variation (adaptation) ! Iterative resolution to converge towards qsat ! with update of temperature, ice fraction and qsat at ! each iteration ! todoan -> sensitivity to iflag_fisrtilp_qsat DO iter = 1, iflag_fisrtilp_qsat + 1 DO i = 1, klon ! keepgoing = .TRUE. while convergence is not satisfied keepgoing(i) = ABS(DT(i))>DDT0 IF ((keepgoing(i) .OR. (n_i(i) == 0)) .AND. lognormale(i)) THEN ! if not convergence: ! 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) ! Rneb, qzn and zcond for lognormal PDFs qtot(i) = zq(i) + zmqc(i) ENDIF ENDDO ! Calculation of saturation specific humidity and ice fraction CALL calc_qsat_ecmwf(klon, Tbef(:), qtot(:), pplay(:, k), RTT, 0, .FALSE., zqs(:), zdqs(:)) CALL calc_gammasat(klon, Tbef(:), qtot(:), pplay(:, k), gammasat(:), dgammasatdt(:)) ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs zdqs(:) = gammasat(:) * zdqs(:) + zqs(:) * dgammasatdt(:) ! cloud phase determination IF (iflag_t_glace>=4) THEN ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top 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(:)) DO i = 1, klon !todoan : check if loop in i is needed IF ((keepgoing(i) .OR. (n_i(i) == 0)) .AND. lognormale(i)) THEN zpdf_sig(i) = ratqs(i, k) * zq(i) zpdf_k(i) = -sqrt(log(1. + (zpdf_sig(i) / zq(i))**2)) zpdf_delta(i) = log(zq(i) / (gammasat(i) * zqs(i))) zpdf_a(i) = zpdf_delta(i) / (zpdf_k(i) * sqrt(2.)) zpdf_b(i) = zpdf_k(i) / (2. * sqrt(2.)) zpdf_e1(i) = zpdf_a(i) - zpdf_b(i) zpdf_e1(i) = sign(min(ABS(zpdf_e1(i)), 5.), zpdf_e1(i)) zpdf_e1(i) = 1. - erf(zpdf_e1(i)) zpdf_e2(i) = zpdf_a(i) + zpdf_b(i) zpdf_e2(i) = sign(min(ABS(zpdf_e2(i)), 5.), zpdf_e2(i)) zpdf_e2(i) = 1. - erf(zpdf_e2(i)) IF ((.NOT.ok_ice_sursat).OR.(Tbef(i)>t_glace_min)) THEN IF (zpdf_e1(i)<1.e-10) THEN rneb(i, k) = 0. zqn(i) = gammasat(i) * zqs(i) ELSE rneb(i, k) = 0.5 * zpdf_e1(i) zqn(i) = zq(i) * zpdf_e2(i) / zpdf_e1(i) ENDIF rnebss(i, k) = 0.0 !--added by OB (needed because of the convergence loop on time) fcontrN(i, k) = 0.0 !--idem fcontrP(i, k) = 0.0 !--idem qss(i, k) = 0.0 !--idem ELSE ! in case of ice supersaturation by Audran !------------------------------------ ! ICE SUPERSATURATION !------------------------------------ CALL ice_sursat(pplay(i, k), paprs(i, k) - paprs(i, k + 1), dtime, i, k, temp(i, k), zq(i), & gamma_ss(i, k), zqs(i), Tbef(i), rneb_seri(i, k), ratqs(i, k), & rneb(i, k), zqn(i), rnebss(i, k), qss(i, k), & Tcontr(i, k), qcontr(i, k), qcontr2(i, k), fcontrN(i, k), fcontrP(i, k)) ENDIF ! ((flag_ice_sursat.EQ.0).OR.(Tbef(i).gt.t_glace_min)) ! If vertical heterogeneity, change fraction by volume as well IF (iflag_cloudth_vert>=3) THEN ctot_vol(i, k) = rneb(i, k) rneblsvol(i, k) = ctot_vol(i, k) ENDIF ! P2.2.2> Approximative calculation of temperature variation DT ! due to condensation. ! Calculated variables: ! dT : temperature change due to condensation !--------------------------------------------------------------- IF (zfice(i)<1) THEN cste = RLVTT ELSE cste = RLSTT ENDIF ! LEA_R : check formule qlbef(i) = max(0., zqn(i) - zqs(i)) num = -Tbef(i) + zt(i) + rneb(i, k) * ((1 - zfice(i)) * RLVTT & + zfice(i) * RLSTT) / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i))) * qlbef(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) & * qlbef(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.3> 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 ! ---------------------------------------------------------------- ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top IF (iflag_t_glace>=4) THEN CALL distance_to_cloud_top(klon, klev, k, temp, pplay, paprs, rneb, zdistcltop, ztemp_cltop) distcltop(:, k) = zdistcltop(:) temp_cltop(:, k) = ztemp_cltop(:) ENDIF ! Partition function depending on temperature CALL icefrac_lscp(klon, Tbef, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice) ! Partition function depending on tke for non shallow-convective clouds IF (iflag_icefrac >= 1) THEN CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:, k), paprs(:, k), paprs(:, k + 1), qice_save(:, k), ziflcld, zqn, & rneb(:, k), tke(:, k), tke_dissip(:, k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:, k), sigma2_icefracturb(:, k), mean_icefracturb(:, k)) ENDIF ! Water vapor update, Phase determination and subsequent latent heat exchange DO i = 1, klon ! Overwrite phase partitioning in boundary layer mixed phase clouds when the ! iflag_cloudth_vert=7 and specific param is activated IF (mpc_bl_points(i, k) > 0) THEN zcond(i) = MAX(0.0, qincloud_mpc(i)) * rneb(i, k) ! following line is very strange and probably wrong rhcl(i, k) = (zqs(i) + zq(i)) / 2. / zqs(i) ! water vapor update and partition function if thermals zq(i) = zq(i) - zcond(i) zfice(i) = zfice_th(i) ELSE ! Checks on rneb, rhcl and zqn IF (rneb(i, k) <= 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) >= 1.0) THEN zqn(i) = zq(i) rneb(i, k) = 1.0 zcond(i) = MAX(0.0, zqn(i) - zqs(i)) rhcl(i, k) = 1.0 ELSE zcond(i) = MAX(0.0, zqn(i) - zqs(i)) * rneb(i, k) ! following line is very strange and probably wrong: rhcl(i, k) = (zqs(i) + zq(i)) / 2. / zqs(i) ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param) IF (iflag_icefrac >= 1) THEN IF (lognormale(i)) THEN zcond(i) = zqliq(i) + zqice(i) zfice(i) = zfice_turb(i) rhcl(i, k) = zqvapcl(i) * rneb(i, k) + (zq(i) - zqn(i)) * (1. - rneb(i, k)) ENDIF ENDIF ENDIF ! water vapor update zq(i) = zq(i) - zcond(i) ENDIF ! c_iso : routine that computes in-cloud supersaturation ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input) ! 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 ! If vertical heterogeneity, change volume fraction IF (iflag_cloudth_vert >= 3) THEN ctot_vol(1:klon, k) = min(max(ctot_vol(1:klon, k), 0.), 1.) rneblsvol(1:klon, k) = ctot_vol(1:klon, k) ENDIF ! ---------------------------------------------------------------- ! P3> Precipitation formation ! ---------------------------------------------------------------- !================================================================ IF (ok_poprecip) THEN DO i = 1, klon zoliql(i) = zcond(i) * (1. - zfice(i)) zoliqi(i) = zcond(i) * zfice(i) ENDDO CALL poprecip_postcld(klon, dtime, paprs(:, k), paprs(:, k + 1), pplay(:, k), & 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) IF (zoliq(i) > 0.) THEN zfice(i) = zoliqi(i) / zoliq(i) ELSE zfice(i) = 0.0 ENDIF ! calculation of specific content of condensates seen by radiative scheme IF (ok_radocond_snow) THEN radocond(i, k) = zoliq(i) radocondl(i, k) = radocond(i, k) * (1. - zfice(i)) radocondi(i, k) = radocond(i, k) * zfice(i) + qsnowdiag(i, k) ELSE radocond(i, k) = zoliq(i) radocondl(i, k) = radocond(i, k) * (1. - zfice(i)) radocondi(i, k) = radocond(i, k) * zfice(i) ENDIF ENDDO !================================================================ ELSE ! LTP: IF (iflag_evap_prec >= 4) THEN !Partitionning between precipitation coming from clouds and that coming from CS !0) Calculate tot_zneb, total cloud fraction above the cloud !assuming a maximum-random overlap (voir Jakob and Klein, 2000) DO i = 1, klon tot_znebn(i) = 1. - (1. - tot_zneb(i)) * (1 - max(rneb(i, k), zneb(i))) & / (1. - min(zneb(i), 1. - smallestreal)) d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i) tot_zneb(i) = tot_znebn(i) !1) Cloudy to clear air d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i, k), znebprecipcld(i)) IF (znebprecipcld(i) > 0.) THEN d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i) / znebprecipcld(i) * zrflcld(i) d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i) / znebprecipcld(i) * ziflcld(i) ELSE d_zrfl_cld_clr(i) = 0. d_zifl_cld_clr(i) = 0. ENDIF !2) Clear to cloudy air d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i, k) - d_tot_zneb(i) - zneb(i))) IF (znebprecipclr(i) > 0) THEN d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i) / znebprecipclr(i) * zrflclr(i) d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i) / znebprecipclr(i) * ziflclr(i) ELSE d_zrfl_clr_cld(i) = 0. d_zifl_clr_cld(i) = 0. ENDIF !Update variables znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i) znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i) zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i) ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i) zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i) ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i) ! c_iso do the same thing for isotopes precip ENDDO ENDIF ! Autoconversion ! ------------------------------------------------------------------------------- ! Initialisation of zoliq and radocond variables DO i = 1, klon zoliq(i) = zcond(i) zoliqi(i) = zoliq(i) * zfice(i) zoliql(i) = zoliq(i) * (1. - zfice(i)) ! c_iso : initialisation of zoliq* also for isotopes zrho(i, k) = pplay(i, k) / zt(i) / RD zdz(i) = (paprs(i, k) - paprs(i, k + 1)) / (zrho(i, k) * RG) iwc(i) = 0. zneb(i) = MAX(rneb(i, k), seuil_neb) radocond(i, k) = zoliq(i) / REAL(niter_lscp + 1) radocondi(i, k) = zoliq(i) * zfice(i) / REAL(niter_lscp + 1) radocondl(i, k) = zoliq(i) * (1. - zfice(i)) / REAL(niter_lscp + 1) ENDDO DO n = 1, niter_lscp DO i = 1, klon IF (rneb(i, k)>0.0) THEN iwc(i) = zrho(i, k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content ENDIF ENDDO CALL fallice_velocity(klon, iwc(:), zt(:), zrho(:, k), pplay(:, k), ptconv(:, k), velo(:, k)) DO i = 1, klon IF (rneb(i, k)>0.0) THEN ! Initialization of zrain, zsnow and zprecip: zrain = 0. zsnow = 0. zprecip = 0. ! c_iso same init for isotopes. Externalisation? IF (zneb(i)==seuil_neb) THEN zprecip = 0.0 zsnow = 0.0 zrain = 0.0 ELSE IF (ptconv(i, k)) THEN ! if convective point zcl = cld_lc_con zct = 1. / cld_tau_con zexpo = cld_expo_con ffallv = ffallv_con ELSE zcl = cld_lc_lsc zct = 1. / cld_tau_lsc zexpo = cld_expo_lsc ffallv = ffallv_lsc ENDIF ! if vertical heterogeneity is taken into account, we use ! the "true" volume fraction instead of a modified ! surface fraction (which is larger and artificially ! reduces the in-cloud water). ! Liquid water quantity to remove: zchau (Sundqvist, 1978) ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2) !......................................................... IF ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) THEN ! if vertical heterogeneity is taken into account, we use ! the "true" volume fraction instead of a modified ! surface fraction (which is larger and artificially ! reduces the in-cloud water). effective_zneb = ctot_vol(i, k) ELSE effective_zneb = zneb(i) ENDIF IF (iflag_autoconversion == 2) THEN ! two-steps resolution with niter_lscp=1 sufficient ! we first treat the second term (with exponential) in an explicit way ! and then treat the first term (-q/tau) in an exact way zchau = zoliql(i) * (1. - exp(-dtime / REAL(niter_lscp) * zct & * (1. - exp(-(zoliql(i) / effective_zneb / zcl)**zexpo)))) ELSE ! old explicit resolution with subtimesteps zchau = zct * dtime / REAL(niter_lscp) * zoliql(i) & * (1.0 - EXP(-(zoliql(i) / effective_zneb / zcl)**zexpo)) ENDIF ! Ice water quantity to remove (Zender & Kiehl, 1997) ! dqice/dt=1/rho*d(rho*wice*qice)/dz !.................................... IF (iflag_autoconversion == 2) THEN ! exact resolution, niter_lscp=1 is sufficient but works only ! with iflag_vice=0 IF (zoliqi(i) > 0.) THEN zfroi = (zoliqi(i) - ((zoliqi(i)**(-dice_velo)) & + dice_velo * dtime / REAL(niter_lscp) * cice_velo / zdz(i) * ffallv)**(-1. / dice_velo)) ELSE zfroi = 0. ENDIF ELSE ! old explicit resolution with subtimesteps zfroi = dtime / REAL(niter_lscp) / zdz(i) * zoliqi(i) * velo(i, k) ENDIF zrain = MIN(MAX(zchau, 0.0), zoliql(i)) zsnow = MIN(MAX(zfroi, 0.0), zoliqi(i)) zprecip = MAX(zrain + zsnow, 0.0) ENDIF IF (iflag_autoconversion >= 1) THEN ! debugged version with phase conservation through the autoconversion process zoliql(i) = MAX(zoliql(i) - 1. * zrain, 0.0) zoliqi(i) = MAX(zoliqi(i) - 1. * zsnow, 0.0) zoliq(i) = MAX(zoliq(i) - zprecip, 0.0) ELSE ! bugged version with phase resetting zoliql(i) = MAX(zoliq(i) * (1. - zfice(i)) - 1. * zrain, 0.0) zoliqi(i) = MAX(zoliq(i) * zfice(i) - 1. * zsnow, 0.0) zoliq(i) = MAX(zoliq(i) - zprecip, 0.0) ENDIF ! c_iso: CALL isotope_conversion (for readibility) or duplicate radocond(i, k) = radocond(i, k) + zoliq(i) / REAL(niter_lscp + 1) radocondl(i, k) = radocondl(i, k) + zoliql(i) / REAL(niter_lscp + 1) radocondi(i, k) = radocondi(i, k) + zoliqi(i) / REAL(niter_lscp + 1) ENDIF ! rneb >0 ENDDO ! i = 1,klon ENDDO ! n = 1,niter ! Precipitation flux calculation DO i = 1, klon IF (iflag_evap_prec>=4) THEN ziflprev(i) = ziflcld(i) ELSE ziflprev(i) = zifl(i) * zneb(i) ENDIF IF (rneb(i, k) > 0.0) THEN ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux: ! If T<0C, liquid precip are converted into ice, which leads to an increase in ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly ! taken into account through a linearization of the equations and by approximating ! the liquid precipitation process with a "threshold" process. We assume that ! condensates are not modified during this operation. Liquid precipitation is ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases. ! Water vapor increases as well ! Note that compared to fisrtilp, we always assume iflag_bergeron=2 zqpreci(i) = (zcond(i) - zoliq(i)) * zfice(i) zqprecl(i) = (zcond(i) - zoliq(i)) * (1. - zfice(i)) zcp = RCPD * (1.0 + RVTMP2 * (zq(i) + zmqc(i) + zcond(i))) coef1 = rneb(i, k) * RLSTT / zcp * zdqsdT_raw(i) ! Computation of DT if all the liquid precip freezes DeltaT = RLMLT * zqprecl(i) / (zcp * (1. + coef1)) ! T should not exceed the freezing point ! that is Delta > RTT-zt(i) DeltaT = max(min(RTT - zt(i), DeltaT), 0.) zt(i) = zt(i) + DeltaT ! water vaporization due to temp. increase Deltaq = rneb(i, k) * zdqsdT_raw(i) * DeltaT ! we add this vaporized water to the total vapor and we remove it from the precip zq(i) = zq(i) + Deltaq ! The three "max" lines herebelow protect from rounding errors zcond(i) = max(zcond(i) - Deltaq, 0.) ! liquid precipitation converted to ice precip Deltaqprecl = -zcp / RLMLT * (1. + coef1) * DeltaT zqprecl(i) = max(zqprecl(i) + Deltaqprecl, 0.) ! iced water budget zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.) ! c_iso : mv here condensation of isotopes + redispatchage en precip IF (iflag_evap_prec>=4) THEN zrflcld(i) = zrflcld(i) + zqprecl(i) & * (paprs(i, k) - paprs(i, k + 1)) / (RG * dtime) ziflcld(i) = ziflcld(i) + zqpreci(i) & * (paprs(i, k) - paprs(i, k + 1)) / (RG * dtime) znebprecipcld(i) = rneb(i, k) zrfl(i) = zrflcld(i) + zrflclr(i) zifl(i) = ziflcld(i) + ziflclr(i) ELSE zrfl(i) = zrfl(i) + zqprecl(i) & * (paprs(i, k) - paprs(i, k + 1)) / (RG * dtime) zifl(i) = zifl(i) + zqpreci(i) & * (paprs(i, k) - paprs(i, k + 1)) / (RG * dtime) ENDIF ! c_iso : same for isotopes ENDIF ! rneb>0 ENDDO ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min ! if iflag_evap_prec>=4 IF (iflag_evap_prec>=4) THEN DO i = 1, klon IF ((zrflclr(i) + ziflclr(i)) > 0.) THEN znebprecipclr(i) = min(znebprecipclr(i), max(zrflclr(i) / & (MAX(znebprecipclr(i), seuil_neb) * rain_int_min), ziflclr(i) / (MAX(znebprecipclr(i), seuil_neb) * rain_int_min))) ELSE znebprecipclr(i) = 0.0 ENDIF IF ((zrflcld(i) + ziflcld(i)) > 0.) THEN znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i) / & (MAX(znebprecipcld(i), seuil_neb) * rain_int_min), ziflcld(i) / (MAX(znebprecipcld(i), seuil_neb) * rain_int_min))) ELSE znebprecipcld(i) = 0.0 ENDIF !IF ( ((1-zfice(i))*zoliq(i) .GT. 0.) .AND. (zt(i) .LE. 233.15) ) THEN !PRINT*,'WARNING LEA OLIQ A <-40°C ' !PRINT*,'zt,Tbef,oliq,oice,cldfraliq,icefrac,rneb',zt(i),Tbef(i),(1-zfice(i))*zoliq(i),zfice(i)*zoliq(i),cldfraliq(i,k),zfice(i),rneb(i,k) !ENDIF ENDDO ENDIF ENDIF ! ok_poprecip ! End of precipitation formation ! -------------------------------- ! Calculation of cloud condensates amount seen by radiative scheme !----------------------------------------------------------------- ! Calculation of the concentration of condensates seen by the radiation scheme ! for liquid phase, we take radocondl ! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true) ! we recalculate radocondi to account for contributions from the precipitation flux ! TODO: for the moment, we deactivate this option when ok_poprecip=.TRUE. IF ((ok_radocond_snow) .AND. (k < klev) .AND. (.NOT. ok_poprecip)) THEN ! for the solid phase (crystals + snowflakes) ! we recalculate radocondi by summing ! the ice content calculated in the mesh ! + the contribution of the non-evaporated snowfall ! from the overlying layer DO i = 1, klon IF (ziflprev(i) /= 0.0) THEN radocondi(i, k) = zoliq(i) * zfice(i) + zqpreci(i) + ziflprev(i) / zrho(i, k + 1) / velo(i, k + 1) ELSE radocondi(i, k) = zoliq(i) * zfice(i) + zqpreci(i) ENDIF radocond(i, k) = radocondl(i, k) + radocondi(i, k) ENDDO ENDIF ! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme DO i = 1, klon IF (radocond(i, k) > 0.) THEN radicefrac(i, k) = MIN(MAX(radocondi(i, k) / radocond(i, k), 0.), 1.) ENDIF ENDDO ! ---------------------------------------------------------------- ! P4> Wet scavenging ! ---------------------------------------------------------------- !Scavenging through nucleation in the layer DO i = 1, klon IF(zcond(i)>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)>0.0.AND.zprec_cond(i)>0.) THEN IF (temp(i, k) >= t_glace_min) THEN zalpha_tr = a_tr_sca(3) ELSE zalpha_tr = a_tr_sca(4) ENDIF zfrac_lessi = 1. - EXP(zalpha_tr * zprec_cond(i) / zneb(i)) frac_nucl(i, k) = 1. - zneb(i) * zfrac_lessi ! Nucleation with a factor of -1 instead of -0.5 zfrac_lessi = 1. - EXP(-zprec_cond(i) / zneb(i)) ENDIF ENDDO ! Scavenging through impaction in the underlying layer DO kk = k - 1, 1, -1 DO i = 1, klon IF (rneb(i, k)>0.0.AND.zprec_cond(i)>0.) THEN IF (temp(i, kk) >= t_glace_min) THEN zalpha_tr = a_tr_sca(1) ELSE zalpha_tr = a_tr_sca(2) ENDIF zfrac_lessi = 1. - EXP(zalpha_tr * zprec_cond(i) / zneb(i)) frac_impa(i, kk) = 1. - zneb(i) * zfrac_lessi ENDIF ENDDO ENDDO !--save some variables for ice supersaturation DO i = 1, klon ! for memory rneb_seri(i, k) = rneb(i, k) ! for diagnostics rnebclr(i, k) = 1.0 - rneb(i, k) - rnebss(i, k) qvc(i, k) = zqs(i) * rneb(i, k) qclr(i, k) = MAX(1.e-10, zq(i) - qvc(i, k) - qss(i, k)) !--added by OB because of pathologic cases with the lognormal qcld(i, k) = qvc(i, k) + zcond(i) ENDDO !q_sat CALL calc_qsat_ecmwf(klon, Tbef(:), qzero(:), pplay(:, k), RTT, 1, .FALSE., qsatl(:, k), zdqs(:)) CALL calc_qsat_ecmwf(klon, Tbef(:), qzero(:), pplay(:, k), RTT, 2, .FALSE., qsats(:, k), zdqs(:)) ! 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) d_ql(i, k) = (1 - zfice(i)) * zoliq(i) d_qi(i, k) = zfice(i) * zoliq(i) d_q(i, k) = zq(i) - qt(i, k) ! c_iso: same for isotopes d_t(i, k) = zt(i) - temp(i, k) ENDDO ENDDO ! 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 (ncoreczq>0) THEN WRITE(lunout, *)'WARNING : ZQ in LSCP ', ncoreczq, ' val < 1.e-15.' ENDIF END SUBROUTINE lscp !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ END MODULE lmdz_lscp