[5614] | 1 | MODULE lmdz_lscp_main |
---|
| 2 | |
---|
| 3 | IMPLICIT NONE |
---|
| 4 | |
---|
| 5 | CONTAINS |
---|
| 6 | |
---|
| 7 | !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 8 | SUBROUTINE lscp(klon, klev, dtime, missing_val, & |
---|
| 9 | paprs, pplay, omega, temp, qt, ql_seri, qi_seri, & |
---|
| 10 | ptconv, ratqs, sigma_qtherm, & |
---|
| 11 | d_t, d_q, d_ql, d_qi, rneb, rneblsvol, & |
---|
| 12 | pfraclr, pfracld, & |
---|
| 13 | cldfraliq, cldfraliqth, & |
---|
| 14 | sigma2_icefracturb,sigma2_icefracturbth, & |
---|
| 15 | mean_icefracturb,mean_icefracturbth, & |
---|
| 16 | radocond, radicefrac, rain, snow, & |
---|
| 17 | frac_impa, frac_nucl, beta, & |
---|
| 18 | prfl, psfl, rhcl, qta, fraca, & |
---|
| 19 | tv, pspsk, tla, thl, wth, iflag_cld_th, & |
---|
| 20 | iflag_ice_thermo, distcltop, temp_cltop, & |
---|
| 21 | tke, tke_dissip, & |
---|
| 22 | entr_therm, detr_therm, & |
---|
| 23 | cell_area, & |
---|
| 24 | cf_seri, rvc_seri, u_seri, v_seri, & |
---|
| 25 | qsub, qissr, qcld, subfra, issrfra, gamma_cond, & |
---|
| 26 | dcf_sub, dcf_con, dcf_mix, & |
---|
| 27 | dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, & |
---|
| 28 | dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati, & |
---|
| 29 | Tcontr, qcontr, qcontr2, fcontrN, fcontrP, dcf_avi,& |
---|
| 30 | dqi_avi, dqvc_avi, flight_dist, flight_h2o, & |
---|
| 31 | cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, & |
---|
| 32 | qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, & |
---|
| 33 | dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,& |
---|
| 34 | dqsmelt, dqsfreez) |
---|
| 35 | |
---|
| 36 | !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 37 | ! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD) |
---|
| 38 | ! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD) |
---|
| 39 | !-------------------------------------------------------------------------------- |
---|
| 40 | ! Date: 01/2021 |
---|
| 41 | !-------------------------------------------------------------------------------- |
---|
| 42 | ! Aim: Large Scale Clouds and Precipitation (LSCP) |
---|
| 43 | ! |
---|
| 44 | ! This code is a new version of the fisrtilp.F90 routine, which is itself a |
---|
| 45 | ! merge of 'first' (superrsaturation physics, P. LeVan K. Laval) |
---|
| 46 | ! and 'ilp' (il pleut, L. Li) |
---|
| 47 | ! |
---|
| 48 | ! Compared to the original fisrtilp code, lscp |
---|
| 49 | ! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.) |
---|
| 50 | ! -> consider always precipitation thermalisation (fl_cor_ebil>0) |
---|
| 51 | ! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T) |
---|
| 52 | ! -> option "oldbug" by JYG has been removed |
---|
| 53 | ! -> iflag_t_glace >0 always |
---|
| 54 | ! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always) |
---|
| 55 | ! -> rectangular distribution from L. Li no longer available |
---|
| 56 | ! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt) |
---|
| 57 | !-------------------------------------------------------------------------------- |
---|
| 58 | ! References: |
---|
| 59 | ! |
---|
| 60 | ! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2 |
---|
| 61 | ! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y |
---|
| 62 | ! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3 |
---|
| 63 | ! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379 |
---|
| 64 | ! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046 |
---|
| 65 | ! - Touzze-Peifert Ludo, PhD thesis, p117-124 |
---|
| 66 | ! ------------------------------------------------------------------------------- |
---|
| 67 | ! Code structure: |
---|
| 68 | ! |
---|
| 69 | ! P0> Thermalization of the precipitation coming from the overlying layer |
---|
| 70 | ! P1> Evaporation of the precipitation (falling from the k+1 level) |
---|
| 71 | ! P2> Cloud formation (at the k level) |
---|
| 72 | ! P2.A.1> With the PDFs, calculation of cloud properties using the inital |
---|
| 73 | ! values of T and Q |
---|
| 74 | ! P2.A.2> Coupling between condensed water and temperature |
---|
| 75 | ! P2.A.3> Calculation of final quantities associated with cloud formation |
---|
| 76 | ! P2.B> Release of Latent heat after cloud formation |
---|
| 77 | ! P3> Autoconversion to precipitation (k-level) |
---|
| 78 | ! P4> Wet scavenging |
---|
| 79 | !------------------------------------------------------------------------------ |
---|
| 80 | ! Some preliminary comments (JBM) : |
---|
| 81 | ! |
---|
| 82 | ! The cloud water that the radiation scheme sees is not the same that the cloud |
---|
| 83 | ! water used in the physics and the dynamics |
---|
| 84 | ! |
---|
| 85 | ! During the autoconversion to precipitation (P3 step), radocond (cloud water used |
---|
| 86 | ! by the radiation scheme) is calculated as an average of the water that remains |
---|
| 87 | ! in the cloud during the precipitation and not the water remaining at the end |
---|
| 88 | ! of the time step. The latter is used in the rest of the physics and advected |
---|
| 89 | ! by the dynamics. |
---|
| 90 | ! |
---|
| 91 | ! In summary: |
---|
| 92 | ! |
---|
| 93 | ! Radiation: |
---|
| 94 | ! xflwc(newmicro)+xfiwc(newmicro) = |
---|
| 95 | ! radocond=lwcon(nc)+iwcon(nc) |
---|
| 96 | ! |
---|
| 97 | ! Notetheless, be aware of: |
---|
| 98 | ! |
---|
| 99 | ! radocond .NE. ocond(nc) |
---|
| 100 | ! i.e.: |
---|
| 101 | ! lwcon(nc)+iwcon(nc) .NE. ocond(nc) |
---|
| 102 | ! but oliq+(ocond-oliq) .EQ. ocond |
---|
| 103 | ! (which is not trivial) |
---|
| 104 | !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 105 | ! |
---|
| 106 | |
---|
| 107 | ! USE de modules contenant des fonctions. |
---|
| 108 | USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat |
---|
| 109 | USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top |
---|
| 110 | USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat |
---|
| 111 | USE lmdz_lscp_condensation, ONLY : cloudth, cloudth_v3, cloudth_v6, condensation_cloudth |
---|
| 112 | USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld |
---|
| 113 | USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld |
---|
| 114 | |
---|
| 115 | ! Use du module lmdz_lscp_ini contenant les constantes |
---|
| 116 | USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps |
---|
| 117 | USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0 |
---|
| 118 | USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca |
---|
| 119 | USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat |
---|
| 120 | USE lmdz_lscp_ini, ONLY : min_frac_th_cld, temp_nowater |
---|
| 121 | USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG |
---|
| 122 | USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp |
---|
| 123 | USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac |
---|
| 124 | USE lmdz_lscp_ini, ONLY : ok_lscp_mergecond, gamma_mixth |
---|
| 125 | |
---|
| 126 | IMPLICIT NONE |
---|
| 127 | |
---|
| 128 | !=============================================================================== |
---|
| 129 | ! VARIABLES DECLARATION |
---|
| 130 | !=============================================================================== |
---|
| 131 | |
---|
| 132 | ! INPUT VARIABLES: |
---|
| 133 | !----------------- |
---|
| 134 | |
---|
| 135 | INTEGER, INTENT(IN) :: klon,klev ! number of horizontal grid points and vertical levels |
---|
| 136 | INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics |
---|
| 137 | |
---|
| 138 | |
---|
| 139 | REAL, INTENT(IN) :: dtime ! time step [s] |
---|
| 140 | REAL, INTENT(IN) :: missing_val ! missing value for output |
---|
| 141 | |
---|
| 142 | REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa] |
---|
| 143 | REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa] |
---|
| 144 | REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K) |
---|
| 145 | REAL, DIMENSION(klon,klev), INTENT(IN) :: omega ! vertical velocity [Pa/s] |
---|
| 146 | REAL, DIMENSION(klon,klev), INTENT(IN) :: qt ! total specific humidity (in vapor phase in input) [kg/kg] |
---|
| 147 | REAL, DIMENSION(klon,klev), INTENT(IN) :: ql_seri ! liquid specific content [kg/kg] |
---|
| 148 | REAL, DIMENSION(klon,klev), INTENT(IN) :: qi_seri ! ice specific content [kg/kg] |
---|
| 149 | REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke ! turbulent kinetic energy [m2/s2] |
---|
| 150 | REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke_dissip ! TKE dissipation [m2/s3] |
---|
[5647] | 151 | REAL, DIMENSION(klon,klev), INTENT(IN) :: entr_therm ! thermal plume entrainment rate * dz [kg/s/m2] |
---|
| 152 | REAL, DIMENSION(klon,klev), INTENT(IN) :: detr_therm ! thermal plume detrainment rate * dz [kg/s/m2] |
---|
[5614] | 153 | |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active |
---|
| 157 | |
---|
| 158 | !Inputs associated with thermal plumes |
---|
| 159 | |
---|
| 160 | INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds |
---|
| 161 | REAL, DIMENSION(klon,klev), INTENT(IN) :: tv ! virtual potential temperature [K] |
---|
| 162 | REAL, DIMENSION(klon,klev), INTENT(IN) :: qta ! specific humidity within thermals [kg/kg] |
---|
| 163 | REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-] |
---|
| 164 | REAL, DIMENSION(klon,klev), INTENT(IN) :: pspsk ! exner potential (p/100000)**(R/cp) |
---|
| 165 | REAL, DIMENSION(klon,klev), INTENT(IN) :: tla ! liquid potential temperature within thermals [K] |
---|
| 166 | REAL, DIMENSION(klon,klev), INTENT(IN) :: wth ! vertical velocity within thermals [m/s] |
---|
| 167 | REAL, DIMENSION(klon,klev), INTENT(IN) :: sigma_qtherm ! controls the saturation deficit distribution width in thermals [-] |
---|
| 168 | |
---|
| 169 | |
---|
| 170 | |
---|
| 171 | ! INPUT/OUTPUT variables |
---|
| 172 | !------------------------ |
---|
| 173 | |
---|
| 174 | REAL, DIMENSION(klon,klev), INTENT(INOUT) :: thl ! liquid potential temperature [K] |
---|
| 175 | REAL, DIMENSION(klon,klev), INTENT(INOUT) :: ratqs ! function that sets the large-scale water distribution width |
---|
| 176 | |
---|
| 177 | |
---|
| 178 | ! INPUT/OUTPUT condensation and ice supersaturation |
---|
| 179 | !-------------------------------------------------- |
---|
| 180 | REAL, DIMENSION(klon,klev), INTENT(INOUT):: cf_seri ! cloud fraction [-] |
---|
| 181 | REAL, DIMENSION(klon,klev), INTENT(INOUT):: rvc_seri ! cloudy water vapor to total water vapor ratio [-] |
---|
| 182 | REAL, DIMENSION(klon,klev), INTENT(IN) :: u_seri ! eastward wind [m/s] |
---|
| 183 | REAL, DIMENSION(klon,klev), INTENT(IN) :: v_seri ! northward wind [m/s] |
---|
| 184 | REAL, DIMENSION(klon), INTENT(IN) :: cell_area ! area of each cell [m2] |
---|
| 185 | |
---|
| 186 | ! INPUT/OUTPUT aviation |
---|
| 187 | !-------------------------------------------------- |
---|
| 188 | REAL, DIMENSION(klon,klev), INTENT(IN) :: flight_dist ! Aviation distance flown within the mesh [m/s/mesh] |
---|
| 189 | REAL, DIMENSION(klon,klev), INTENT(IN) :: flight_h2o ! Aviation H2O emitted within the mesh [kg H2O/s/mesh] |
---|
| 190 | |
---|
| 191 | ! OUTPUT variables |
---|
| 192 | !----------------- |
---|
| 193 | |
---|
| 194 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K] |
---|
| 195 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg] |
---|
| 196 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg] |
---|
| 197 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg] |
---|
| 198 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-] |
---|
| 199 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneblsvol ! cloud fraction per unit volume [-] |
---|
| 200 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-] |
---|
| 201 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-] |
---|
| 202 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: cldfraliq ! liquid fraction of cloud fraction [-] |
---|
| 203 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: cldfraliqth ! liquid fraction of cloud fraction in thermals [-] |
---|
| 204 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: sigma2_icefracturb ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-] |
---|
| 205 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: mean_icefracturb ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-] |
---|
| 206 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: sigma2_icefracturbth ! Variance of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-] |
---|
| 207 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: mean_icefracturbth ! Mean of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-] |
---|
| 208 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: radocond ! condensed water used in the radiation scheme [kg/kg] |
---|
| 209 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme |
---|
| 210 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-] |
---|
| 211 | REAL, DIMENSION(klon), INTENT(OUT) :: rain ! surface large-scale rainfall [kg/s/m2] |
---|
| 212 | REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2] |
---|
| 213 | REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2] |
---|
| 214 | REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2] |
---|
| 215 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: distcltop ! distance to cloud top [m] |
---|
| 216 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: temp_cltop ! temperature of cloud top [K] |
---|
| 217 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: beta ! conversion rate of condensed water |
---|
| 218 | |
---|
| 219 | ! fraction of aerosol scavenging through impaction and nucleation (for on-line) |
---|
| 220 | |
---|
| 221 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-] |
---|
| 222 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-] |
---|
| 223 | |
---|
| 224 | ! for condensation and ice supersaturation |
---|
| 225 | |
---|
| 226 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsub !--specific total water content in sub-saturated clear sky region [kg/kg] |
---|
| 227 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qissr !--specific total water content in supersat region [kg/kg] |
---|
| 228 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld !--specific total water content in cloudy region [kg/kg] |
---|
| 229 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: subfra !--mesh fraction of subsaturated clear sky [-] |
---|
| 230 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: issrfra !--mesh fraction of ISSR [-] |
---|
| 231 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_cond !--coefficient governing the ice nucleation RHi threshold [-] |
---|
| 232 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_sub !--cloud fraction tendency because of sublimation [s-1] |
---|
| 233 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_con !--cloud fraction tendency because of condensation [s-1] |
---|
| 234 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_mix !--cloud fraction tendency because of cloud mixing [s-1] |
---|
| 235 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_adj !--specific ice content tendency because of temperature adjustment [kg/kg/s] |
---|
| 236 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_sub !--specific ice content tendency because of sublimation [kg/kg/s] |
---|
| 237 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_con !--specific ice content tendency because of condensation [kg/kg/s] |
---|
| 238 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_mix !--specific ice content tendency because of cloud mixing [kg/kg/s] |
---|
| 239 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_adj !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s] |
---|
| 240 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_sub !--specific cloud water vapor tendency because of sublimation [kg/kg/s] |
---|
| 241 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_con !--specific cloud water vapor tendency because of condensation [kg/kg/s] |
---|
| 242 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_mix !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s] |
---|
| 243 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl !--saturation specific humidity wrt liquid [kg/kg] |
---|
| 244 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsati !--saturation specific humidity wrt ice [kg/kg] |
---|
| 245 | |
---|
| 246 | ! for contrails and aviation |
---|
| 247 | |
---|
| 248 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcontr !--threshold temperature for contrail formation [K] |
---|
| 249 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr !--threshold humidity for contrail formation [kg/kg] |
---|
| 250 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr2 !--// (2nd expression more consistent with LMDZ expression of q) |
---|
| 251 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrN !--fraction of grid favourable to non-persistent contrails |
---|
| 252 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrP !--fraction of grid favourable to persistent contrails |
---|
| 253 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_avi !--cloud fraction tendency because of aviation [s-1] |
---|
| 254 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_avi !--specific ice content tendency because of aviation [kg/kg/s] |
---|
| 255 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_avi !--specific cloud water vapor tendency because of aviation [kg/kg/s] |
---|
| 256 | |
---|
| 257 | |
---|
| 258 | ! for POPRECIP |
---|
| 259 | |
---|
| 260 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg] |
---|
| 261 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg] |
---|
| 262 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqreva !--rain tendendy due to evaporation [kg/kg/s] |
---|
| 263 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqssub !--snow tendency due to sublimation [kg/kg/s] |
---|
| 264 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s] |
---|
| 265 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s] |
---|
| 266 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s] |
---|
| 267 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s] |
---|
| 268 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsrim !--snow tendency due to riming [kg/kg/s] |
---|
| 269 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsmelt !--snow tendency due to melting [kg/kg/s] |
---|
| 270 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrmelt !--rain tendency due to melting [kg/kg/s] |
---|
| 271 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsfreez !--snow tendency due to freezing [kg/kg/s] |
---|
| 272 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrfreez !--rain tendency due to freezing [kg/kg/s] |
---|
| 273 | |
---|
| 274 | ! for thermals |
---|
| 275 | |
---|
| 276 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth !--mean saturation deficit in thermals |
---|
| 277 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv !--mean saturation deficit in environment |
---|
| 278 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath !--std of saturation deficit in thermals |
---|
| 279 | REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv !--std of saturation deficit in environment |
---|
| 280 | |
---|
| 281 | |
---|
| 282 | ! LOCAL VARIABLES: |
---|
| 283 | !---------------- |
---|
| 284 | REAL, DIMENSION(klon,klev) :: ctot, rnebth, ctot_vol |
---|
| 285 | REAL, DIMENSION(klon,klev) :: wls !-- large scalce vertical velocity [m/s] |
---|
| 286 | REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi |
---|
| 287 | REAL, DIMENSION(klon) :: zqsth, zqslth, zqsith, zdqsth, zdqslth, zdqsith |
---|
| 288 | REAL :: zdelta |
---|
| 289 | REAL, DIMENSION(klon) :: zdqsdT_raw |
---|
| 290 | REAL, DIMENSION(klon) :: gammasat,dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T |
---|
[5633] | 291 | REAL, DIMENSION(klon) :: Tbef,Tbefth,Tbefthm1,qlibef,DT ! temperature, humidity and temp. variation during condensation iteration |
---|
[5614] | 292 | REAL :: num,denom |
---|
| 293 | REAL :: cste |
---|
| 294 | REAL, DIMENSION(klon) :: qincloud |
---|
| 295 | REAL, DIMENSION(klon) :: zrfl, zifl |
---|
| 296 | REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqth, zqn, zqnl |
---|
| 297 | REAL, DIMENSION(klon) :: zoliql, zoliqi |
---|
| 298 | REAL, DIMENSION(klon) :: zt, zp |
---|
| 299 | REAL, DIMENSION(klon) :: zfice, zficeth, zficeenv, zneb, zcf, zqi_ini, zsnow |
---|
| 300 | REAL, DIMENSION(klon) :: dzfice, dzficeth, dzficeenv |
---|
| 301 | REAL, DIMENSION(klon) :: qtot, zeroklon |
---|
| 302 | ! Variables precipitation energy conservation |
---|
| 303 | REAL, DIMENSION(klon) :: zmqc |
---|
| 304 | REAL :: zalpha_tr |
---|
| 305 | REAL :: zfrac_lessi |
---|
| 306 | REAL, DIMENSION(klon) :: zprec_cond |
---|
| 307 | REAL, DIMENSION(klon) :: zlh_solid |
---|
| 308 | REAL, DIMENSION(klon) :: ztupnew |
---|
| 309 | REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl |
---|
| 310 | REAL, DIMENSION(klon) :: zrflclr, zrflcld |
---|
| 311 | REAL, DIMENSION(klon) :: ziflclr, ziflcld |
---|
| 312 | REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld |
---|
| 313 | REAL, DIMENSION(klon) :: tot_zneb |
---|
[5647] | 314 | REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop, zdeltaz |
---|
[5614] | 315 | REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl, zqliqth, zqiceth, zqvapclth, sursat_e, invtau_e ! for icefrac_lscp_turb |
---|
| 316 | |
---|
| 317 | ! for quantity of condensates seen by radiation |
---|
| 318 | REAL, DIMENSION(klon) :: zradocond, zradoice |
---|
| 319 | REAL, DIMENSION(klon) :: zrho_up, zvelo_up |
---|
| 320 | |
---|
| 321 | ! for condensation and ice supersaturation |
---|
| 322 | REAL, DIMENSION(klon) :: qvc, qvcl, shear |
---|
| 323 | REAL :: delta_z |
---|
| 324 | !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails) |
---|
| 325 | ! Constants used for calculating ratios that are advected (using a parent-child |
---|
| 326 | ! formalism). This is not done in the dynamical core because at this moment, |
---|
| 327 | ! only isotopes can use this parent-child formalism. Note that the two constants |
---|
| 328 | ! are the same as the one use in the dynamical core, being also defined in |
---|
| 329 | ! dyn3d_common/infotrac.F90 |
---|
| 330 | REAL :: min_qParent, min_ratio |
---|
| 331 | |
---|
| 332 | INTEGER i, k, kk, iter |
---|
| 333 | INTEGER, DIMENSION(klon) :: n_i |
---|
| 334 | INTEGER ncoreczq |
---|
| 335 | LOGICAL iftop |
---|
| 336 | |
---|
| 337 | LOGICAL, DIMENSION(klon) :: stratiform_or_all_distrib,pticefracturb |
---|
| 338 | LOGICAL, DIMENSION(klon) :: keepgoing |
---|
| 339 | |
---|
| 340 | CHARACTER (len = 20) :: modname = 'lscp' |
---|
| 341 | CHARACTER (len = 80) :: abort_message |
---|
| 342 | |
---|
| 343 | |
---|
| 344 | !=============================================================================== |
---|
| 345 | ! INITIALISATION |
---|
| 346 | !=============================================================================== |
---|
| 347 | |
---|
| 348 | ! Few initial checks |
---|
| 349 | |
---|
| 350 | |
---|
| 351 | IF (iflag_fisrtilp_qsat .LT. 0) THEN |
---|
| 352 | abort_message = 'lscp cannot be used with iflag_fisrtilp<0' |
---|
| 353 | CALL abort_physic(modname,abort_message,1) |
---|
| 354 | ENDIF |
---|
| 355 | |
---|
| 356 | ! AA for 'safety' reasons |
---|
| 357 | zalpha_tr = 0. |
---|
| 358 | zfrac_lessi = 0. |
---|
| 359 | beta(:,:)= 0. |
---|
| 360 | |
---|
| 361 | ! Initialisation of variables: |
---|
| 362 | |
---|
| 363 | prfl(:,:) = 0.0 |
---|
| 364 | psfl(:,:) = 0.0 |
---|
| 365 | d_t(:,:) = 0.0 |
---|
| 366 | d_q(:,:) = 0.0 |
---|
| 367 | d_ql(:,:) = 0.0 |
---|
| 368 | d_qi(:,:) = 0.0 |
---|
| 369 | rneb(:,:) = 0.0 |
---|
| 370 | rnebth(:,:)=0.0 |
---|
| 371 | pfraclr(:,:)=0.0 |
---|
| 372 | pfracld(:,:)=0.0 |
---|
| 373 | cldfraliq(:,:)=0. |
---|
| 374 | sigma2_icefracturb(:,:)=0. |
---|
| 375 | mean_icefracturb(:,:)=0. |
---|
| 376 | cldfraliqth(:,:)=0. |
---|
| 377 | sigma2_icefracturbth(:,:)=0. |
---|
| 378 | mean_icefracturbth(:,:)=0. |
---|
| 379 | radocond(:,:) = 0.0 |
---|
| 380 | radicefrac(:,:) = 0.0 |
---|
| 381 | frac_nucl(:,:) = 1.0 |
---|
| 382 | frac_impa(:,:) = 1.0 |
---|
| 383 | rain(:) = 0.0 |
---|
| 384 | snow(:) = 0.0 |
---|
| 385 | zrfl(:) = 0.0 |
---|
| 386 | zifl(:) = 0.0 |
---|
| 387 | zneb(:) = seuil_neb |
---|
| 388 | zrflclr(:) = 0.0 |
---|
| 389 | ziflclr(:) = 0.0 |
---|
| 390 | zrflcld(:) = 0.0 |
---|
| 391 | ziflcld(:) = 0.0 |
---|
| 392 | tot_zneb(:) = 0.0 |
---|
| 393 | zeroklon(:) = 0.0 |
---|
| 394 | zdistcltop(:)=0.0 |
---|
| 395 | ztemp_cltop(:) = 0.0 |
---|
| 396 | ztupnew(:)=0.0 |
---|
| 397 | ctot_vol(:,:)=0.0 |
---|
| 398 | rneblsvol(:,:)=0.0 |
---|
| 399 | znebprecip(:)=0.0 |
---|
| 400 | znebprecipclr(:)=0.0 |
---|
| 401 | znebprecipcld(:)=0.0 |
---|
| 402 | distcltop(:,:)=0. |
---|
| 403 | temp_cltop(:,:)=0. |
---|
| 404 | |
---|
| 405 | |
---|
| 406 | !--Ice supersaturation |
---|
| 407 | gamma_cond(:,:) = 1. |
---|
| 408 | qissr(:,:) = 0. |
---|
| 409 | issrfra(:,:) = 0. |
---|
| 410 | dcf_sub(:,:) = 0. |
---|
| 411 | dcf_con(:,:) = 0. |
---|
| 412 | dcf_mix(:,:) = 0. |
---|
| 413 | dqi_adj(:,:) = 0. |
---|
| 414 | dqi_sub(:,:) = 0. |
---|
| 415 | dqi_con(:,:) = 0. |
---|
| 416 | dqi_mix(:,:) = 0. |
---|
| 417 | dqvc_adj(:,:) = 0. |
---|
| 418 | dqvc_sub(:,:) = 0. |
---|
| 419 | dqvc_con(:,:) = 0. |
---|
| 420 | dqvc_mix(:,:) = 0. |
---|
| 421 | fcontrN(:,:) = 0. |
---|
| 422 | fcontrP(:,:) = 0. |
---|
| 423 | Tcontr(:,:) = missing_val |
---|
| 424 | qcontr(:,:) = missing_val |
---|
| 425 | qcontr2(:,:) = missing_val |
---|
| 426 | dcf_avi(:,:) = 0. |
---|
| 427 | dqi_avi(:,:) = 0. |
---|
| 428 | dqvc_avi(:,:) = 0. |
---|
| 429 | qvc(:) = 0. |
---|
| 430 | shear(:) = 0. |
---|
| 431 | min_qParent = 1.e-30 |
---|
| 432 | min_ratio = 1.e-16 |
---|
| 433 | |
---|
| 434 | !-- poprecip |
---|
| 435 | qraindiag(:,:)= 0. |
---|
| 436 | qsnowdiag(:,:)= 0. |
---|
| 437 | dqreva(:,:) = 0. |
---|
| 438 | dqrauto(:,:) = 0. |
---|
| 439 | dqrmelt(:,:) = 0. |
---|
| 440 | dqrfreez(:,:) = 0. |
---|
| 441 | dqrcol(:,:) = 0. |
---|
| 442 | dqssub(:,:) = 0. |
---|
| 443 | dqsauto(:,:) = 0. |
---|
| 444 | dqsrim(:,:) = 0. |
---|
| 445 | dqsagg(:,:) = 0. |
---|
| 446 | dqsfreez(:,:) = 0. |
---|
| 447 | dqsmelt(:,:) = 0. |
---|
| 448 | zqupnew(:) = 0. |
---|
| 449 | zqvapclr(:) = 0. |
---|
| 450 | |
---|
| 451 | !-- cloud phase useful variables |
---|
| 452 | wls(:,:)=-omega(:,:) / RG / (pplay(:,:)/RD/temp(:,:)) |
---|
| 453 | zqliq(:)=0. |
---|
| 454 | zqice(:)=0. |
---|
| 455 | zqvapcl(:)=0. |
---|
| 456 | zqliqth(:)=0. |
---|
| 457 | zqiceth(:)=0. |
---|
| 458 | zqvapclth(:)=0. |
---|
| 459 | sursat_e(:)=0. |
---|
| 460 | invtau_e(:)=0. |
---|
| 461 | pticefracturb(:)=.FALSE. |
---|
| 462 | |
---|
| 463 | |
---|
| 464 | !=============================================================================== |
---|
| 465 | ! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM |
---|
| 466 | !=============================================================================== |
---|
| 467 | |
---|
| 468 | ncoreczq=0 |
---|
| 469 | |
---|
| 470 | DO k = klev, 1, -1 |
---|
| 471 | |
---|
| 472 | IF (k.LE.klev-1) THEN |
---|
| 473 | iftop=.false. |
---|
| 474 | ELSE |
---|
| 475 | iftop=.true. |
---|
| 476 | ENDIF |
---|
| 477 | |
---|
| 478 | ! Initialisation temperature and specific humidity |
---|
| 479 | ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt |
---|
| 480 | ! at the end of the klon loop, a temperature incremtent d_t due to all processes |
---|
| 481 | ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated |
---|
| 482 | ! d_t = temperature tendency due to lscp |
---|
| 483 | ! The temperature of the overlying layer is updated here because needed for thermalization |
---|
| 484 | DO i = 1, klon |
---|
| 485 | zt(i)=temp(i,k) |
---|
| 486 | zq(i)=qt(i,k) |
---|
| 487 | zp(i)=pplay(i,k) |
---|
| 488 | zqi_ini(i)=qi_seri(i,k) |
---|
| 489 | zcf(i) = 0. |
---|
| 490 | zfice(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation |
---|
| 491 | dzfice(i) = 0.0 |
---|
| 492 | zficeth(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation |
---|
| 493 | dzficeth(i) = 0.0 |
---|
| 494 | zficeenv(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation |
---|
| 495 | dzficeenv(i) = 0.0 |
---|
| 496 | |
---|
| 497 | |
---|
| 498 | IF (.not. iftop) THEN |
---|
| 499 | ztupnew(i) = temp(i,k+1) + d_t(i,k+1) |
---|
| 500 | zqupnew(i) = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1) |
---|
| 501 | !--zqs(i) is the saturation specific humidity in the layer above |
---|
| 502 | zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i)) |
---|
| 503 | ENDIF |
---|
| 504 | !c_iso init of iso |
---|
| 505 | ENDDO |
---|
| 506 | |
---|
| 507 | ! -------------------------------------------------------------------- |
---|
| 508 | ! P1> Precipitation processes, before cloud formation: |
---|
| 509 | ! Thermalization of precipitation falling from the overlying layer AND |
---|
| 510 | ! Precipitation evaporation/sublimation/melting |
---|
| 511 | !--------------------------------------------------------------------- |
---|
| 512 | |
---|
| 513 | !================================================================ |
---|
| 514 | ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R) |
---|
| 515 | IF ( ok_poprecip ) THEN |
---|
| 516 | |
---|
| 517 | CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, & |
---|
| 518 | zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, & |
---|
| 519 | zqvapclr, zqupnew, & |
---|
| 520 | cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), & |
---|
| 521 | zrfl, zrflclr, zrflcld, & |
---|
| 522 | zifl, ziflclr, ziflcld, & |
---|
| 523 | dqreva(:,k), dqssub(:,k) & |
---|
| 524 | ) |
---|
| 525 | |
---|
| 526 | !================================================================ |
---|
| 527 | ELSE |
---|
| 528 | |
---|
| 529 | CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, & |
---|
| 530 | zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, & |
---|
| 531 | zrfl, zrflclr, zrflcld, & |
---|
| 532 | zifl, ziflclr, ziflcld, & |
---|
| 533 | dqreva(:,k), dqssub(:,k) & |
---|
| 534 | ) |
---|
| 535 | |
---|
| 536 | ENDIF ! (ok_poprecip) |
---|
| 537 | |
---|
[5633] | 538 | ! Calculation of qsat,L/cp*dqsat/dT and ncoreczq counter |
---|
[5614] | 539 | !------------------------------------------------------- |
---|
| 540 | |
---|
| 541 | qtot(:)=zq(:)+zmqc(:) |
---|
| 542 | CALL calc_qsat_ecmwf(klon,zt,qtot,zp,RTT,0,.false.,zqs,zdqs) |
---|
| 543 | DO i = 1, klon |
---|
| 544 | zdelta = MAX(0.,SIGN(1.,RTT-zt(i))) |
---|
| 545 | zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta) |
---|
| 546 | IF (zq(i) .LT. 1.e-15) THEN |
---|
| 547 | ncoreczq=ncoreczq+1 |
---|
| 548 | zq(i)=1.e-15 |
---|
| 549 | ENDIF |
---|
| 550 | ! c_iso: do something similar for isotopes |
---|
| 551 | |
---|
| 552 | ENDDO |
---|
| 553 | |
---|
| 554 | ! ------------------------------------------------------------------------- |
---|
| 555 | ! P2> Cloud formation including condensation and cloud phase determination |
---|
| 556 | !-------------------------------------------------------------------------- |
---|
| 557 | ! |
---|
| 558 | ! We always assume a 'fractional cloud' approach |
---|
| 559 | ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution |
---|
| 560 | ! is prescribed and depends on large scale variables and boundary layer |
---|
| 561 | ! properties) |
---|
| 562 | ! The decrease in condensed part due tu latent heating is taken into |
---|
| 563 | ! account |
---|
| 564 | ! ------------------------------------------------------------------- |
---|
| 565 | |
---|
| 566 | ! P2.1> With the PDFs (log-normal, bigaussian) |
---|
| 567 | ! cloud properties calculation with the initial values of t and q |
---|
| 568 | ! ---------------------------------------------------------------- |
---|
| 569 | |
---|
| 570 | ! initialise gammasat and stratiform_or_all_distrib |
---|
| 571 | stratiform_or_all_distrib(:)=.TRUE. |
---|
| 572 | gammasat(:)=1. |
---|
| 573 | |
---|
| 574 | ! part of the code that is supposed to become obsolete soon |
---|
| 575 | !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 576 | IF (.NOT. ok_lscp_mergecond) THEN |
---|
| 577 | IF (iflag_cld_th.GE.5) THEN |
---|
| 578 | ! Cloud cover and content in meshes affected by shallow convection, |
---|
| 579 | ! are retrieved from a bi-gaussian distribution of the saturation deficit |
---|
| 580 | ! following Jam et al. 2013 |
---|
| 581 | |
---|
| 582 | IF (iflag_cloudth_vert.LE.2) THEN |
---|
| 583 | ! Old version of Arnaud Jam |
---|
| 584 | |
---|
| 585 | CALL cloudth(klon,klev,k,tv, & |
---|
| 586 | zq,qta,fraca, & |
---|
| 587 | qincloud,ctot,pspsk,paprs,pplay,tla,thl, & |
---|
| 588 | ratqs,zqs,temp, & |
---|
| 589 | cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) |
---|
| 590 | |
---|
| 591 | |
---|
| 592 | ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN |
---|
| 593 | ! Default version of Arnaud Jam |
---|
| 594 | |
---|
| 595 | CALL cloudth_v3(klon,klev,k,tv, & |
---|
| 596 | zq,qta,fraca, & |
---|
| 597 | qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, & |
---|
| 598 | ratqs,sigma_qtherm,zqs,temp, & |
---|
| 599 | cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) |
---|
| 600 | |
---|
| 601 | |
---|
| 602 | ELSEIF (iflag_cloudth_vert.EQ.6) THEN |
---|
| 603 | ! Jean Jouhaud's version, with specific separation between surface and volume |
---|
| 604 | ! cloud fraction Decembre 2018 |
---|
| 605 | |
---|
| 606 | CALL cloudth_v6(klon,klev,k,tv, & |
---|
| 607 | zq,qta,fraca, & |
---|
| 608 | qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, & |
---|
| 609 | ratqs,zqs,temp, & |
---|
| 610 | cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) |
---|
| 611 | |
---|
| 612 | ENDIF |
---|
| 613 | |
---|
| 614 | |
---|
| 615 | DO i=1,klon |
---|
| 616 | rneb(i,k)=ctot(i,k) |
---|
| 617 | ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.) |
---|
| 618 | rneblsvol(i,k)=ctot_vol(i,k) |
---|
| 619 | zqn(i)=qincloud(i) |
---|
| 620 | !--AB grid-mean vapor in the cloud - we assume saturation adjustment |
---|
| 621 | qvc(i) = rneb(i,k) * zqs(i) |
---|
| 622 | ENDDO |
---|
| 623 | |
---|
| 624 | ! Cloud phase final determination for clouds after "old" cloudth calls |
---|
| 625 | ! for those clouds, only temperature dependent phase partitioning (eventually modulated by |
---|
| 626 | ! distance to cloud top) is available |
---|
| 627 | ! compute distance to cloud top when cloud phase is computed depending on temperature |
---|
| 628 | ! and distance to cloud top |
---|
| 629 | IF (iflag_t_glace .GE. 4) THEN |
---|
| 630 | CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop) |
---|
| 631 | ENDIF |
---|
| 632 | CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice) |
---|
| 633 | |
---|
| 634 | ENDIF |
---|
| 635 | |
---|
| 636 | IF (iflag_cld_th .EQ. 5) THEN |
---|
| 637 | ! When iflag_cld_th=5, we always assume |
---|
| 638 | ! bi-gaussian distribution |
---|
| 639 | stratiform_or_all_distrib(:) = .FALSE. |
---|
| 640 | |
---|
| 641 | ELSEIF (iflag_cld_th .GE. 6) THEN |
---|
| 642 | ! stratiform distribution only when no thermals |
---|
| 643 | stratiform_or_all_distrib(:) = fraca(:,k) < min_frac_th_cld |
---|
| 644 | |
---|
| 645 | ENDIF |
---|
| 646 | |
---|
| 647 | ENDIF ! .not. ok_lscp_mergecond |
---|
| 648 | !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 649 | |
---|
| 650 | DT(:) = 0. |
---|
| 651 | n_i(:)=0 |
---|
| 652 | Tbef(:)=zt(:) |
---|
| 653 | qlibef(:)=0. |
---|
| 654 | Tbefth(:)=tla(:,k)*pspsk(:,k) |
---|
[5633] | 655 | IF (k .GT. 1) THEN |
---|
| 656 | Tbefthm1(:)=tla(:,k-1)*pspsk(:,k-1) |
---|
| 657 | ELSE |
---|
| 658 | Tbefthm1(:)=Tbefth(:) |
---|
| 659 | ENDIF |
---|
[5647] | 660 | zqth(:)=qta(:,k) |
---|
| 661 | zdeltaz(:)=(paprs(:,k)-paprs(:,k+1))/RG/zp(:)*RD*zt(:) |
---|
[5614] | 662 | |
---|
| 663 | ! Treatment of stratiform clouds (lognormale or ice-sursat) or all clouds (including cloudth |
---|
| 664 | ! in case of ok_lscp_mergecond) |
---|
| 665 | ! We iterate here to take into account the change in qsat(T) and ice fraction |
---|
| 666 | ! during the condensation process |
---|
| 667 | ! the increment in temperature is calculated using the first principle of |
---|
| 668 | ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction |
---|
| 669 | ! and a clear sky fraction) |
---|
| 670 | ! note that we assume that the vapor in the cloud is at saturation for this calculation |
---|
| 671 | |
---|
| 672 | DO iter=1,iflag_fisrtilp_qsat+1 |
---|
| 673 | |
---|
| 674 | keepgoing(:) = .FALSE. |
---|
| 675 | |
---|
| 676 | DO i=1,klon |
---|
| 677 | |
---|
| 678 | ! keepgoing = .true. while convergence is not satisfied |
---|
| 679 | |
---|
| 680 | IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. stratiform_or_all_distrib(i)) THEN |
---|
| 681 | |
---|
| 682 | ! if not convergence: |
---|
| 683 | ! we calculate a new iteration |
---|
| 684 | keepgoing(i) = .TRUE. |
---|
| 685 | |
---|
| 686 | ! P2.2.1> cloud fraction and condensed water mass calculation |
---|
| 687 | ! Calculated variables: |
---|
| 688 | ! rneb : cloud fraction |
---|
| 689 | ! zqn : total water within the cloud |
---|
| 690 | ! zcond: mean condensed water within the mesh |
---|
| 691 | ! rhcl: clear-sky relative humidity |
---|
| 692 | !--------------------------------------------------------------- |
---|
| 693 | |
---|
| 694 | ! new temperature that only serves in the iteration process: |
---|
| 695 | Tbef(i)=Tbef(i)+DT(i) |
---|
| 696 | |
---|
| 697 | ! total water including mass of precip |
---|
| 698 | qtot(i)=zq(i)+zmqc(i) |
---|
| 699 | |
---|
| 700 | ENDIF |
---|
| 701 | |
---|
| 702 | ENDDO |
---|
| 703 | |
---|
| 704 | ! Calculation of saturation specific humidity |
---|
| 705 | CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,0,.false.,zqs,zdqs) |
---|
| 706 | ! also in thermals |
---|
| 707 | CALL calc_qsat_ecmwf(klon,Tbefth,zqth,zp,RTT,0,.false.,zqsth,zdqsth) |
---|
| 708 | |
---|
| 709 | IF (iflag_icefrac .GE. 1) THEN |
---|
| 710 | ! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction |
---|
| 711 | ! and cloud condensed water content. Principle from Dietlitcher et al. 2018, GMD |
---|
| 712 | ! we make this option works only for the physically-based and tke-dependent param from Raillard et al. |
---|
| 713 | ! (iflag_icefrac>=1) |
---|
| 714 | CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,1,.false.,zqsl,zdqsl) |
---|
| 715 | CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,2,.false.,zqsi,zdqsi) |
---|
| 716 | DO i=1,klon |
---|
| 717 | zqs(i)=zfice(i)*zqsi(i)+(1.-zfice(i))*zqsl(i) |
---|
| 718 | zdqs(i)=zfice(i)*zdqsi(i)+zqsi(i)*dzfice(i)+(1.-zfice(i))*zdqsl(i)-zqsl(i)*dzfice(i) |
---|
| 719 | ENDDO |
---|
| 720 | ENDIF |
---|
| 721 | IF (iflag_icefrac .GE. 2) THEN |
---|
| 722 | ! same adjustment for thermals |
---|
| 723 | CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,1,.false.,zqslth,zdqslth) |
---|
| 724 | CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,2,.false.,zqsith,zdqsith) |
---|
| 725 | DO i=1,klon |
---|
| 726 | zqsth(i)=zficeth(i)*zqsith(i)+(1.-zficeth(i))*zqslth(i) |
---|
| 727 | zdqsth(i)=zficeth(i)*zdqsith(i)+zqsith(i)*dzficeth(i)+(1.-zficeth(i))*zdqslth(i)-zqslth(i)*dzficeth(i) |
---|
| 728 | ENDDO |
---|
| 729 | ENDIF |
---|
| 730 | |
---|
| 731 | CALL calc_gammasat(klon,Tbef,qtot,zp,gammasat,dgammasatdt) |
---|
| 732 | ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs |
---|
| 733 | zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:) |
---|
| 734 | |
---|
| 735 | ! Cloud condensation based on subgrid pdf |
---|
| 736 | !--AB Activates a condensation scheme that allows for |
---|
| 737 | !--ice supersaturation and contrails evolution from aviation |
---|
| 738 | IF (ok_ice_supersat) THEN |
---|
| 739 | |
---|
| 740 | !--Calculate the shear value (input for condensation and ice supersat) |
---|
| 741 | DO i = 1, klon |
---|
| 742 | !--Cell thickness [m] |
---|
| 743 | delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD |
---|
| 744 | IF ( iftop ) THEN |
---|
| 745 | shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. & |
---|
| 746 | + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. ) |
---|
| 747 | ELSEIF ( k .EQ. 1 ) THEN |
---|
| 748 | ! surface |
---|
| 749 | shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. & |
---|
| 750 | + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. ) |
---|
| 751 | ELSE |
---|
| 752 | ! other layers |
---|
| 753 | shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. & |
---|
| 754 | - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. & |
---|
| 755 | + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. & |
---|
| 756 | - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. ) |
---|
| 757 | ENDIF |
---|
| 758 | ENDDO |
---|
| 759 | |
---|
| 760 | !--------------------------------------------- |
---|
| 761 | !-- CONDENSATION AND ICE SUPERSATURATION -- |
---|
| 762 | !--------------------------------------------- |
---|
| 763 | |
---|
| 764 | CALL condensation_ice_supersat( & |
---|
| 765 | klon, dtime, missing_val, & |
---|
| 766 | zp, paprs(:,k), paprs(:,k+1), & |
---|
| 767 | cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), & |
---|
| 768 | shear, tke_dissip(:,k), cell_area, & |
---|
| 769 | Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, & |
---|
| 770 | rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), & |
---|
| 771 | dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), & |
---|
| 772 | dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), & |
---|
| 773 | dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), & |
---|
| 774 | Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), & |
---|
| 775 | flight_dist(:,k), flight_h2o(:,k), & |
---|
| 776 | dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k)) |
---|
| 777 | |
---|
| 778 | |
---|
| 779 | ELSE |
---|
| 780 | !--generalised lognormal condensation scheme (Bony and Emanuel 2001) |
---|
| 781 | |
---|
| 782 | CALL condensation_lognormal( & |
---|
| 783 | klon, Tbef, zq, zqs, gammasat, ratqs(:,k), & |
---|
| 784 | keepgoing, rneb(:,k), zqn, qvc) |
---|
| 785 | |
---|
| 786 | |
---|
| 787 | ENDIF ! .NOT. ok_ice_supersat |
---|
| 788 | |
---|
| 789 | ! Volume cloud fraction |
---|
| 790 | rneblsvol(:,k)=rneb(:,k) |
---|
| 791 | |
---|
| 792 | |
---|
| 793 | IF (ok_lscp_mergecond) THEN |
---|
| 794 | ! in that case we overwrite rneb, rneblsvol and zqn for shallow convective clouds only |
---|
| 795 | CALL condensation_cloudth(klon,Tbef,zq,zqth,fraca(:,k), & |
---|
| 796 | pspsk(:,k),zp,tla(:,k), & |
---|
| 797 | ratqs(:,k),sigma_qtherm(:,k),zqsth,zqs,zqn,rneb(:,k),rnebth(:,k),rneblsvol(:,k), & |
---|
| 798 | cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k)) |
---|
| 799 | ENDIF |
---|
| 800 | |
---|
| 801 | |
---|
| 802 | |
---|
| 803 | ! Cloud phase determination |
---|
| 804 | |
---|
| 805 | |
---|
| 806 | IF (iflag_icefrac .LE. 1) THEN |
---|
| 807 | ! "old" phase partitioning depending on temperature and eventually distance to cloud top everywhere |
---|
| 808 | IF (iflag_t_glace .GE. 4) THEN |
---|
| 809 | CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop) |
---|
| 810 | ENDIF |
---|
| 811 | CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice) |
---|
| 812 | ENDIF |
---|
| 813 | |
---|
| 814 | IF (iflag_icefrac .EQ. 1) THEN |
---|
| 815 | ! phase partitioning depending on turbulence, vertical velocity and ice crystal microphysics |
---|
| 816 | ! only in stratiform clouds in the mixed phase regime (Raillard et al. 2025) |
---|
| 817 | ! it overwrites temperature-dependent phase partitioning except for convective boundary layer clouds |
---|
| 818 | pticefracturb(:) = (fraca(:,k) < min_frac_th_cld) .AND. (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT) |
---|
| 819 | DO i=1,klon |
---|
| 820 | qincloud(i)=zqn(i) |
---|
| 821 | zcf(i)=MIN(MAX(rneb(i,k), 0.),1.) |
---|
| 822 | sursat_e(i) = 0. |
---|
| 823 | invtau_e(i) = 0. |
---|
| 824 | ENDDO |
---|
[5633] | 825 | CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), zqi_ini, & |
---|
| 826 | ziflcld, znebprecipcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, & |
---|
| 827 | zficeenv, dzficeenv, cldfraliq(:,k),sigma2_icefracturb(:,k),mean_icefracturb(:,k)) |
---|
[5614] | 828 | DO i=1,klon |
---|
| 829 | IF (pticefracturb(i)) THEN |
---|
| 830 | zfice(i)=zficeenv(i) |
---|
| 831 | dzfice(i)=dzficeenv(i) |
---|
| 832 | ENDIF |
---|
| 833 | ENDDO |
---|
| 834 | |
---|
| 835 | |
---|
| 836 | ELSEIF (iflag_icefrac .GE. 2) THEN |
---|
| 837 | ! compute also phase partitioning also in thermal clouds (neglecting tke in thermals as first assumption) |
---|
| 838 | ! moreover, given the upward velocity of thermals, we assume that precipitation falls in the environment |
---|
| 839 | ! hence we assume that no snow falls in thermals. |
---|
| 840 | ! we activate only in the mixed phase regime not to distrub the cirrus param at very cold T |
---|
| 841 | pticefracturb(:) = (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT) |
---|
| 842 | !Thermals |
---|
| 843 | DO i=1,klon |
---|
| 844 | IF (fraca(i,k) .GT. min_frac_th_cld) THEN |
---|
| 845 | zcf(i)=MIN(MAX(rnebth(i,k),0.), 1.)*fraca(i,k) |
---|
| 846 | qincloud(i)=zqn(i)/fraca(i,k) |
---|
| 847 | ELSE |
---|
| 848 | zcf(i) = 0. |
---|
| 849 | qincloud(i) = 0. |
---|
| 850 | ENDIF |
---|
| 851 | sursat_e(i)=cloudth_senv(i,k)/zqsi(i) |
---|
[5647] | 852 | invtau_e(i)=gamma_mixth*MAX(entr_therm(i,k)-detr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i) |
---|
[5614] | 853 | ENDDO |
---|
[5633] | 854 | CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbefth, zp, paprs(:,k), paprs(:,k+1), wth(:,k), zqi_ini, & |
---|
| 855 | zeroklon, znebprecipcld, qincloud, zcf, zeroklon, zeroklon, sursat_e, invtau_e, zqliqth, zqvapclth, zqiceth, & |
---|
| 856 | zficeth, dzficeth,cldfraliqth(:,k), sigma2_icefracturbth(:,k), mean_icefracturbth(:,k)) |
---|
[5614] | 857 | !Environment |
---|
| 858 | DO i=1,klon |
---|
| 859 | qincloud(i)=zqn(i)/(1.-fraca(i,k)) |
---|
| 860 | zcf(i)=MIN(MAX(rneb(i,k)-rnebth(i,k), 0.),1.)*(1.-fraca(i,k)) |
---|
[5633] | 861 | IF (k .GT. 1) THEN |
---|
| 862 | ! evaluate the mixing sursaturation using saturation deficit at level below |
---|
[5647] | 863 | ! as air pacels detraining into clouds have not (less) seen yet entrainement from above |
---|
[5633] | 864 | sursat_e(i)=cloudth_sth(i,k-1)/(zqsith(i)+zdqsith(i)*RCPD/RLSTT*(Tbefthm1(i)-Tbefth(i))) |
---|
[5647] | 865 | ! mixing is assumed to scales with intensity of net detrainment/entrainment rate (D/dz-E/dz) / rho |
---|
| 866 | invtau_e(i)=gamma_mixth*MAX(detr_therm(i,k)-entr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i) |
---|
[5633] | 867 | ELSE |
---|
| 868 | sursat_e(i)=0. |
---|
| 869 | invtau_e(i)=0. |
---|
| 870 | ENDIF |
---|
[5614] | 871 | ENDDO |
---|
[5633] | 872 | CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), zqi_ini, & |
---|
| 873 | ziflcld, znebprecipcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, & |
---|
| 874 | zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k)) |
---|
[5614] | 875 | |
---|
| 876 | ! adjust zfice to account for condensates in thermals'fraction |
---|
| 877 | DO i=1,klon |
---|
| 878 | IF ( zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i) .GT. 0.) THEN |
---|
| 879 | zfice(i)=MIN(1., (zqiceth(i)+zqice(i))/(zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i))) |
---|
| 880 | dzfice(i)=0. ! dxice/dT=0. always when using icefrac_lscp_turb |
---|
| 881 | ENDIF |
---|
| 882 | ENDDO |
---|
| 883 | |
---|
| 884 | ENDIF ! iflag_icefrac |
---|
| 885 | |
---|
| 886 | |
---|
| 887 | |
---|
| 888 | |
---|
| 889 | DO i=1,klon |
---|
| 890 | IF (keepgoing(i)) THEN |
---|
| 891 | |
---|
| 892 | ! P2.2.2> Approximative calculation of temperature variation DT |
---|
| 893 | ! due to condensation. |
---|
| 894 | ! Calculated variables: |
---|
| 895 | ! dT : temperature change due to condensation |
---|
| 896 | !--------------------------------------------------------------- |
---|
| 897 | |
---|
| 898 | |
---|
| 899 | IF (zfice(i).LT.1) THEN |
---|
| 900 | cste=RLVTT |
---|
| 901 | ELSE |
---|
| 902 | cste=RLSTT |
---|
| 903 | ENDIF |
---|
| 904 | |
---|
| 905 | IF ( ok_unadjusted_clouds ) THEN |
---|
| 906 | !--AB We relax the saturation adjustment assumption |
---|
| 907 | !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme |
---|
| 908 | IF ( rneb(i,k) .GT. eps ) THEN |
---|
| 909 | qlibef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k)) |
---|
| 910 | ELSE |
---|
| 911 | qlibef(i) = 0. |
---|
| 912 | ENDIF |
---|
| 913 | ELSE |
---|
| 914 | qlibef(i)=max(0.,zqn(i)-zqs(i)) |
---|
| 915 | ENDIF |
---|
| 916 | |
---|
| 917 | num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT & |
---|
| 918 | +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlibef(i) |
---|
| 919 | denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) & |
---|
| 920 | -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k) & |
---|
| 921 | *qlibef(i)*dzfice(i) |
---|
| 922 | ! here we update a provisory temperature variable that only serves in the iteration |
---|
| 923 | ! process |
---|
| 924 | DT(i)=num/denom |
---|
| 925 | n_i(i)=n_i(i)+1 |
---|
| 926 | |
---|
| 927 | ENDIF ! end keepgoing |
---|
| 928 | |
---|
| 929 | ENDDO ! end loop on i |
---|
| 930 | |
---|
| 931 | ENDDO ! iter=1,iflag_fisrtilp_qsat+1 |
---|
| 932 | |
---|
| 933 | ! P2.2> Final quantities calculation |
---|
| 934 | ! Calculated variables: |
---|
| 935 | ! rneb : cloud fraction |
---|
| 936 | ! zcond: mean condensed water in the mesh |
---|
| 937 | ! zqn : mean water vapor in the mesh |
---|
| 938 | ! zfice: ice fraction in clouds |
---|
| 939 | ! zt : temperature |
---|
| 940 | ! rhcl : clear-sky relative humidity |
---|
| 941 | ! ---------------------------------------------------------------- |
---|
| 942 | |
---|
| 943 | |
---|
| 944 | ! Water vapor and condensed water update, subsequent latent heat exchange for each cloud type |
---|
| 945 | !--------------------------------------------------------------------------------------------- |
---|
| 946 | DO i=1, klon |
---|
| 947 | ! Checks on rneb, rhcl and zqn |
---|
| 948 | IF (rneb(i,k) .LE. 0.0) THEN |
---|
| 949 | zqn(i) = 0.0 |
---|
| 950 | rneb(i,k) = 0.0 |
---|
| 951 | zcond(i) = 0.0 |
---|
| 952 | rhcl(i,k)=zq(i)/zqs(i) |
---|
| 953 | ELSE IF (rneb(i,k) .GE. 1.0) THEN |
---|
| 954 | zqn(i) = zq(i) |
---|
| 955 | rneb(i,k) = 1.0 |
---|
| 956 | IF ( ok_unadjusted_clouds ) THEN |
---|
| 957 | !--AB We relax the saturation adjustment assumption |
---|
| 958 | !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme |
---|
| 959 | zcond(i) = MAX(0., zqn(i) - qvc(i)) |
---|
| 960 | ELSE |
---|
| 961 | zcond(i) = MAX(0.0,zqn(i)-zqs(i)) |
---|
| 962 | ENDIF |
---|
| 963 | rhcl(i,k)=1.0 |
---|
| 964 | ELSE |
---|
| 965 | IF ( ok_unadjusted_clouds ) THEN |
---|
| 966 | !--AB We relax the saturation adjustment assumption |
---|
| 967 | |
---|
| 968 | !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme |
---|
| 969 | zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i)) |
---|
| 970 | ELSE |
---|
| 971 | zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k) |
---|
| 972 | ENDIF |
---|
| 973 | |
---|
| 974 | ! following line is very strange and probably wrong |
---|
| 975 | rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i) |
---|
| 976 | ! Correct calculation of clear-sky relative humidity. To activate once |
---|
| 977 | ! issues related to zqn>zq in bi-gaussian clouds are corrected |
---|
| 978 | !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat |
---|
| 979 | !--This is slighly approximated, the actual formula is |
---|
| 980 | !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q) |
---|
| 981 | !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat |
---|
| 982 | !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content |
---|
| 983 | ! rhcl(i,k) = ( zq(i) - zqn(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i) |
---|
| 984 | ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param) |
---|
| 985 | |
---|
| 986 | ENDIF |
---|
| 987 | |
---|
| 988 | ! water vapor update |
---|
| 989 | zq(i) = zq(i) - zcond(i) |
---|
| 990 | |
---|
| 991 | |
---|
| 992 | ! temperature update due to phase change |
---|
| 993 | zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & |
---|
| 994 | & * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) & |
---|
| 995 | +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) |
---|
| 996 | ENDDO |
---|
| 997 | |
---|
| 998 | ! ---------------------------------------------------------------- |
---|
| 999 | ! P3> Precipitation processes, after cloud formation |
---|
| 1000 | ! - precipitation formation, melting/freezing |
---|
| 1001 | ! ---------------------------------------------------------------- |
---|
| 1002 | |
---|
| 1003 | ! Initiate the quantity of liquid and solid condensates |
---|
| 1004 | ! Note that in the following, zcond is the total amount of condensates |
---|
| 1005 | ! including newly formed precipitations (i.e., condensates formed by the |
---|
| 1006 | ! condensation process), while zoliq is the total amount of condensates in |
---|
| 1007 | ! the cloud (i.e., on which precipitation processes have applied) |
---|
| 1008 | DO i = 1, klon |
---|
| 1009 | zoliq(i) = zcond(i) |
---|
| 1010 | zoliql(i) = zcond(i) * ( 1. - zfice(i) ) |
---|
| 1011 | zoliqi(i) = zcond(i) * zfice(i) |
---|
| 1012 | ENDDO |
---|
| 1013 | |
---|
| 1014 | !================================================================ |
---|
| 1015 | ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R) |
---|
| 1016 | IF (ok_poprecip) THEN |
---|
| 1017 | |
---|
| 1018 | CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), zp, & |
---|
| 1019 | ctot_vol(:,k), ptconv(:,k), & |
---|
| 1020 | zt, zq, zoliql, zoliqi, zfice, & |
---|
| 1021 | rneb(:,k), znebprecipclr, znebprecipcld, & |
---|
| 1022 | zrfl, zrflclr, zrflcld, & |
---|
| 1023 | zifl, ziflclr, ziflcld, & |
---|
| 1024 | qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), & |
---|
| 1025 | dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), & |
---|
| 1026 | dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), & |
---|
| 1027 | dqsmelt(:,k), dqsfreez(:,k) & |
---|
| 1028 | ) |
---|
| 1029 | DO i = 1, klon |
---|
| 1030 | zoliq(i) = zoliql(i) + zoliqi(i) |
---|
| 1031 | ENDDO |
---|
| 1032 | |
---|
| 1033 | !================================================================ |
---|
| 1034 | ELSE |
---|
| 1035 | |
---|
| 1036 | CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, & |
---|
| 1037 | ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, & |
---|
| 1038 | zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, & |
---|
| 1039 | rneb(:,k), znebprecipclr, znebprecipcld, & |
---|
| 1040 | zneb, tot_zneb, zrho_up, zvelo_up, & |
---|
| 1041 | zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, & |
---|
| 1042 | zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) & |
---|
| 1043 | ) |
---|
| 1044 | |
---|
| 1045 | ENDIF ! ok_poprecip |
---|
| 1046 | |
---|
| 1047 | ! End of precipitation processes after cloud formation |
---|
| 1048 | ! ---------------------------------------------------- |
---|
| 1049 | |
---|
| 1050 | !---------------------------------------------------------------------- |
---|
| 1051 | ! P4> Calculation of cloud condensates amount seen by radiative scheme |
---|
| 1052 | !---------------------------------------------------------------------- |
---|
| 1053 | |
---|
| 1054 | DO i=1,klon |
---|
| 1055 | |
---|
| 1056 | IF (ok_poprecip) THEN |
---|
| 1057 | IF (ok_radocond_snow) THEN |
---|
| 1058 | radocond(i,k) = zoliq(i) |
---|
| 1059 | zradoice(i) = zoliqi(i) + qsnowdiag(i,k) |
---|
| 1060 | ELSE |
---|
| 1061 | radocond(i,k) = zoliq(i) |
---|
| 1062 | zradoice(i) = zoliqi(i) |
---|
| 1063 | ENDIF |
---|
| 1064 | ELSE |
---|
| 1065 | radocond(i,k) = zradocond(i) |
---|
| 1066 | ENDIF |
---|
| 1067 | |
---|
| 1068 | ! calculate the percentage of ice in "radocond" so cloud+precip seen |
---|
| 1069 | ! by the radiation scheme |
---|
| 1070 | IF (radocond(i,k) .GT. 0.) THEN |
---|
| 1071 | radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.) |
---|
| 1072 | ENDIF |
---|
| 1073 | ENDDO |
---|
| 1074 | |
---|
| 1075 | ! ---------------------------------------------------------------- |
---|
| 1076 | ! P5> Wet scavenging |
---|
| 1077 | ! ---------------------------------------------------------------- |
---|
| 1078 | |
---|
| 1079 | !Scavenging through nucleation in the layer |
---|
| 1080 | |
---|
| 1081 | DO i = 1,klon |
---|
| 1082 | |
---|
| 1083 | IF(zcond(i).GT.zoliq(i)+1.e-10) THEN |
---|
| 1084 | beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime |
---|
| 1085 | ELSE |
---|
| 1086 | beta(i,k) = 0. |
---|
| 1087 | ENDIF |
---|
| 1088 | |
---|
| 1089 | zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG |
---|
| 1090 | |
---|
| 1091 | IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN |
---|
| 1092 | |
---|
| 1093 | IF (temp(i,k) .GE. temp_nowater) THEN |
---|
| 1094 | zalpha_tr = a_tr_sca(3) |
---|
| 1095 | ELSE |
---|
| 1096 | zalpha_tr = a_tr_sca(4) |
---|
| 1097 | ENDIF |
---|
| 1098 | |
---|
| 1099 | zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb)) |
---|
| 1100 | frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi |
---|
| 1101 | |
---|
| 1102 | ! Nucleation with a factor of -1 instead of -0.5 |
---|
| 1103 | zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb)) |
---|
| 1104 | |
---|
| 1105 | ENDIF |
---|
| 1106 | |
---|
| 1107 | ENDDO |
---|
| 1108 | |
---|
| 1109 | ! Scavenging through impaction in the underlying layer |
---|
| 1110 | |
---|
| 1111 | DO kk = k-1, 1, -1 |
---|
| 1112 | |
---|
| 1113 | DO i = 1, klon |
---|
| 1114 | |
---|
| 1115 | IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN |
---|
| 1116 | |
---|
| 1117 | IF (temp(i,kk) .GE. temp_nowater) THEN |
---|
| 1118 | zalpha_tr = a_tr_sca(1) |
---|
| 1119 | ELSE |
---|
| 1120 | zalpha_tr = a_tr_sca(2) |
---|
| 1121 | ENDIF |
---|
| 1122 | |
---|
| 1123 | zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb)) |
---|
| 1124 | frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi |
---|
| 1125 | |
---|
| 1126 | ENDIF |
---|
| 1127 | |
---|
| 1128 | ENDDO |
---|
| 1129 | |
---|
| 1130 | ENDDO |
---|
| 1131 | |
---|
| 1132 | !------------------------------------------------------------ |
---|
| 1133 | ! P6 > write diagnostics and outputs |
---|
| 1134 | !------------------------------------------------------------ |
---|
| 1135 | |
---|
| 1136 | !--AB Write diagnostics and tracers for ice supersaturation |
---|
| 1137 | IF ( ok_ice_supersat ) THEN |
---|
| 1138 | CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,1,.false.,qsatl(:,k),zdqs) |
---|
| 1139 | CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,2,.false.,qsati(:,k),zdqs) |
---|
| 1140 | |
---|
| 1141 | DO i = 1, klon |
---|
| 1142 | |
---|
| 1143 | IF ( zoliq(i) .LE. 0. ) THEN |
---|
| 1144 | !--If everything was precipitated, the remaining empty cloud is dissipated |
---|
| 1145 | !--and everything is transfered to the subsaturated clear sky region |
---|
| 1146 | rneb(i,k) = 0. |
---|
| 1147 | qvc(i) = 0. |
---|
| 1148 | ENDIF |
---|
| 1149 | |
---|
| 1150 | cf_seri(i,k) = rneb(i,k) |
---|
| 1151 | |
---|
| 1152 | IF ( .NOT. ok_unadjusted_clouds ) THEN |
---|
| 1153 | qvc(i) = zqs(i) * rneb(i,k) |
---|
| 1154 | ENDIF |
---|
| 1155 | IF ( zq(i) .GT. min_qParent ) THEN |
---|
| 1156 | rvc_seri(i,k) = qvc(i) / zq(i) |
---|
| 1157 | ELSE |
---|
| 1158 | rvc_seri(i,k) = min_ratio |
---|
| 1159 | ENDIF |
---|
| 1160 | !--The MIN barrier is NEEDED because of: |
---|
| 1161 | !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes) |
---|
| 1162 | !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot) |
---|
| 1163 | !--The MAX barrier is a safeguard that should not be activated |
---|
| 1164 | rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.) |
---|
| 1165 | |
---|
| 1166 | !--Diagnostics |
---|
| 1167 | gamma_cond(i,k) = gammasat(i) |
---|
| 1168 | subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k) |
---|
| 1169 | qsub(i,k) = zq(i) - qvc(i) - qissr(i,k) |
---|
| 1170 | qcld(i,k) = qvc(i) + zoliq(i) |
---|
| 1171 | ENDDO |
---|
| 1172 | ENDIF |
---|
| 1173 | |
---|
| 1174 | ! Outputs: |
---|
| 1175 | !------------------------------- |
---|
| 1176 | ! Precipitation fluxes at layer interfaces |
---|
| 1177 | ! + precipitation fractions + |
---|
| 1178 | ! temperature and water species tendencies |
---|
| 1179 | DO i = 1, klon |
---|
| 1180 | psfl(i,k)=zifl(i) |
---|
| 1181 | prfl(i,k)=zrfl(i) |
---|
| 1182 | pfraclr(i,k)=znebprecipclr(i) |
---|
| 1183 | pfracld(i,k)=znebprecipcld(i) |
---|
| 1184 | distcltop(i,k)=zdistcltop(i) |
---|
| 1185 | temp_cltop(i,k)=ztemp_cltop(i) |
---|
| 1186 | d_q(i,k) = zq(i) - qt(i,k) |
---|
| 1187 | d_t(i,k) = zt(i) - temp(i,k) |
---|
| 1188 | |
---|
| 1189 | IF (ok_bug_phase_lscp) THEN |
---|
| 1190 | d_ql(i,k) = (1-zfice(i))*zoliq(i) |
---|
| 1191 | d_qi(i,k) = zfice(i)*zoliq(i) |
---|
| 1192 | ELSE |
---|
| 1193 | d_ql(i,k) = zoliql(i) |
---|
| 1194 | d_qi(i,k) = zoliqi(i) |
---|
| 1195 | ENDIF |
---|
| 1196 | |
---|
| 1197 | ENDDO |
---|
| 1198 | |
---|
| 1199 | |
---|
| 1200 | ENDDO ! loop on k from top to bottom |
---|
| 1201 | |
---|
| 1202 | |
---|
| 1203 | ! Rain or snow at the surface (depending on the first layer temperature) |
---|
| 1204 | DO i = 1, klon |
---|
| 1205 | snow(i) = zifl(i) |
---|
| 1206 | rain(i) = zrfl(i) |
---|
| 1207 | ! c_iso final output |
---|
| 1208 | ENDDO |
---|
| 1209 | |
---|
| 1210 | IF (ncoreczq>0) THEN |
---|
| 1211 | WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.' |
---|
| 1212 | ENDIF |
---|
| 1213 | |
---|
| 1214 | |
---|
| 1215 | RETURN |
---|
| 1216 | |
---|
| 1217 | END SUBROUTINE lscp |
---|
| 1218 | !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 1219 | |
---|
| 1220 | END MODULE lmdz_lscp_main |
---|