[4664] | 1 | MODULE lmdz_lscp |
---|
[3999] | 2 | |
---|
| 3 | IMPLICIT NONE |
---|
| 4 | |
---|
| 5 | CONTAINS |
---|
| 6 | |
---|
| 7 | !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
[4380] | 8 | SUBROUTINE lscp(klon,klev,dtime,missing_val, & |
---|
[4686] | 9 | paprs,pplay,temp,qt,ptconv,ratqs, & |
---|
[4380] | 10 | d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri, & |
---|
[4530] | 11 | pfraclr,pfracld, & |
---|
| 12 | radocond, radicefrac, rain, snow, & |
---|
[3999] | 13 | frac_impa, frac_nucl, beta, & |
---|
| 14 | prfl, psfl, rhcl, zqta, fraca, & |
---|
| 15 | ztv, zpspsk, ztla, zthl, iflag_cld_th, & |
---|
[4380] | 16 | iflag_ice_thermo, ok_ice_sursat, qsatl, qsats, & |
---|
[4830] | 17 | distcltop,temp_cltop, & |
---|
[4380] | 18 | qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & |
---|
[4651] | 19 | Tcontr, qcontr, qcontr2, fcontrN, fcontrP, & |
---|
[4803] | 20 | cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, & |
---|
[4830] | 21 | qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, & |
---|
| 22 | dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,& |
---|
[4819] | 23 | dqsmelt, dqsfreez) |
---|
[3999] | 24 | |
---|
| 25 | !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 26 | ! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD) |
---|
| 27 | ! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD) |
---|
| 28 | !-------------------------------------------------------------------------------- |
---|
| 29 | ! Date: 01/2021 |
---|
| 30 | !-------------------------------------------------------------------------------- |
---|
| 31 | ! Aim: Large Scale Clouds and Precipitation (LSCP) |
---|
[4226] | 32 | ! |
---|
[3999] | 33 | ! This code is a new version of the fisrtilp.F90 routine, which is itself a |
---|
[4072] | 34 | ! merge of 'first' (superrsaturation physics, P. LeVan K. Laval) |
---|
[3999] | 35 | ! and 'ilp' (il pleut, L. Li) |
---|
| 36 | ! |
---|
[4380] | 37 | ! Compared to the original fisrtilp code, lscp |
---|
[3999] | 38 | ! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.) |
---|
| 39 | ! -> consider always precipitation thermalisation (fl_cor_ebil>0) |
---|
| 40 | ! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T) |
---|
| 41 | ! -> option "oldbug" by JYG has been removed |
---|
| 42 | ! -> iflag_t_glace >0 always |
---|
| 43 | ! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always) |
---|
| 44 | ! -> rectangular distribution from L. Li no longer available |
---|
| 45 | ! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt) |
---|
| 46 | !-------------------------------------------------------------------------------- |
---|
| 47 | ! References: |
---|
| 48 | ! |
---|
[4412] | 49 | ! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2 |
---|
[3999] | 50 | ! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y |
---|
| 51 | ! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3 |
---|
[4412] | 52 | ! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379 |
---|
[3999] | 53 | ! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046 |
---|
[4412] | 54 | ! - Touzze-Peifert Ludo, PhD thesis, p117-124 |
---|
[3999] | 55 | ! ------------------------------------------------------------------------------- |
---|
| 56 | ! Code structure: |
---|
| 57 | ! |
---|
[4226] | 58 | ! P0> Thermalization of the precipitation coming from the overlying layer |
---|
[3999] | 59 | ! P1> Evaporation of the precipitation (falling from the k+1 level) |
---|
| 60 | ! P2> Cloud formation (at the k level) |
---|
[4412] | 61 | ! P2.A.1> With the PDFs, calculation of cloud properties using the inital |
---|
[3999] | 62 | ! values of T and Q |
---|
| 63 | ! P2.A.2> Coupling between condensed water and temperature |
---|
| 64 | ! P2.A.3> Calculation of final quantities associated with cloud formation |
---|
[4412] | 65 | ! P2.B> Release of Latent heat after cloud formation |
---|
[3999] | 66 | ! P3> Autoconversion to precipitation (k-level) |
---|
| 67 | ! P4> Wet scavenging |
---|
| 68 | !------------------------------------------------------------------------------ |
---|
| 69 | ! Some preliminary comments (JBM) : |
---|
| 70 | ! |
---|
| 71 | ! The cloud water that the radiation scheme sees is not the same that the cloud |
---|
| 72 | ! water used in the physics and the dynamics |
---|
| 73 | ! |
---|
[4412] | 74 | ! During the autoconversion to precipitation (P3 step), radocond (cloud water used |
---|
[3999] | 75 | ! by the radiation scheme) is calculated as an average of the water that remains |
---|
| 76 | ! in the cloud during the precipitation and not the water remaining at the end |
---|
| 77 | ! of the time step. The latter is used in the rest of the physics and advected |
---|
| 78 | ! by the dynamics. |
---|
| 79 | ! |
---|
| 80 | ! In summary: |
---|
| 81 | ! |
---|
| 82 | ! Radiation: |
---|
| 83 | ! xflwc(newmicro)+xfiwc(newmicro) = |
---|
[4412] | 84 | ! radocond=lwcon(nc)+iwcon(nc) |
---|
[3999] | 85 | ! |
---|
| 86 | ! Notetheless, be aware of: |
---|
| 87 | ! |
---|
[4412] | 88 | ! radocond .NE. ocond(nc) |
---|
[3999] | 89 | ! i.e.: |
---|
| 90 | ! lwcon(nc)+iwcon(nc) .NE. ocond(nc) |
---|
[4412] | 91 | ! but oliq+(ocond-oliq) .EQ. ocond |
---|
[3999] | 92 | ! (which is not trivial) |
---|
| 93 | !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
[4226] | 94 | ! |
---|
[4654] | 95 | |
---|
| 96 | ! USE de modules contenant des fonctions. |
---|
[4651] | 97 | USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc |
---|
[4664] | 98 | USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat |
---|
| 99 | USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top |
---|
[4380] | 100 | USE ice_sursat_mod, ONLY : ice_sursat |
---|
[4803] | 101 | USE lmdz_lscp_poprecip, ONLY : poprecip_evapsub, poprecip_postcld |
---|
[3999] | 102 | |
---|
[4664] | 103 | ! Use du module lmdz_lscp_ini contenant les constantes |
---|
[4666] | 104 | USE lmdz_lscp_ini, ONLY : prt_level, lunout |
---|
[4664] | 105 | USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min |
---|
| 106 | USE lmdz_lscp_ini, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc |
---|
| 107 | USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min |
---|
[4830] | 108 | USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con |
---|
[4664] | 109 | USE lmdz_lscp_ini, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo |
---|
| 110 | USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc |
---|
| 111 | USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG |
---|
[4869] | 112 | USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_fonte_lscp |
---|
[4380] | 113 | |
---|
[3999] | 114 | IMPLICIT NONE |
---|
| 115 | |
---|
| 116 | !=============================================================================== |
---|
[4380] | 117 | ! VARIABLES DECLARATION |
---|
[3999] | 118 | !=============================================================================== |
---|
| 119 | |
---|
| 120 | ! INPUT VARIABLES: |
---|
| 121 | !----------------- |
---|
| 122 | |
---|
[4380] | 123 | INTEGER, INTENT(IN) :: klon,klev ! number of horizontal grid points and vertical levels |
---|
[3999] | 124 | REAL, INTENT(IN) :: dtime ! time step [s] |
---|
[4380] | 125 | REAL, INTENT(IN) :: missing_val ! missing value for output |
---|
[4059] | 126 | |
---|
[3999] | 127 | REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa] |
---|
| 128 | REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa] |
---|
[4654] | 129 | REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K) |
---|
[4686] | 130 | REAL, DIMENSION(klon,klev), INTENT(IN) :: qt ! total specific humidity (in vapor phase in input) [kg/kg] |
---|
[3999] | 131 | INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds |
---|
| 132 | INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics |
---|
[4226] | 133 | ! CR: if iflag_ice_thermo=2, only convection is active |
---|
[4072] | 134 | LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated |
---|
[4059] | 135 | |
---|
[3999] | 136 | LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active |
---|
| 137 | |
---|
| 138 | !Inputs associated with thermal plumes |
---|
[4226] | 139 | |
---|
[3999] | 140 | REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! virtual potential temperature [K] |
---|
| 141 | REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg] |
---|
| 142 | REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-] |
---|
| 143 | REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk ! exner potential (p/100000)**(R/cp) |
---|
| 144 | REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! liquid temperature within thermals [K] |
---|
| 145 | |
---|
| 146 | ! INPUT/OUTPUT variables |
---|
| 147 | !------------------------ |
---|
[4562] | 148 | |
---|
| 149 | REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! liquid potential temperature [K] |
---|
[4686] | 150 | REAL, DIMENSION(klon,klev), INTENT(INOUT) :: ratqs ! function of pressure that sets the large-scale |
---|
[4059] | 151 | |
---|
[4380] | 152 | |
---|
[4059] | 153 | ! Input sursaturation en glace |
---|
| 154 | REAL, DIMENSION(klon,klev), INTENT(INOUT):: rneb_seri ! fraction nuageuse en memoire |
---|
[3999] | 155 | |
---|
| 156 | ! OUTPUT variables |
---|
| 157 | !----------------- |
---|
| 158 | |
---|
| 159 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K] |
---|
| 160 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg] |
---|
| 161 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg] |
---|
| 162 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg] |
---|
| 163 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-] |
---|
[4380] | 164 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneblsvol ! cloud fraction per unit volume [-] |
---|
[4530] | 165 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-] |
---|
| 166 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-] |
---|
[4412] | 167 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: radocond ! condensed water used in the radiation scheme [kg/kg] |
---|
[3999] | 168 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme |
---|
| 169 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-] |
---|
[4412] | 170 | REAL, DIMENSION(klon), INTENT(OUT) :: rain ! surface large-scale rainfall [kg/s/m2] |
---|
| 171 | REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2] |
---|
[4380] | 172 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl ! saturation specific humidity wrt liquid [kg/kg] |
---|
| 173 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsats ! saturation specific humidity wrt ice [kg/kg] |
---|
[3999] | 174 | REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2] |
---|
| 175 | REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2] |
---|
[4562] | 176 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: distcltop ! distance to cloud top [m] |
---|
[4639] | 177 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: temp_cltop ! temperature of cloud top [K] |
---|
[4686] | 178 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: beta ! conversion rate of condensed water |
---|
| 179 | |
---|
[3999] | 180 | ! fraction of aerosol scavenging through impaction and nucleation (for on-line) |
---|
| 181 | |
---|
| 182 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-] |
---|
| 183 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-] |
---|
| 184 | |
---|
[4380] | 185 | ! for supersaturation and contrails parameterisation |
---|
[3999] | 186 | |
---|
[4380] | 187 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qclr ! specific total water content in clear sky region [kg/kg] |
---|
| 188 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld ! specific total water content in cloudy region [kg/kg] |
---|
| 189 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qss ! specific total water content in supersat region [kg/kg] |
---|
| 190 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qvc ! specific vapor content in clouds [kg/kg] |
---|
| 191 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebclr ! mesh fraction of clear sky [-] |
---|
| 192 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebss ! mesh fraction of ISSR [-] |
---|
| 193 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_ss ! coefficient governing the ice nucleation RHi threshold [-] |
---|
[4425] | 194 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcontr ! threshold temperature for contrail formation [K] |
---|
| 195 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr ! threshold humidity for contrail formation [kg/kg] |
---|
| 196 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr2 ! // (2nd expression more consistent with LMDZ expression of q) |
---|
| 197 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrN ! fraction of grid favourable to non-persistent contrails |
---|
| 198 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrP ! fraction of grid favourable to persistent contrails |
---|
[4686] | 199 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals |
---|
| 200 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment |
---|
| 201 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals |
---|
| 202 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment |
---|
[3999] | 203 | |
---|
[4803] | 204 | ! for POPRECIP |
---|
[4686] | 205 | |
---|
[4830] | 206 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg] |
---|
| 207 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg] |
---|
[4819] | 208 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqreva !--rain tendendy due to evaporation [kg/kg/s] |
---|
| 209 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqssub !--snow tendency due to sublimation [kg/kg/s] |
---|
| 210 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s] |
---|
| 211 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s] |
---|
| 212 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s] |
---|
| 213 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s] |
---|
| 214 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsrim !--snow tendency due to riming [kg/kg/s] |
---|
| 215 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsmelt !--snow tendency due to melting [kg/kg/s] |
---|
| 216 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrmelt !--rain tendency due to melting [kg/kg/s] |
---|
| 217 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsfreez !--snow tendency due to freezing [kg/kg/s] |
---|
| 218 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrfreez !--rain tendency due to freezing [kg/kg/s] |
---|
[4686] | 219 | |
---|
[4803] | 220 | |
---|
[3999] | 221 | ! LOCAL VARIABLES: |
---|
| 222 | !---------------- |
---|
| 223 | |
---|
[4654] | 224 | REAL,DIMENSION(klon) :: qsl, qsi |
---|
| 225 | REAL :: zct, zcl,zexpo |
---|
| 226 | REAL, DIMENSION(klon,klev) :: ctot |
---|
| 227 | REAL, DIMENSION(klon,klev) :: ctot_vol |
---|
| 228 | REAL, DIMENSION(klon) :: zqs, zdqs |
---|
| 229 | REAL :: zdelta, zcor, zcvm5 |
---|
| 230 | REAL, DIMENSION(klon) :: zdqsdT_raw |
---|
| 231 | REAL, DIMENSION(klon) :: gammasat,dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T |
---|
| 232 | REAL, DIMENSION(klon) :: Tbef,qlbef,DT |
---|
| 233 | REAL :: num,denom |
---|
| 234 | REAL :: cste |
---|
| 235 | REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta |
---|
| 236 | REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2 |
---|
| 237 | REAL :: erf |
---|
| 238 | REAL, DIMENSION(klon,klev) :: icefrac_mpc |
---|
| 239 | REAL, DIMENSION(klon) :: qcloud, qincloud_mpc |
---|
| 240 | REAL, DIMENSION(klon) :: zrfl, zrfln |
---|
| 241 | REAL :: zqev, zqevt |
---|
| 242 | REAL, DIMENSION(klon) :: zifl, zifln, ziflprev |
---|
| 243 | REAL :: zqev0,zqevi, zqevti |
---|
| 244 | REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn |
---|
| 245 | REAL, DIMENSION(klon) :: zoliql, zoliqi |
---|
| 246 | REAL, DIMENSION(klon) :: zt |
---|
| 247 | REAL, DIMENSION(klon,klev) :: zrho |
---|
| 248 | REAL, DIMENSION(klon) :: zdz,iwc |
---|
| 249 | REAL :: zchau,zfroi |
---|
| 250 | REAL, DIMENSION(klon) :: zfice,zneb,znebprecip |
---|
| 251 | REAL :: zmelt,zrain,zsnow,zprecip |
---|
| 252 | REAL, DIMENSION(klon) :: dzfice |
---|
| 253 | REAL :: zsolid |
---|
| 254 | REAL, DIMENSION(klon) :: qtot, qzero |
---|
| 255 | REAL, DIMENSION(klon) :: dqsl,dqsi |
---|
| 256 | REAL :: smallestreal |
---|
[3999] | 257 | ! Variables for Bergeron process |
---|
[4654] | 258 | REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl |
---|
| 259 | REAL, DIMENSION(klon) :: zqpreci, zqprecl |
---|
[3999] | 260 | ! Variables precipitation energy conservation |
---|
[4654] | 261 | REAL, DIMENSION(klon) :: zmqc |
---|
| 262 | REAL :: zalpha_tr |
---|
| 263 | REAL :: zfrac_lessi |
---|
| 264 | REAL, DIMENSION(klon) :: zprec_cond |
---|
[4803] | 265 | REAL :: zmair |
---|
[4654] | 266 | REAL :: zcpair, zcpeau |
---|
| 267 | REAL, DIMENSION(klon) :: zlh_solid |
---|
[4803] | 268 | REAL, DIMENSION(klon) :: ztupnew |
---|
[4654] | 269 | REAL :: zm_solid ! for liquid -> solid conversion |
---|
| 270 | REAL, DIMENSION(klon) :: zrflclr, zrflcld |
---|
| 271 | REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld |
---|
| 272 | REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr |
---|
| 273 | REAL, DIMENSION(klon) :: ziflclr, ziflcld |
---|
| 274 | REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld |
---|
| 275 | REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb |
---|
| 276 | REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr |
---|
| 277 | REAL, DIMENSION(klon,klev) :: velo |
---|
| 278 | REAL :: vr, ffallv |
---|
| 279 | REAL :: qlmpc, qimpc, rnebmpc |
---|
| 280 | REAL, DIMENSION(klon,klev) :: radocondi, radocondl |
---|
| 281 | REAL :: effective_zneb |
---|
| 282 | REAL, DIMENSION(klon) :: distcltop1D, temp_cltop1D |
---|
[3999] | 283 | |
---|
[4562] | 284 | |
---|
[4654] | 285 | INTEGER i, k, n, kk, iter |
---|
| 286 | INTEGER, DIMENSION(klon) :: n_i |
---|
[4226] | 287 | INTEGER ncoreczq |
---|
[4654] | 288 | INTEGER, DIMENSION(klon,klev) :: mpc_bl_points |
---|
[4803] | 289 | LOGICAL iftop |
---|
[3999] | 290 | |
---|
[4654] | 291 | LOGICAL, DIMENSION(klon) :: lognormale |
---|
| 292 | LOGICAL, DIMENSION(klon) :: keepgoing |
---|
[3999] | 293 | |
---|
| 294 | CHARACTER (len = 20) :: modname = 'lscp' |
---|
| 295 | CHARACTER (len = 80) :: abort_message |
---|
| 296 | |
---|
| 297 | |
---|
| 298 | !=============================================================================== |
---|
| 299 | ! INITIALISATION |
---|
| 300 | !=============================================================================== |
---|
| 301 | |
---|
[4226] | 302 | ! Few initial checks |
---|
[3999] | 303 | |
---|
[4654] | 304 | |
---|
[3999] | 305 | IF (iflag_fisrtilp_qsat .LT. 0) THEN |
---|
[4226] | 306 | abort_message = 'lscp cannot be used with iflag_fisrtilp<0' |
---|
| 307 | CALL abort_physic(modname,abort_message,1) |
---|
[3999] | 308 | ENDIF |
---|
| 309 | |
---|
| 310 | ! Few initialisations |
---|
| 311 | |
---|
[4226] | 312 | znebprecip(:)=0.0 |
---|
[3999] | 313 | ctot_vol(1:klon,1:klev)=0.0 |
---|
| 314 | rneblsvol(1:klon,1:klev)=0.0 |
---|
[4226] | 315 | smallestreal=1.e-9 |
---|
| 316 | znebprecipclr(:)=0.0 |
---|
| 317 | znebprecipcld(:)=0.0 |
---|
[3999] | 318 | mpc_bl_points(:,:)=0 |
---|
| 319 | |
---|
| 320 | IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM' |
---|
| 321 | |
---|
| 322 | ! AA for 'safety' reasons |
---|
| 323 | zalpha_tr = 0. |
---|
| 324 | zfrac_lessi = 0. |
---|
[4380] | 325 | beta(:,:)= 0. |
---|
[3999] | 326 | |
---|
[4380] | 327 | ! Initialisation of variables: |
---|
| 328 | |
---|
[3999] | 329 | prfl(:,:) = 0.0 |
---|
| 330 | psfl(:,:) = 0.0 |
---|
| 331 | d_t(:,:) = 0.0 |
---|
| 332 | d_q(:,:) = 0.0 |
---|
| 333 | d_ql(:,:) = 0.0 |
---|
| 334 | d_qi(:,:) = 0.0 |
---|
| 335 | rneb(:,:) = 0.0 |
---|
[4530] | 336 | pfraclr(:,:)=0.0 |
---|
| 337 | pfracld(:,:)=0.0 |
---|
[4412] | 338 | radocond(:,:) = 0.0 |
---|
[3999] | 339 | radicefrac(:,:) = 0.0 |
---|
[4226] | 340 | frac_nucl(:,:) = 1.0 |
---|
| 341 | frac_impa(:,:) = 1.0 |
---|
[3999] | 342 | rain(:) = 0.0 |
---|
| 343 | snow(:) = 0.0 |
---|
[4114] | 344 | zoliq(:)=0.0 |
---|
| 345 | zfice(:)=0.0 |
---|
| 346 | dzfice(:)=0.0 |
---|
| 347 | zqprecl(:)=0.0 |
---|
| 348 | zqpreci(:)=0.0 |
---|
[3999] | 349 | zrfl(:) = 0.0 |
---|
| 350 | zifl(:) = 0.0 |
---|
[4114] | 351 | ziflprev(:)=0.0 |
---|
[3999] | 352 | zneb(:) = seuil_neb |
---|
[4226] | 353 | zrflclr(:) = 0.0 |
---|
| 354 | ziflclr(:) = 0.0 |
---|
| 355 | zrflcld(:) = 0.0 |
---|
| 356 | ziflcld(:) = 0.0 |
---|
| 357 | tot_zneb(:) = 0.0 |
---|
| 358 | tot_znebn(:) = 0.0 |
---|
| 359 | d_tot_zneb(:) = 0.0 |
---|
| 360 | qzero(:) = 0.0 |
---|
[4562] | 361 | distcltop1D(:)=0.0 |
---|
[4639] | 362 | temp_cltop1D(:) = 0.0 |
---|
[4803] | 363 | ztupnew(:)=0.0 |
---|
[4380] | 364 | !--ice supersaturation |
---|
[4226] | 365 | gamma_ss(:,:) = 1.0 |
---|
| 366 | qss(:,:) = 0.0 |
---|
| 367 | rnebss(:,:) = 0.0 |
---|
[4059] | 368 | Tcontr(:,:) = missing_val |
---|
| 369 | qcontr(:,:) = missing_val |
---|
| 370 | qcontr2(:,:) = missing_val |
---|
| 371 | fcontrN(:,:) = 0.0 |
---|
| 372 | fcontrP(:,:) = 0.0 |
---|
[4654] | 373 | distcltop(:,:)=0. |
---|
| 374 | temp_cltop(:,:)=0. |
---|
[4803] | 375 | !-- poprecip |
---|
[4830] | 376 | qraindiag(:,:)= 0. |
---|
| 377 | qsnowdiag(:,:)= 0. |
---|
[4819] | 378 | dqreva(:,:) = 0. |
---|
| 379 | dqrauto(:,:) = 0. |
---|
| 380 | dqrmelt(:,:) = 0. |
---|
| 381 | dqrfreez(:,:) = 0. |
---|
| 382 | dqrcol(:,:) = 0. |
---|
| 383 | dqssub(:,:) = 0. |
---|
| 384 | dqsauto(:,:) = 0. |
---|
| 385 | dqsrim(:,:) = 0. |
---|
| 386 | dqsagg(:,:) = 0. |
---|
| 387 | dqsfreez(:,:) = 0. |
---|
| 388 | dqsmelt(:,:) = 0. |
---|
[3999] | 389 | |
---|
[4803] | 390 | |
---|
| 391 | |
---|
| 392 | |
---|
[4392] | 393 | !c_iso: variable initialisation for iso |
---|
| 394 | |
---|
| 395 | |
---|
[3999] | 396 | !=============================================================================== |
---|
| 397 | ! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM |
---|
| 398 | !=============================================================================== |
---|
| 399 | |
---|
| 400 | ncoreczq=0 |
---|
| 401 | |
---|
| 402 | DO k = klev, 1, -1 |
---|
| 403 | |
---|
[4803] | 404 | IF (k.LE.klev-1) THEN |
---|
| 405 | iftop=.false. |
---|
| 406 | ELSE |
---|
| 407 | iftop=.true. |
---|
| 408 | ENDIF |
---|
| 409 | |
---|
[3999] | 410 | ! Initialisation temperature and specific humidity |
---|
| 411 | DO i = 1, klon |
---|
[4654] | 412 | zt(i)=temp(i,k) |
---|
[4686] | 413 | zq(i)=qt(i,k) |
---|
[4803] | 414 | IF (.not. iftop) THEN |
---|
| 415 | ztupnew(i)=temp(i,k+1)+d_t(i,k+1) |
---|
| 416 | ENDIF |
---|
[4392] | 417 | !c_iso init of iso |
---|
[3999] | 418 | ENDDO |
---|
[4803] | 419 | |
---|
| 420 | !================================================================ |
---|
| 421 | ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R) |
---|
| 422 | IF (ok_poprecip) THEN |
---|
| 423 | |
---|
| 424 | CALL poprecip_evapsub(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), & |
---|
| 425 | zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, & |
---|
| 426 | zrfl, zrflclr, zrflcld, & |
---|
| 427 | zifl, ziflclr, ziflcld, & |
---|
| 428 | dqreva(:,k),dqssub(:,k) & |
---|
| 429 | ) |
---|
| 430 | |
---|
| 431 | !================================================================ |
---|
| 432 | ELSE |
---|
| 433 | |
---|
[3999] | 434 | ! -------------------------------------------------------------------- |
---|
[4803] | 435 | ! P1> Thermalization of precipitation falling from the overlying layer |
---|
[3999] | 436 | ! -------------------------------------------------------------------- |
---|
[4412] | 437 | ! Computes air temperature variation due to enthalpy transported by |
---|
[3999] | 438 | ! precipitation. Precipitation is then thermalized with the air in the |
---|
| 439 | ! layer. |
---|
| 440 | ! The precipitation should remain thermalized throughout the different |
---|
[4412] | 441 | ! thermodynamical transformations. |
---|
| 442 | ! The corresponding water mass should |
---|
[3999] | 443 | ! be added when calculating the layer's enthalpy change with |
---|
| 444 | ! temperature |
---|
[4412] | 445 | ! See lmdzpedia page todoan |
---|
| 446 | ! todoan: check consistency with ice phase |
---|
| 447 | ! todoan: understand why several steps |
---|
[3999] | 448 | ! --------------------------------------------------------------------- |
---|
| 449 | |
---|
[4803] | 450 | IF (iftop) THEN |
---|
[3999] | 451 | |
---|
| 452 | DO i = 1, klon |
---|
[4803] | 453 | zmqc(i) = 0. |
---|
| 454 | ENDDO |
---|
| 455 | |
---|
| 456 | ELSE |
---|
| 457 | |
---|
| 458 | DO i = 1, klon |
---|
[3999] | 459 | |
---|
[4803] | 460 | zmair=(paprs(i,k)-paprs(i,k+1))/RG |
---|
[3999] | 461 | ! no condensed water so cp=cp(vapor+dry air) |
---|
| 462 | ! RVTMP2=rcpv/rcpd-1 |
---|
| 463 | zcpair=RCPD*(1.0+RVTMP2*zq(i)) |
---|
[4226] | 464 | zcpeau=RCPD*RVTMP2 |
---|
| 465 | |
---|
[3999] | 466 | ! zmqc: precipitation mass that has to be thermalized with |
---|
| 467 | ! layer's air so that precipitation at the ground has the |
---|
| 468 | ! same temperature as the lowermost layer |
---|
[4803] | 469 | zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair |
---|
[3999] | 470 | ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer |
---|
[4803] | 471 | zt(i) = ( ztupnew(i)*zmqc(i)*zcpeau + zcpair*zt(i) ) & |
---|
[3999] | 472 | / (zcpair + zmqc(i)*zcpeau) |
---|
[4226] | 473 | |
---|
[3999] | 474 | ENDDO |
---|
| 475 | |
---|
| 476 | ENDIF |
---|
| 477 | |
---|
| 478 | ! -------------------------------------------------------------------- |
---|
[4803] | 479 | ! P2> Precipitation evaporation/sublimation/melting |
---|
[3999] | 480 | ! -------------------------------------------------------------------- |
---|
[4803] | 481 | ! A part of the precipitation coming from above is evaporated/sublimated/melted. |
---|
[3999] | 482 | ! -------------------------------------------------------------------- |
---|
[4803] | 483 | |
---|
[4397] | 484 | IF (iflag_evap_prec.GE.1) THEN |
---|
[3999] | 485 | |
---|
[4072] | 486 | ! Calculation of saturation specific humidity |
---|
| 487 | ! depending on temperature: |
---|
[4380] | 488 | CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) |
---|
[4072] | 489 | ! wrt liquid water |
---|
[4380] | 490 | CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:)) |
---|
[4072] | 491 | ! wrt ice |
---|
[4380] | 492 | CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:)) |
---|
[3999] | 493 | |
---|
| 494 | DO i = 1, klon |
---|
[4226] | 495 | |
---|
[3999] | 496 | ! if precipitation |
---|
| 497 | IF (zrfl(i)+zifl(i).GT.0.) THEN |
---|
[4226] | 498 | |
---|
[4563] | 499 | ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4). |
---|
[4392] | 500 | ! c_iso: likely important to distinguish cs from neb isotope precipitation |
---|
| 501 | |
---|
[4563] | 502 | IF (iflag_evap_prec.GE.4) THEN |
---|
[4226] | 503 | zrfl(i) = zrflclr(i) |
---|
| 504 | zifl(i) = ziflclr(i) |
---|
| 505 | ENDIF |
---|
[3999] | 506 | |
---|
| 507 | IF (iflag_evap_prec.EQ.1) THEN |
---|
| 508 | znebprecip(i)=zneb(i) |
---|
| 509 | ELSE |
---|
| 510 | znebprecip(i)=MAX(zneb(i),znebprecip(i)) |
---|
| 511 | ENDIF |
---|
| 512 | |
---|
[4563] | 513 | IF (iflag_evap_prec.GT.4) THEN |
---|
| 514 | ! Max evaporation not to saturate the clear sky precip fraction |
---|
| 515 | ! i.e. the fraction where evaporation occurs |
---|
| 516 | zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i)) |
---|
| 517 | ELSEIF (iflag_evap_prec .EQ. 4) THEN |
---|
[4226] | 518 | ! Max evaporation not to saturate the whole mesh |
---|
[4563] | 519 | ! Pay attention -> lead to unrealistic and excessive evaporation |
---|
[4226] | 520 | zqev0 = MAX(0.0, zqs(i)-zq(i)) |
---|
| 521 | ELSE |
---|
| 522 | ! Max evap not to saturate the fraction below the cloud |
---|
| 523 | zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i)) |
---|
| 524 | ENDIF |
---|
[3999] | 525 | |
---|
| 526 | ! Evaporation of liquid precipitation coming from above |
---|
| 527 | ! dP/dz=beta*(1-q/qsat)*sqrt(P) |
---|
| 528 | ! formula from Sundquist 1988, Klemp & Wilhemson 1978 |
---|
[4563] | 529 | ! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4) |
---|
[3999] | 530 | |
---|
| 531 | IF (iflag_evap_prec.EQ.3) THEN |
---|
[4072] | 532 | zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) & |
---|
[3999] | 533 | *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) & |
---|
| 534 | *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
---|
[4563] | 535 | ELSE IF (iflag_evap_prec.GE.4) THEN |
---|
[4072] | 536 | zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) & |
---|
[3999] | 537 | *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) & |
---|
| 538 | *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
---|
| 539 | ELSE |
---|
[4072] | 540 | zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) & |
---|
[3999] | 541 | *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
---|
| 542 | ENDIF |
---|
| 543 | |
---|
[4226] | 544 | zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) & |
---|
| 545 | *RG*dtime/(paprs(i,k)-paprs(i,k+1)) |
---|
[3999] | 546 | |
---|
| 547 | ! sublimation of the solid precipitation coming from above |
---|
| 548 | IF (iflag_evap_prec.EQ.3) THEN |
---|
[4830] | 549 | zqevti = znebprecip(i)*coef_sub*(1.0-zq(i)/qsi(i)) & |
---|
[3999] | 550 | *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) & |
---|
| 551 | *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
---|
[4563] | 552 | ELSE IF (iflag_evap_prec.GE.4) THEN |
---|
[4830] | 553 | zqevti = znebprecipclr(i)*coef_sub*(1.0-zq(i)/qsi(i)) & |
---|
[3999] | 554 | *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) & |
---|
| 555 | *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
---|
| 556 | ELSE |
---|
[4830] | 557 | zqevti = 1.*coef_sub*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) & |
---|
[3999] | 558 | *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
---|
| 559 | ENDIF |
---|
[4226] | 560 | |
---|
[3999] | 561 | zqevti = MAX(0.0,MIN(zqevti,zifl(i))) & |
---|
[4226] | 562 | *RG*dtime/(paprs(i,k)-paprs(i,k+1)) |
---|
[3999] | 563 | |
---|
| 564 | ! A. JAM |
---|
| 565 | ! Evaporation limit: we ensure that the layer's fraction below |
---|
| 566 | ! the cloud or the whole mesh (depending on iflag_evap_prec) |
---|
| 567 | ! does not reach saturation. In this case, we |
---|
| 568 | ! redistribute zqev0 conserving the ratio liquid/ice |
---|
[4226] | 569 | |
---|
[3999] | 570 | IF (zqevt+zqevti.GT.zqev0) THEN |
---|
| 571 | zqev=zqev0*zqevt/(zqevt+zqevti) |
---|
| 572 | zqevi=zqev0*zqevti/(zqevt+zqevti) |
---|
| 573 | ELSE |
---|
| 574 | zqev=zqevt |
---|
| 575 | zqevi=zqevti |
---|
| 576 | ENDIF |
---|
| 577 | |
---|
| 578 | |
---|
[4392] | 579 | ! New solid and liquid precipitation fluxes after evap and sublimation |
---|
[3999] | 580 | zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & |
---|
| 581 | /RG/dtime) |
---|
| 582 | zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) & |
---|
| 583 | /RG/dtime) |
---|
| 584 | |
---|
[4392] | 585 | |
---|
[3999] | 586 | ! vapor, temperature, precip fluxes update |
---|
[4803] | 587 | ! vapor is updated after evaporation/sublimation (it is increased) |
---|
[3999] | 588 | zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & |
---|
| 589 | * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
---|
[4803] | 590 | ! zmqc is the total condensed water in the precip flux (it is decreased) |
---|
[3999] | 591 | zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & |
---|
| 592 | * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
---|
[4803] | 593 | ! air and precip temperature (i.e., gridbox temperature) |
---|
| 594 | ! is updated due to latent heat cooling |
---|
[3999] | 595 | zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & |
---|
| 596 | * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
---|
| 597 | * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) & |
---|
| 598 | + (zifln(i)-zifl(i)) & |
---|
| 599 | * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
---|
| 600 | * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) |
---|
| 601 | |
---|
| 602 | ! New values of liquid and solid precipitation |
---|
| 603 | zrfl(i) = zrfln(i) |
---|
| 604 | zifl(i) = zifln(i) |
---|
| 605 | |
---|
[4392] | 606 | ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout) |
---|
| 607 | ! due to evap + sublim |
---|
| 608 | |
---|
| 609 | |
---|
[4563] | 610 | IF (iflag_evap_prec.GE.4) THEN |
---|
[4226] | 611 | zrflclr(i) = zrfl(i) |
---|
| 612 | ziflclr(i) = zifl(i) |
---|
| 613 | IF(zrflclr(i) + ziflclr(i).LE.0) THEN |
---|
| 614 | znebprecipclr(i) = 0.0 |
---|
| 615 | ENDIF |
---|
| 616 | zrfl(i) = zrflclr(i) + zrflcld(i) |
---|
| 617 | zifl(i) = ziflclr(i) + ziflcld(i) |
---|
| 618 | ENDIF |
---|
[3999] | 619 | |
---|
[4392] | 620 | ! c_iso duplicate for isotopes or loop on isotopes |
---|
[3999] | 621 | |
---|
[4392] | 622 | ! Melting: |
---|
[4412] | 623 | zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG |
---|
[3999] | 624 | ! precip fraction that is melted |
---|
| 625 | zmelt = MIN(MAX(zmelt,0.),1.) |
---|
[4392] | 626 | |
---|
[4412] | 627 | ! update of rainfall and snowfall due to melting |
---|
[4563] | 628 | IF (iflag_evap_prec.GE.4) THEN |
---|
[4226] | 629 | zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i) |
---|
| 630 | zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i) |
---|
| 631 | zrfl(i)=zrflclr(i)+zrflcld(i) |
---|
[4869] | 632 | IF (ok_bug_fonte_lscp) THEN |
---|
[4412] | 633 | ziflclr(i)=ziflclr(i)*(1.-zmelt) |
---|
| 634 | ziflcld(i)=ziflcld(i)*(1.-zmelt) |
---|
| 635 | zifl(i)=ziflclr(i)+ziflcld(i) |
---|
[4869] | 636 | ENDIF |
---|
[4226] | 637 | ELSE |
---|
| 638 | zrfl(i)=zrfl(i)+zmelt*zifl(i) |
---|
[4869] | 639 | IF (ok_bug_fonte_lscp) THEN |
---|
[4412] | 640 | zifl(i)=zifl(i)*(1.-zmelt) |
---|
[4869] | 641 | ENDIF |
---|
[4226] | 642 | ENDIF |
---|
[4412] | 643 | |
---|
| 644 | |
---|
[4392] | 645 | ! c_iso: melting of isotopic precipi with zmelt (no fractionation) |
---|
[3999] | 646 | |
---|
[4803] | 647 | ! Latent heat of melting because of precipitation melting |
---|
| 648 | ! NB: the air + precip temperature is simultaneously updated |
---|
[3999] | 649 | zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) & |
---|
| 650 | *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) |
---|
[4869] | 651 | |
---|
| 652 | IF (.NOT. ok_bug_fonte_lscp) THEN |
---|
| 653 | IF (iflag_evap_prec.GE.4) THEN |
---|
| 654 | ziflclr(i)=ziflclr(i)*(1.-zmelt) |
---|
| 655 | ziflcld(i)=ziflcld(i)*(1.-zmelt) |
---|
| 656 | zifl(i)=ziflclr(i)+ziflcld(i) |
---|
| 657 | ELSE |
---|
| 658 | zifl(i)=zifl(i)*(1.-zmelt) |
---|
| 659 | ENDIF |
---|
| 660 | ENDIF |
---|
[3999] | 661 | |
---|
| 662 | ELSE |
---|
| 663 | ! if no precip, we reinitialize the cloud fraction used for the precip to 0 |
---|
| 664 | znebprecip(i)=0. |
---|
| 665 | |
---|
| 666 | ENDIF ! (zrfl(i)+zifl(i).GT.0.) |
---|
| 667 | |
---|
[4803] | 668 | ENDDO ! loop on klon |
---|
[4226] | 669 | |
---|
[3999] | 670 | ENDIF ! (iflag_evap_prec>=1) |
---|
| 671 | |
---|
[4803] | 672 | ENDIF ! (ok_poprecip) |
---|
| 673 | |
---|
[3999] | 674 | ! -------------------------------------------------------------------- |
---|
| 675 | ! End precip evaporation |
---|
| 676 | ! -------------------------------------------------------------------- |
---|
| 677 | |
---|
| 678 | ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter |
---|
| 679 | !------------------------------------------------------- |
---|
[4226] | 680 | |
---|
[4072] | 681 | qtot(:)=zq(:)+zmqc(:) |
---|
[4380] | 682 | CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) |
---|
[4072] | 683 | DO i = 1, klon |
---|
[3999] | 684 | zdelta = MAX(0.,SIGN(1.,RTT-zt(i))) |
---|
[4059] | 685 | zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta) |
---|
[3999] | 686 | IF (zq(i) .LT. 1.e-15) THEN |
---|
| 687 | ncoreczq=ncoreczq+1 |
---|
| 688 | zq(i)=1.e-15 |
---|
| 689 | ENDIF |
---|
[4392] | 690 | ! c_iso: do something similar for isotopes |
---|
[3999] | 691 | |
---|
[4072] | 692 | ENDDO |
---|
[3999] | 693 | |
---|
| 694 | ! -------------------------------------------------------------------- |
---|
| 695 | ! P2> Cloud formation |
---|
| 696 | !--------------------------------------------------------------------- |
---|
| 697 | ! |
---|
| 698 | ! Unlike fisrtilp, we always assume a 'fractional cloud' approach |
---|
| 699 | ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution |
---|
| 700 | ! is prescribed and depends on large scale variables and boundary layer |
---|
| 701 | ! properties) |
---|
| 702 | ! The decrease in condensed part due tu latent heating is taken into |
---|
| 703 | ! account |
---|
| 704 | ! ------------------------------------------------------------------- |
---|
| 705 | |
---|
| 706 | ! P2.1> With the PDFs (log-normal, bigaussian) |
---|
[4424] | 707 | ! cloud properties calculation with the initial values of t and q |
---|
[3999] | 708 | ! ---------------------------------------------------------------- |
---|
| 709 | |
---|
| 710 | ! initialise gammasat and qincloud_mpc |
---|
| 711 | gammasat(:)=1. |
---|
| 712 | qincloud_mpc(:)=0. |
---|
| 713 | |
---|
| 714 | IF (iflag_cld_th.GE.5) THEN |
---|
| 715 | |
---|
| 716 | IF (iflag_mpc_bl .LT. 1) THEN |
---|
| 717 | |
---|
| 718 | IF (iflag_cloudth_vert.LE.2) THEN |
---|
| 719 | |
---|
| 720 | CALL cloudth(klon,klev,k,ztv, & |
---|
| 721 | zq,zqta,fraca, & |
---|
| 722 | qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & |
---|
[4654] | 723 | ratqs,zqs,temp, & |
---|
[4651] | 724 | cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) |
---|
[3999] | 725 | |
---|
[4651] | 726 | |
---|
[3999] | 727 | ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN |
---|
| 728 | |
---|
| 729 | CALL cloudth_v3(klon,klev,k,ztv, & |
---|
| 730 | zq,zqta,fraca, & |
---|
| 731 | qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & |
---|
[4654] | 732 | ratqs,zqs,temp, & |
---|
[4651] | 733 | cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) |
---|
[3999] | 734 | |
---|
[4226] | 735 | !Jean Jouhaud's version, Decembre 2018 |
---|
| 736 | !------------------------------------- |
---|
[3999] | 737 | |
---|
| 738 | ELSEIF (iflag_cloudth_vert.EQ.6) THEN |
---|
| 739 | |
---|
| 740 | CALL cloudth_v6(klon,klev,k,ztv, & |
---|
| 741 | zq,zqta,fraca, & |
---|
| 742 | qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & |
---|
[4654] | 743 | ratqs,zqs,temp, & |
---|
[4651] | 744 | cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) |
---|
[3999] | 745 | |
---|
| 746 | ENDIF |
---|
| 747 | |
---|
| 748 | ELSE |
---|
| 749 | ! cloudth_v3 + cold microphysical considerations |
---|
| 750 | ! + stationary mixed-phase cloud model |
---|
| 751 | |
---|
[4380] | 752 | CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, & |
---|
[4654] | 753 | temp,ztv,zq,zqta,fraca, zpspsk, & |
---|
[3999] | 754 | paprs,pplay,ztla,zthl,ratqs,zqs,psfl, & |
---|
[4651] | 755 | qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol, & |
---|
| 756 | cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) |
---|
[3999] | 757 | |
---|
| 758 | ENDIF ! iflag_mpc_bl |
---|
| 759 | |
---|
| 760 | DO i=1,klon |
---|
| 761 | rneb(i,k)=ctot(i,k) |
---|
| 762 | rneblsvol(i,k)=ctot_vol(i,k) |
---|
| 763 | zqn(i)=qcloud(i) |
---|
| 764 | ENDDO |
---|
| 765 | |
---|
| 766 | ENDIF |
---|
| 767 | |
---|
| 768 | IF (iflag_cld_th .LE. 4) THEN |
---|
| 769 | |
---|
| 770 | ! lognormal |
---|
| 771 | lognormale = .TRUE. |
---|
| 772 | |
---|
| 773 | ELSEIF (iflag_cld_th .GE. 6) THEN |
---|
| 774 | |
---|
| 775 | ! lognormal distribution when no thermals |
---|
| 776 | lognormale = fraca(:,k) < 1e-10 |
---|
| 777 | |
---|
| 778 | ELSE |
---|
[4226] | 779 | ! When iflag_cld_th=5, we always assume |
---|
[3999] | 780 | ! bi-gaussian distribution |
---|
| 781 | lognormale = .FALSE. |
---|
| 782 | |
---|
| 783 | ENDIF |
---|
| 784 | |
---|
| 785 | DT(:) = 0. |
---|
| 786 | n_i(:)=0 |
---|
| 787 | Tbef(:)=zt(:) |
---|
| 788 | qlbef(:)=0. |
---|
| 789 | |
---|
| 790 | ! Treatment of non-boundary layer clouds (lognormale) |
---|
| 791 | ! condensation with qsat(T) variation (adaptation) |
---|
[4424] | 792 | ! Iterative resolution to converge towards qsat |
---|
| 793 | ! with update of temperature, ice fraction and qsat at |
---|
| 794 | ! each iteration |
---|
[4226] | 795 | |
---|
[4424] | 796 | ! todoan -> sensitivity to iflag_fisrtilp_qsat |
---|
[3999] | 797 | DO iter=1,iflag_fisrtilp_qsat+1 |
---|
[4226] | 798 | |
---|
[3999] | 799 | DO i=1,klon |
---|
| 800 | |
---|
[4424] | 801 | ! keepgoing = .true. while convergence is not satisfied |
---|
| 802 | keepgoing(i)=ABS(DT(i)).GT.DDT0 |
---|
[4226] | 803 | |
---|
[4424] | 804 | IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN |
---|
[4226] | 805 | |
---|
[3999] | 806 | ! if not convergence: |
---|
| 807 | |
---|
| 808 | ! P2.2.1> cloud fraction and condensed water mass calculation |
---|
| 809 | ! Calculated variables: |
---|
| 810 | ! rneb : cloud fraction |
---|
| 811 | ! zqn : total water within the cloud |
---|
| 812 | ! zcond: mean condensed water within the mesh |
---|
| 813 | ! rhcl: clear-sky relative humidity |
---|
| 814 | !--------------------------------------------------------------- |
---|
| 815 | |
---|
[4424] | 816 | ! new temperature that only serves in the iteration process: |
---|
[3999] | 817 | Tbef(i)=Tbef(i)+DT(i) |
---|
| 818 | |
---|
| 819 | ! Rneb, qzn and zcond for lognormal PDFs |
---|
[4072] | 820 | qtot(i)=zq(i)+zmqc(i) |
---|
[4226] | 821 | |
---|
[4072] | 822 | ENDIF |
---|
[3999] | 823 | |
---|
[4072] | 824 | ENDDO |
---|
| 825 | |
---|
[4380] | 826 | ! Calculation of saturation specific humidity and ice fraction |
---|
[4226] | 827 | CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) |
---|
| 828 | CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:)) |
---|
[4072] | 829 | ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs |
---|
| 830 | zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:) |
---|
[4562] | 831 | ! cloud phase determination |
---|
| 832 | IF (iflag_t_glace.GE.4) THEN |
---|
| 833 | ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top |
---|
[4654] | 834 | CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D) |
---|
[4562] | 835 | ENDIF |
---|
[4639] | 836 | CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),temp_cltop1D(:),zfice(:),dzfice(:)) |
---|
[4072] | 837 | |
---|
[4424] | 838 | DO i=1,klon !todoan : check if loop in i is needed |
---|
[4072] | 839 | |
---|
[4424] | 840 | IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN |
---|
[4072] | 841 | |
---|
[3999] | 842 | zpdf_sig(i)=ratqs(i,k)*zq(i) |
---|
| 843 | zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2)) |
---|
| 844 | zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i))) |
---|
| 845 | zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.)) |
---|
| 846 | zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.)) |
---|
| 847 | zpdf_e1(i)=zpdf_a(i)-zpdf_b(i) |
---|
| 848 | zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i)) |
---|
| 849 | zpdf_e1(i)=1.-erf(zpdf_e1(i)) |
---|
| 850 | zpdf_e2(i)=zpdf_a(i)+zpdf_b(i) |
---|
| 851 | zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i)) |
---|
| 852 | zpdf_e2(i)=1.-erf(zpdf_e2(i)) |
---|
[4226] | 853 | |
---|
[4062] | 854 | IF ((.NOT.ok_ice_sursat).OR.(Tbef(i).GT.t_glace_min)) THEN |
---|
[4059] | 855 | |
---|
| 856 | IF (zpdf_e1(i).LT.1.e-10) THEN |
---|
| 857 | rneb(i,k)=0. |
---|
| 858 | zqn(i)=gammasat(i)*zqs(i) |
---|
| 859 | ELSE |
---|
| 860 | rneb(i,k)=0.5*zpdf_e1(i) |
---|
| 861 | zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i) |
---|
| 862 | ENDIF |
---|
| 863 | |
---|
[4072] | 864 | rnebss(i,k)=0.0 !--added by OB (needed because of the convergence loop on time) |
---|
[4059] | 865 | fcontrN(i,k)=0.0 !--idem |
---|
| 866 | fcontrP(i,k)=0.0 !--idem |
---|
| 867 | qss(i,k)=0.0 !--idem |
---|
[4226] | 868 | |
---|
[4424] | 869 | ELSE ! in case of ice supersaturation by Audran |
---|
[4226] | 870 | |
---|
[4059] | 871 | !------------------------------------ |
---|
[4072] | 872 | ! ICE SUPERSATURATION |
---|
[4059] | 873 | !------------------------------------ |
---|
[3999] | 874 | |
---|
[4654] | 875 | CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(i,k), zq(i), & |
---|
[4226] | 876 | gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), & |
---|
| 877 | rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), & |
---|
[4059] | 878 | Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) ) |
---|
| 879 | |
---|
| 880 | ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min)) |
---|
| 881 | |
---|
[3999] | 882 | ! If vertical heterogeneity, change fraction by volume as well |
---|
| 883 | IF (iflag_cloudth_vert.GE.3) THEN |
---|
| 884 | ctot_vol(i,k)=rneb(i,k) |
---|
| 885 | rneblsvol(i,k)=ctot_vol(i,k) |
---|
| 886 | ENDIF |
---|
| 887 | |
---|
| 888 | |
---|
[4072] | 889 | ! P2.2.2> Approximative calculation of temperature variation DT |
---|
| 890 | ! due to condensation. |
---|
| 891 | ! Calculated variables: |
---|
| 892 | ! dT : temperature change due to condensation |
---|
| 893 | !--------------------------------------------------------------- |
---|
[3999] | 894 | |
---|
| 895 | |
---|
| 896 | IF (zfice(i).LT.1) THEN |
---|
| 897 | cste=RLVTT |
---|
| 898 | ELSE |
---|
| 899 | cste=RLSTT |
---|
| 900 | ENDIF |
---|
| 901 | |
---|
| 902 | qlbef(i)=max(0.,zqn(i)-zqs(i)) |
---|
| 903 | num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT & |
---|
[4226] | 904 | +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i) |
---|
[3999] | 905 | denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) & |
---|
| 906 | -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k) & |
---|
[4226] | 907 | *qlbef(i)*dzfice(i) |
---|
[4424] | 908 | ! here we update a provisory temperature variable that only serves in the iteration |
---|
| 909 | ! process |
---|
[3999] | 910 | DT(i)=num/denom |
---|
| 911 | n_i(i)=n_i(i)+1 |
---|
| 912 | |
---|
[4424] | 913 | ENDIF ! end keepgoing |
---|
[3999] | 914 | |
---|
| 915 | ENDDO ! end loop on i |
---|
| 916 | |
---|
| 917 | ENDDO ! iter=1,iflag_fisrtilp_qsat+1 |
---|
| 918 | |
---|
| 919 | ! P2.3> Final quantities calculation |
---|
| 920 | ! Calculated variables: |
---|
| 921 | ! rneb : cloud fraction |
---|
| 922 | ! zcond: mean condensed water in the mesh |
---|
| 923 | ! zqn : mean water vapor in the mesh |
---|
[4562] | 924 | ! zfice: ice fraction in clouds |
---|
[3999] | 925 | ! zt : temperature |
---|
| 926 | ! rhcl : clear-sky relative humidity |
---|
| 927 | ! ---------------------------------------------------------------- |
---|
| 928 | |
---|
| 929 | |
---|
[4562] | 930 | ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top |
---|
| 931 | IF (iflag_t_glace.GE.4) THEN |
---|
[4654] | 932 | CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D) |
---|
[4562] | 933 | distcltop(:,k)=distcltop1D(:) |
---|
[4639] | 934 | temp_cltop(:,k)=temp_cltop1D(:) |
---|
[4562] | 935 | ENDIF |
---|
| 936 | |
---|
[3999] | 937 | ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs) |
---|
[4803] | 938 | CALL icefrac_lscp(klon,zt,iflag_ice_thermo,distcltop1D,temp_cltop1D,zfice, dzfice) |
---|
[3999] | 939 | |
---|
[4562] | 940 | |
---|
| 941 | ! Water vapor update, Phase determination and subsequent latent heat exchange |
---|
[3999] | 942 | DO i=1, klon |
---|
| 943 | |
---|
| 944 | IF (mpc_bl_points(i,k) .GT. 0) THEN |
---|
| 945 | zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k) |
---|
| 946 | ! following line is very strange and probably wrong |
---|
| 947 | rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i) |
---|
| 948 | ! water vapor update and partition function if thermals |
---|
| 949 | zq(i) = zq(i) - zcond(i) |
---|
| 950 | zfice(i)=icefrac_mpc(i,k) |
---|
| 951 | |
---|
| 952 | ELSE |
---|
| 953 | |
---|
| 954 | ! Checks on rneb, rhcl and zqn |
---|
| 955 | IF (rneb(i,k) .LE. 0.0) THEN |
---|
| 956 | zqn(i) = 0.0 |
---|
| 957 | rneb(i,k) = 0.0 |
---|
| 958 | zcond(i) = 0.0 |
---|
| 959 | rhcl(i,k)=zq(i)/zqs(i) |
---|
| 960 | ELSE IF (rneb(i,k) .GE. 1.0) THEN |
---|
| 961 | zqn(i) = zq(i) |
---|
| 962 | rneb(i,k) = 1.0 |
---|
[4818] | 963 | zcond(i) = MAX(0.0,zqn(i)-zqs(i)) |
---|
[3999] | 964 | rhcl(i,k)=1.0 |
---|
| 965 | ELSE |
---|
[4818] | 966 | zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k) |
---|
[3999] | 967 | ! following line is very strange and probably wrong: |
---|
| 968 | rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i) |
---|
| 969 | ENDIF |
---|
| 970 | |
---|
[4226] | 971 | ! water vapor update |
---|
| 972 | zq(i) = zq(i) - zcond(i) |
---|
[3999] | 973 | |
---|
| 974 | ENDIF |
---|
| 975 | |
---|
[4392] | 976 | ! c_iso : routine that computes in-cloud supersaturation |
---|
| 977 | ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input) |
---|
| 978 | |
---|
[3999] | 979 | ! temperature update due to phase change |
---|
| 980 | zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & |
---|
| 981 | & * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) & |
---|
| 982 | +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) |
---|
| 983 | ENDDO |
---|
| 984 | |
---|
| 985 | ! If vertical heterogeneity, change volume fraction |
---|
| 986 | IF (iflag_cloudth_vert .GE. 3) THEN |
---|
| 987 | ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.) |
---|
| 988 | rneblsvol(1:klon,k)=ctot_vol(1:klon,k) |
---|
| 989 | ENDIF |
---|
| 990 | |
---|
| 991 | ! ---------------------------------------------------------------- |
---|
[4114] | 992 | ! P3> Precipitation formation |
---|
[3999] | 993 | ! ---------------------------------------------------------------- |
---|
[4226] | 994 | |
---|
[4803] | 995 | !================================================================ |
---|
| 996 | IF (ok_poprecip) THEN |
---|
| 997 | |
---|
| 998 | DO i = 1, klon |
---|
| 999 | zoliql(i) = zcond(i) * ( 1. - zfice(i) ) |
---|
| 1000 | zoliqi(i) = zcond(i) * zfice(i) |
---|
| 1001 | ENDDO |
---|
| 1002 | |
---|
| 1003 | CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), & |
---|
| 1004 | ctot_vol(:,k), ptconv(:,k), & |
---|
| 1005 | zt, zq, zoliql, zoliqi, zfice, & |
---|
| 1006 | rneb(:,k), znebprecipclr, znebprecipcld, & |
---|
| 1007 | zrfl, zrflclr, zrflcld, & |
---|
| 1008 | zifl, ziflclr, ziflcld, & |
---|
[4830] | 1009 | qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), & |
---|
[4819] | 1010 | dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), & |
---|
| 1011 | dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), & |
---|
| 1012 | dqsmelt(:,k), dqsfreez(:,k) & |
---|
[4803] | 1013 | ) |
---|
| 1014 | |
---|
| 1015 | DO i = 1, klon |
---|
| 1016 | zoliq(i) = zoliql(i) + zoliqi(i) |
---|
| 1017 | IF ( zoliq(i) .GT. 0. ) THEN |
---|
| 1018 | zfice(i) = zoliqi(i)/zoliq(i) |
---|
| 1019 | ELSE |
---|
| 1020 | zfice(i) = 0.0 |
---|
| 1021 | ENDIF |
---|
[4819] | 1022 | |
---|
| 1023 | ! calculation of specific content of condensates seen by radiative scheme |
---|
| 1024 | IF (ok_radocond_snow) THEN |
---|
| 1025 | radocond(i,k) = zoliq(i) |
---|
| 1026 | radocondl(i,k)= radocond(i,k)*(1.-zfice(i)) |
---|
[4830] | 1027 | radocondi(i,k)= radocond(i,k)*zfice(i)+qsnowdiag(i,k) |
---|
[4819] | 1028 | ELSE |
---|
| 1029 | radocond(i,k) = zoliq(i) |
---|
| 1030 | radocondl(i,k)= radocond(i,k)*(1.-zfice(i)) |
---|
| 1031 | radocondi(i,k)= radocond(i,k)*zfice(i) |
---|
| 1032 | ENDIF |
---|
[4803] | 1033 | ENDDO |
---|
| 1034 | |
---|
| 1035 | !================================================================ |
---|
| 1036 | ELSE |
---|
| 1037 | |
---|
[3999] | 1038 | ! LTP: |
---|
[4563] | 1039 | IF (iflag_evap_prec .GE. 4) THEN |
---|
[3999] | 1040 | |
---|
| 1041 | !Partitionning between precipitation coming from clouds and that coming from CS |
---|
| 1042 | |
---|
| 1043 | !0) Calculate tot_zneb, total cloud fraction above the cloud |
---|
| 1044 | !assuming a maximum-random overlap (voir Jakob and Klein, 2000) |
---|
| 1045 | |
---|
| 1046 | DO i=1, klon |
---|
[4424] | 1047 | tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) & |
---|
| 1048 | /(1.-min(zneb(i),1.-smallestreal)) |
---|
[3999] | 1049 | d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i) |
---|
| 1050 | tot_zneb(i) = tot_znebn(i) |
---|
| 1051 | |
---|
| 1052 | |
---|
| 1053 | !1) Cloudy to clear air |
---|
| 1054 | d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i)) |
---|
[4424] | 1055 | IF (znebprecipcld(i) .GT. 0.) THEN |
---|
[3999] | 1056 | d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i) |
---|
| 1057 | d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i) |
---|
| 1058 | ELSE |
---|
| 1059 | d_zrfl_cld_clr(i) = 0. |
---|
| 1060 | d_zifl_cld_clr(i) = 0. |
---|
| 1061 | ENDIF |
---|
| 1062 | |
---|
| 1063 | !2) Clear to cloudy air |
---|
[4226] | 1064 | d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i))) |
---|
[3999] | 1065 | IF (znebprecipclr(i) .GT. 0) THEN |
---|
| 1066 | d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i) |
---|
| 1067 | d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i) |
---|
| 1068 | ELSE |
---|
| 1069 | d_zrfl_clr_cld(i) = 0. |
---|
| 1070 | d_zifl_clr_cld(i) = 0. |
---|
| 1071 | ENDIF |
---|
| 1072 | |
---|
| 1073 | !Update variables |
---|
[4226] | 1074 | znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i) |
---|
[3999] | 1075 | znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i) |
---|
| 1076 | zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i) |
---|
| 1077 | ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i) |
---|
| 1078 | zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i) |
---|
| 1079 | ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i) |
---|
| 1080 | |
---|
[4392] | 1081 | ! c_iso do the same thing for isotopes precip |
---|
[3999] | 1082 | ENDDO |
---|
| 1083 | ENDIF |
---|
| 1084 | |
---|
[4559] | 1085 | |
---|
| 1086 | ! Autoconversion |
---|
| 1087 | ! ------------------------------------------------------------------------------- |
---|
| 1088 | |
---|
| 1089 | |
---|
[4412] | 1090 | ! Initialisation of zoliq and radocond variables |
---|
[3999] | 1091 | |
---|
| 1092 | DO i = 1, klon |
---|
| 1093 | zoliq(i) = zcond(i) |
---|
[4072] | 1094 | zoliqi(i)= zoliq(i)*zfice(i) |
---|
| 1095 | zoliql(i)= zoliq(i)*(1.-zfice(i)) |
---|
[4392] | 1096 | ! c_iso : initialisation of zoliq* also for isotopes |
---|
[4114] | 1097 | zrho(i,k) = pplay(i,k) / zt(i) / RD |
---|
| 1098 | zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG) |
---|
[4072] | 1099 | iwc(i) = 0. |
---|
| 1100 | zneb(i) = MAX(rneb(i,k),seuil_neb) |
---|
[4559] | 1101 | radocond(i,k) = zoliq(i)/REAL(niter_lscp+1) |
---|
| 1102 | radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1) |
---|
| 1103 | radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1) |
---|
[3999] | 1104 | ENDDO |
---|
| 1105 | |
---|
[4559] | 1106 | |
---|
| 1107 | DO n = 1, niter_lscp |
---|
[3999] | 1108 | |
---|
| 1109 | DO i=1,klon |
---|
| 1110 | IF (rneb(i,k).GT.0.0) THEN |
---|
[4114] | 1111 | iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content |
---|
[3999] | 1112 | ENDIF |
---|
| 1113 | ENDDO |
---|
| 1114 | |
---|
[4226] | 1115 | CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k)) |
---|
[3999] | 1116 | |
---|
| 1117 | DO i = 1, klon |
---|
| 1118 | |
---|
| 1119 | IF (rneb(i,k).GT.0.0) THEN |
---|
| 1120 | |
---|
[4072] | 1121 | ! Initialization of zrain, zsnow and zprecip: |
---|
| 1122 | zrain=0. |
---|
| 1123 | zsnow=0. |
---|
| 1124 | zprecip=0. |
---|
[4392] | 1125 | ! c_iso same init for isotopes. Externalisation? |
---|
[3999] | 1126 | |
---|
| 1127 | IF (zneb(i).EQ.seuil_neb) THEN |
---|
[4072] | 1128 | zprecip = 0.0 |
---|
| 1129 | zsnow = 0.0 |
---|
| 1130 | zrain= 0.0 |
---|
[3999] | 1131 | ELSE |
---|
[4420] | 1132 | |
---|
[3999] | 1133 | IF (ptconv(i,k)) THEN ! if convective point |
---|
| 1134 | zcl=cld_lc_con |
---|
| 1135 | zct=1./cld_tau_con |
---|
[4559] | 1136 | zexpo=cld_expo_con |
---|
| 1137 | ffallv=ffallv_con |
---|
[3999] | 1138 | ELSE |
---|
| 1139 | zcl=cld_lc_lsc |
---|
| 1140 | zct=1./cld_tau_lsc |
---|
[4559] | 1141 | zexpo=cld_expo_lsc |
---|
| 1142 | ffallv=ffallv_lsc |
---|
[3999] | 1143 | ENDIF |
---|
| 1144 | |
---|
[4559] | 1145 | |
---|
[3999] | 1146 | ! if vertical heterogeneity is taken into account, we use |
---|
| 1147 | ! the "true" volume fraction instead of a modified |
---|
| 1148 | ! surface fraction (which is larger and artificially |
---|
| 1149 | ! reduces the in-cloud water). |
---|
[4072] | 1150 | |
---|
[4559] | 1151 | ! Liquid water quantity to remove: zchau (Sundqvist, 1978) |
---|
| 1152 | ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2) |
---|
| 1153 | !......................................................... |
---|
[3999] | 1154 | IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN |
---|
| 1155 | |
---|
[4420] | 1156 | ! if vertical heterogeneity is taken into account, we use |
---|
| 1157 | ! the "true" volume fraction instead of a modified |
---|
| 1158 | ! surface fraction (which is larger and artificially |
---|
| 1159 | ! reduces the in-cloud water). |
---|
[4559] | 1160 | effective_zneb=ctot_vol(i,k) |
---|
| 1161 | ELSE |
---|
| 1162 | effective_zneb=zneb(i) |
---|
| 1163 | ENDIF |
---|
[4420] | 1164 | |
---|
| 1165 | |
---|
[4559] | 1166 | IF (iflag_autoconversion .EQ. 2) THEN |
---|
| 1167 | ! two-steps resolution with niter_lscp=1 sufficient |
---|
| 1168 | ! we first treat the second term (with exponential) in an explicit way |
---|
| 1169 | ! and then treat the first term (-q/tau) in an exact way |
---|
| 1170 | zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct & |
---|
| 1171 | *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo)))) |
---|
[4420] | 1172 | ELSE |
---|
[4559] | 1173 | ! old explicit resolution with subtimesteps |
---|
| 1174 | zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) & |
---|
| 1175 | *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo)) |
---|
[4420] | 1176 | ENDIF |
---|
| 1177 | |
---|
| 1178 | |
---|
[4559] | 1179 | ! Ice water quantity to remove (Zender & Kiehl, 1997) |
---|
| 1180 | ! dqice/dt=1/rho*d(rho*wice*qice)/dz |
---|
| 1181 | !.................................... |
---|
| 1182 | IF (iflag_autoconversion .EQ. 2) THEN |
---|
| 1183 | ! exact resolution, niter_lscp=1 is sufficient but works only |
---|
| 1184 | ! with iflag_vice=0 |
---|
| 1185 | IF (zoliqi(i) .GT. 0.) THEN |
---|
| 1186 | zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) & |
---|
| 1187 | +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo)) |
---|
| 1188 | ELSE |
---|
| 1189 | zfroi=0. |
---|
| 1190 | ENDIF |
---|
| 1191 | ELSE |
---|
| 1192 | ! old explicit resolution with subtimesteps |
---|
| 1193 | zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*velo(i,k) |
---|
| 1194 | ENDIF |
---|
| 1195 | |
---|
[4397] | 1196 | zrain = MIN(MAX(zchau,0.0),zoliql(i)) |
---|
| 1197 | zsnow = MIN(MAX(zfroi,0.0),zoliqi(i)) |
---|
[4072] | 1198 | zprecip = MAX(zrain + zsnow,0.0) |
---|
[3999] | 1199 | |
---|
| 1200 | ENDIF |
---|
[4226] | 1201 | |
---|
[4559] | 1202 | IF (iflag_autoconversion .GE. 1) THEN |
---|
| 1203 | ! debugged version with phase conservation through the autoconversion process |
---|
| 1204 | zoliql(i) = MAX(zoliql(i)-1.*zrain , 0.0) |
---|
| 1205 | zoliqi(i) = MAX(zoliqi(i)-1.*zsnow , 0.0) |
---|
| 1206 | zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) |
---|
| 1207 | ELSE |
---|
| 1208 | ! bugged version with phase resetting |
---|
| 1209 | zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain , 0.0) |
---|
| 1210 | zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow , 0.0) |
---|
| 1211 | zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) |
---|
| 1212 | ENDIF |
---|
[3999] | 1213 | |
---|
[4392] | 1214 | ! c_iso: call isotope_conversion (for readibility) or duplicate |
---|
| 1215 | |
---|
[4559] | 1216 | radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1) |
---|
| 1217 | radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1) |
---|
| 1218 | radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1) |
---|
[4226] | 1219 | |
---|
[4072] | 1220 | ENDIF ! rneb >0 |
---|
[3999] | 1221 | |
---|
| 1222 | ENDDO ! i = 1,klon |
---|
| 1223 | |
---|
| 1224 | ENDDO ! n = 1,niter |
---|
| 1225 | |
---|
[4072] | 1226 | ! Precipitation flux calculation |
---|
[3999] | 1227 | |
---|
| 1228 | DO i = 1, klon |
---|
[4380] | 1229 | |
---|
[4563] | 1230 | IF (iflag_evap_prec.GE.4) THEN |
---|
[4380] | 1231 | ziflprev(i)=ziflcld(i) |
---|
| 1232 | ELSE |
---|
| 1233 | ziflprev(i)=zifl(i)*zneb(i) |
---|
| 1234 | ENDIF |
---|
[3999] | 1235 | |
---|
| 1236 | IF (rneb(i,k) .GT. 0.0) THEN |
---|
| 1237 | |
---|
[4072] | 1238 | ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux: |
---|
| 1239 | ! If T<0C, liquid precip are converted into ice, which leads to an increase in |
---|
| 1240 | ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly |
---|
| 1241 | ! taken into account through a linearization of the equations and by approximating |
---|
| 1242 | ! the liquid precipitation process with a "threshold" process. We assume that |
---|
| 1243 | ! condensates are not modified during this operation. Liquid precipitation is |
---|
| 1244 | ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases. |
---|
| 1245 | ! Water vapor increases as well |
---|
| 1246 | ! Note that compared to fisrtilp, we always assume iflag_bergeron=2 |
---|
| 1247 | |
---|
| 1248 | zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i) |
---|
| 1249 | zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i)) |
---|
| 1250 | zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) |
---|
| 1251 | coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i) |
---|
| 1252 | ! Computation of DT if all the liquid precip freezes |
---|
| 1253 | DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1)) |
---|
| 1254 | ! T should not exceed the freezing point |
---|
| 1255 | ! that is Delta > RTT-zt(i) |
---|
| 1256 | DeltaT = max( min( RTT-zt(i), DeltaT) , 0. ) |
---|
| 1257 | zt(i) = zt(i) + DeltaT |
---|
| 1258 | ! water vaporization due to temp. increase |
---|
| 1259 | Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT |
---|
| 1260 | ! we add this vaporized water to the total vapor and we remove it from the precip |
---|
| 1261 | zq(i) = zq(i) + Deltaq |
---|
| 1262 | ! The three "max" lines herebelow protect from rounding errors |
---|
| 1263 | zcond(i) = max( zcond(i)- Deltaq, 0. ) |
---|
| 1264 | ! liquid precipitation converted to ice precip |
---|
| 1265 | Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT |
---|
| 1266 | zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. ) |
---|
| 1267 | ! iced water budget |
---|
| 1268 | zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.) |
---|
[4226] | 1269 | |
---|
[4392] | 1270 | ! c_iso : mv here condensation of isotopes + redispatchage en precip |
---|
| 1271 | |
---|
[4563] | 1272 | IF (iflag_evap_prec.GE.4) THEN |
---|
[4226] | 1273 | zrflcld(i) = zrflcld(i)+zqprecl(i) & |
---|
| 1274 | *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
---|
| 1275 | ziflcld(i) = ziflcld(i)+ zqpreci(i) & |
---|
| 1276 | *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
---|
| 1277 | znebprecipcld(i) = rneb(i,k) |
---|
| 1278 | zrfl(i) = zrflcld(i) + zrflclr(i) |
---|
| 1279 | zifl(i) = ziflcld(i) + ziflclr(i) |
---|
| 1280 | ELSE |
---|
[3999] | 1281 | zrfl(i) = zrfl(i)+ zqprecl(i) & |
---|
| 1282 | *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
---|
| 1283 | zifl(i) = zifl(i)+ zqpreci(i) & |
---|
[4226] | 1284 | *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
---|
[3999] | 1285 | ENDIF |
---|
[4392] | 1286 | ! c_iso : same for isotopes |
---|
[3999] | 1287 | |
---|
| 1288 | ENDIF ! rneb>0 |
---|
| 1289 | |
---|
| 1290 | ENDDO |
---|
| 1291 | |
---|
[4226] | 1292 | ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min |
---|
[4563] | 1293 | ! if iflag_evap_prec>=4 |
---|
| 1294 | IF (iflag_evap_prec.GE.4) THEN |
---|
[4114] | 1295 | |
---|
[4226] | 1296 | DO i=1,klon |
---|
[4114] | 1297 | |
---|
[4226] | 1298 | IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN |
---|
[3999] | 1299 | znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ & |
---|
| 1300 | (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min))) |
---|
| 1301 | ELSE |
---|
[4226] | 1302 | znebprecipclr(i)=0.0 |
---|
| 1303 | ENDIF |
---|
| 1304 | |
---|
| 1305 | IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN |
---|
[3999] | 1306 | znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ & |
---|
| 1307 | (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min))) |
---|
[4226] | 1308 | ELSE |
---|
| 1309 | znebprecipcld(i)=0.0 |
---|
| 1310 | ENDIF |
---|
[3999] | 1311 | |
---|
[4226] | 1312 | ENDDO |
---|
[3999] | 1313 | |
---|
[4226] | 1314 | ENDIF |
---|
[3999] | 1315 | |
---|
[4803] | 1316 | ENDIF ! ok_poprecip |
---|
| 1317 | |
---|
[3999] | 1318 | ! End of precipitation formation |
---|
| 1319 | ! -------------------------------- |
---|
| 1320 | |
---|
| 1321 | |
---|
[4803] | 1322 | ! Calculation of cloud condensates amount seen by radiative scheme |
---|
| 1323 | !----------------------------------------------------------------- |
---|
| 1324 | |
---|
[4114] | 1325 | ! Calculation of the concentration of condensates seen by the radiation scheme |
---|
[4412] | 1326 | ! for liquid phase, we take radocondl |
---|
| 1327 | ! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true) |
---|
[4803] | 1328 | ! we recalculate radocondi to account for contributions from the precipitation flux |
---|
| 1329 | ! TODO: for the moment, we deactivate this option when ok_poprecip=.true. |
---|
[3999] | 1330 | |
---|
[4803] | 1331 | IF ((ok_radocond_snow) .AND. (k .LT. klev) .AND. (.NOT. ok_poprecip)) THEN |
---|
[4114] | 1332 | ! for the solid phase (crystals + snowflakes) |
---|
[4412] | 1333 | ! we recalculate radocondi by summing |
---|
[4114] | 1334 | ! the ice content calculated in the mesh |
---|
| 1335 | ! + the contribution of the non-evaporated snowfall |
---|
| 1336 | ! from the overlying layer |
---|
| 1337 | DO i=1,klon |
---|
| 1338 | IF (ziflprev(i) .NE. 0.0) THEN |
---|
[4412] | 1339 | radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1) |
---|
[4114] | 1340 | ELSE |
---|
[4412] | 1341 | radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i) |
---|
[4114] | 1342 | ENDIF |
---|
[4412] | 1343 | radocond(i,k)=radocondl(i,k)+radocondi(i,k) |
---|
[4114] | 1344 | ENDDO |
---|
| 1345 | ENDIF |
---|
| 1346 | |
---|
[4412] | 1347 | ! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme |
---|
[4114] | 1348 | DO i=1,klon |
---|
[4412] | 1349 | IF (radocond(i,k) .GT. 0.) THEN |
---|
| 1350 | radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.) |
---|
[4114] | 1351 | ENDIF |
---|
| 1352 | ENDDO |
---|
| 1353 | |
---|
[3999] | 1354 | ! ---------------------------------------------------------------- |
---|
| 1355 | ! P4> Wet scavenging |
---|
| 1356 | ! ---------------------------------------------------------------- |
---|
| 1357 | |
---|
| 1358 | !Scavenging through nucleation in the layer |
---|
| 1359 | |
---|
| 1360 | DO i = 1,klon |
---|
| 1361 | |
---|
| 1362 | IF(zcond(i).GT.zoliq(i)+1.e-10) THEN |
---|
| 1363 | beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime |
---|
| 1364 | ELSE |
---|
| 1365 | beta(i,k) = 0. |
---|
| 1366 | ENDIF |
---|
| 1367 | |
---|
[4059] | 1368 | zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG |
---|
[3999] | 1369 | |
---|
| 1370 | IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN |
---|
| 1371 | |
---|
[4654] | 1372 | IF (temp(i,k) .GE. t_glace_min) THEN |
---|
[3999] | 1373 | zalpha_tr = a_tr_sca(3) |
---|
| 1374 | ELSE |
---|
| 1375 | zalpha_tr = a_tr_sca(4) |
---|
| 1376 | ENDIF |
---|
| 1377 | |
---|
| 1378 | zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) |
---|
| 1379 | frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi |
---|
[4226] | 1380 | |
---|
[3999] | 1381 | ! Nucleation with a factor of -1 instead of -0.5 |
---|
| 1382 | zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i)) |
---|
| 1383 | |
---|
| 1384 | ENDIF |
---|
| 1385 | |
---|
[4226] | 1386 | ENDDO |
---|
| 1387 | |
---|
[3999] | 1388 | ! Scavenging through impaction in the underlying layer |
---|
| 1389 | |
---|
| 1390 | DO kk = k-1, 1, -1 |
---|
| 1391 | |
---|
| 1392 | DO i = 1, klon |
---|
| 1393 | |
---|
| 1394 | IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN |
---|
| 1395 | |
---|
[4654] | 1396 | IF (temp(i,kk) .GE. t_glace_min) THEN |
---|
[3999] | 1397 | zalpha_tr = a_tr_sca(1) |
---|
| 1398 | ELSE |
---|
| 1399 | zalpha_tr = a_tr_sca(2) |
---|
| 1400 | ENDIF |
---|
| 1401 | |
---|
| 1402 | zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) |
---|
| 1403 | frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi |
---|
| 1404 | |
---|
| 1405 | ENDIF |
---|
| 1406 | |
---|
| 1407 | ENDDO |
---|
| 1408 | |
---|
| 1409 | ENDDO |
---|
[4226] | 1410 | |
---|
[4072] | 1411 | !--save some variables for ice supersaturation |
---|
[4059] | 1412 | ! |
---|
| 1413 | DO i = 1, klon |
---|
[4072] | 1414 | ! for memory |
---|
[4059] | 1415 | rneb_seri(i,k) = rneb(i,k) |
---|
[3999] | 1416 | |
---|
[4072] | 1417 | ! for diagnostics |
---|
[4059] | 1418 | rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k) |
---|
[3999] | 1419 | |
---|
[4059] | 1420 | qvc(i,k) = zqs(i) * rneb(i,k) |
---|
[4072] | 1421 | 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 |
---|
[4059] | 1422 | qcld(i,k) = qvc(i,k) + zcond(i) |
---|
[4226] | 1423 | ENDDO |
---|
[4072] | 1424 | !q_sat |
---|
[4380] | 1425 | CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs(:)) |
---|
| 1426 | CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:)) |
---|
[4059] | 1427 | |
---|
[4803] | 1428 | |
---|
| 1429 | |
---|
| 1430 | ! Outputs: |
---|
| 1431 | !------------------------------- |
---|
| 1432 | ! Precipitation fluxes at layer interfaces |
---|
| 1433 | ! + precipitation fractions + |
---|
| 1434 | ! temperature and water species tendencies |
---|
| 1435 | DO i = 1, klon |
---|
| 1436 | psfl(i,k)=zifl(i) |
---|
| 1437 | prfl(i,k)=zrfl(i) |
---|
| 1438 | pfraclr(i,k)=znebprecipclr(i) |
---|
| 1439 | pfracld(i,k)=znebprecipcld(i) |
---|
| 1440 | d_ql(i,k) = (1-zfice(i))*zoliq(i) |
---|
| 1441 | d_qi(i,k) = zfice(i)*zoliq(i) |
---|
| 1442 | d_q(i,k) = zq(i) - qt(i,k) |
---|
| 1443 | ! c_iso: same for isotopes |
---|
| 1444 | d_t(i,k) = zt(i) - temp(i,k) |
---|
| 1445 | ENDDO |
---|
| 1446 | |
---|
| 1447 | |
---|
[4226] | 1448 | ENDDO |
---|
[4059] | 1449 | |
---|
[3999] | 1450 | |
---|
| 1451 | ! Rain or snow at the surface (depending on the first layer temperature) |
---|
| 1452 | DO i = 1, klon |
---|
| 1453 | snow(i) = zifl(i) |
---|
| 1454 | rain(i) = zrfl(i) |
---|
[4392] | 1455 | ! c_iso final output |
---|
[3999] | 1456 | ENDDO |
---|
| 1457 | |
---|
| 1458 | IF (ncoreczq>0) THEN |
---|
| 1459 | WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.' |
---|
| 1460 | ENDIF |
---|
| 1461 | |
---|
[4654] | 1462 | |
---|
| 1463 | RETURN |
---|
| 1464 | |
---|
[4380] | 1465 | END SUBROUTINE lscp |
---|
[3999] | 1466 | !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 1467 | |
---|
[4664] | 1468 | END MODULE lmdz_lscp |
---|