[4803] | 1 | MODULE lmdz_lscp_poprecip |
---|
| 2 | !---------------------------------------------------------------- |
---|
| 3 | ! Module for the process-oriented treament of precipitation |
---|
| 4 | ! that are called in LSCP |
---|
| 5 | ! Authors: Atelier Nuage (G. Riviere, L. Raillard, M. Wimmer, |
---|
| 6 | ! N. Dutrievoz, E. Vignon, A. Borella, et al.) |
---|
| 7 | ! Jan. 2024 |
---|
| 8 | |
---|
| 9 | |
---|
| 10 | IMPLICIT NONE |
---|
| 11 | |
---|
| 12 | CONTAINS |
---|
| 13 | |
---|
| 14 | !---------------------------------------------------------------- |
---|
| 15 | ! Computes the processes-oriented precipitation formulations for |
---|
| 16 | ! evaporation and sublimation |
---|
| 17 | ! |
---|
[4879] | 18 | SUBROUTINE poprecip_precld( & |
---|
[4803] | 19 | klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, & |
---|
| 20 | qprecip, precipfracclr, precipfraccld, & |
---|
| 21 | rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub & |
---|
| 22 | ) |
---|
| 23 | |
---|
| 24 | USE lmdz_lscp_ini, ONLY : prt_level, lunout |
---|
[4830] | 25 | USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac |
---|
[4803] | 26 | USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG |
---|
| 27 | USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf |
---|
| 28 | |
---|
| 29 | IMPLICIT NONE |
---|
| 30 | |
---|
| 31 | |
---|
| 32 | INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-] |
---|
| 33 | REAL, INTENT(IN) :: dtime !--time step [s] |
---|
| 34 | LOGICAL, INTENT(IN) :: iftop !--if top of the column |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa] |
---|
| 38 | REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa] |
---|
| 39 | REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa] |
---|
| 40 | |
---|
| 41 | REAL, INTENT(INOUT), DIMENSION(klon) :: temp !--current temperature [K] |
---|
| 42 | REAL, INTENT(INOUT), DIMENSION(klon) :: tempupnew !--updated temperature of the overlying layer [K] |
---|
| 43 | |
---|
| 44 | REAL, INTENT(INOUT), DIMENSION(klon) :: qvap !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg] |
---|
| 45 | REAL, INTENT(INOUT), DIMENSION(klon) :: qprecip !--specific humidity in the precipitation falling from the upper layer [kg/kg] |
---|
| 46 | |
---|
| 47 | REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-] |
---|
| 48 | REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-] |
---|
| 49 | |
---|
[4809] | 50 | REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2] |
---|
| 51 | REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2] |
---|
| 52 | REAL, INTENT(IN), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2] |
---|
| 53 | REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2] |
---|
| 54 | REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2] |
---|
| 55 | REAL, INTENT(IN), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2] |
---|
[4803] | 56 | |
---|
| 57 | REAL, INTENT(OUT), DIMENSION(klon) :: dqreva !--rain tendency due to evaporation [kg/kg/s] |
---|
| 58 | REAL, INTENT(OUT), DIMENSION(klon) :: dqssub !--snow tendency due to sublimation [kg/kg/s] |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | |
---|
[4819] | 63 | !--Integer for interating over klon |
---|
[4803] | 64 | INTEGER :: i |
---|
[4833] | 65 | !--dhum_to_dflux: coef to convert a specific quantity to a flux |
---|
| 66 | REAL, DIMENSION(klon) :: dhum_to_dflux |
---|
| 67 | !-- |
---|
| 68 | REAL, DIMENSION(klon) :: rho, dz |
---|
[4803] | 69 | |
---|
[4819] | 70 | !--Saturation values |
---|
[4803] | 71 | REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati |
---|
[4819] | 72 | !--Fluxes tendencies because of evaporation and sublimation |
---|
| 73 | REAL :: dprecip_evasub_max, draineva, dsnowsub, dprecip_evasub_tot |
---|
| 74 | !--Specific humidity tendencies because of evaporation and sublimation |
---|
| 75 | REAL :: dqrevap, dqssubl |
---|
| 76 | !--Specific heat constant |
---|
[4803] | 77 | REAL :: cpair, cpw |
---|
| 78 | |
---|
[4819] | 79 | !--Initialisation |
---|
| 80 | qzero(:) = 0. |
---|
| 81 | dqreva(:) = 0. |
---|
| 82 | dqssub(:) = 0. |
---|
| 83 | dqrevap = 0. |
---|
| 84 | dqssubl = 0. |
---|
[4803] | 85 | |
---|
[4833] | 86 | !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt |
---|
| 87 | dhum_to_dflux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime |
---|
| 88 | rho(:) = pplay(:) / temp(:) / RD |
---|
| 89 | dz(:) = ( paprsdn(:) - paprsup(:) ) / RG / rho(:) |
---|
[4819] | 90 | |
---|
| 91 | !--Calculation of saturation specific humidity |
---|
| 92 | !--depending on temperature: |
---|
[4803] | 93 | CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:)) |
---|
[4819] | 94 | !--wrt liquid water |
---|
[4803] | 95 | CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:)) |
---|
[4819] | 96 | !--wrt ice |
---|
[4803] | 97 | CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:)) |
---|
| 98 | |
---|
| 99 | |
---|
| 100 | |
---|
[4819] | 101 | !--First step consists in "thermalizing" the layer: |
---|
| 102 | !--as the flux of precip from layer above "advects" some heat (as the precip is at the temperature |
---|
| 103 | !--of the overlying layer) we recalculate a mean temperature that both the air and the precip in the |
---|
| 104 | !--layer have. |
---|
[4803] | 105 | |
---|
| 106 | IF (iftop) THEN |
---|
| 107 | |
---|
[4819] | 108 | DO i = 1, klon |
---|
| 109 | qprecip(i) = 0. |
---|
| 110 | ENDDO |
---|
[4803] | 111 | |
---|
| 112 | ELSE |
---|
| 113 | |
---|
[4819] | 114 | DO i = 1, klon |
---|
| 115 | !--No condensed water so cp=cp(vapor+dry air) |
---|
| 116 | !-- RVTMP2=rcpv/rcpd-1 |
---|
| 117 | cpair = RCPD * ( 1. + RVTMP2 * qvap(i) ) |
---|
| 118 | cpw = RCPD * RVTMP2 |
---|
| 119 | !--qprecip has to be thermalized with |
---|
| 120 | !--layer's air so that precipitation at the ground has the |
---|
| 121 | !--same temperature as the lowermost layer |
---|
| 122 | !--we convert the flux into a specific quantity qprecip |
---|
[4833] | 123 | qprecip(i) = ( rain(i) + snow(i) ) / dhum_to_dflux(i) |
---|
[4819] | 124 | !-- t(i,k+1) + d_t(i,k+1): new temperature of the overlying layer |
---|
| 125 | temp(i) = ( tempupnew(i) * qprecip(i) * cpw + cpair * temp(i) ) & |
---|
| 126 | / ( cpair + qprecip(i) * cpw ) |
---|
| 127 | ENDDO |
---|
[4803] | 128 | |
---|
| 129 | ENDIF |
---|
| 130 | |
---|
| 131 | |
---|
| 132 | DO i = 1, klon |
---|
| 133 | |
---|
[4819] | 134 | !--If there is precipitation from the layer above |
---|
[4803] | 135 | IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN |
---|
| 136 | |
---|
[4819] | 137 | !--Evaporation of liquid precipitation coming from above |
---|
| 138 | !--in the clear sky only |
---|
[4833] | 139 | !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva) |
---|
[4819] | 140 | !--formula from Sundqvist 1988, Klemp & Wilhemson 1978 |
---|
[4833] | 141 | !--Explicit formulation |
---|
| 142 | !draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) * dz(i) & |
---|
| 143 | ! * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva & |
---|
| 144 | !draineva = MAX( - rainclr(i), draineva) |
---|
| 145 | !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly) |
---|
| 146 | !--which does not need a barrier on rainclr, because included in the formula |
---|
| 147 | draineva = precipfracclr(i) * ( MAX(0., & |
---|
| 148 | - coef_eva * ( 1. - expo_eva ) * (1. - qvap(i) / qsatl(i)) * dz(i) & |
---|
| 149 | + ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_eva ) & |
---|
| 150 | ) )**( 1. / ( 1. - expo_eva ) ) - rainclr(i) |
---|
| 151 | |
---|
| 152 | !--Evaporation is limited by 0 |
---|
| 153 | draineva = MIN(0., draineva) |
---|
[4803] | 154 | |
---|
| 155 | |
---|
[4819] | 156 | !--Sublimation of the solid precipitation coming from above |
---|
| 157 | !--(same formula as for liquid precip) |
---|
[4833] | 158 | !--Explicit formulation |
---|
| 159 | !dsnowsub = - precipfracclr(i) * coef_sub * (1. - qvap(i) / qsatl(i)) * dz(i) & |
---|
| 160 | ! * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_sub & |
---|
| 161 | !dsnowsub = MAX( - snowclr(i), dsnowsub) |
---|
| 162 | !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly) |
---|
| 163 | !--which does not need a barrier on snowclr, because included in the formula |
---|
| 164 | dsnowsub = precipfracclr(i) * ( MAX(0., & |
---|
[4878] | 165 | - coef_sub * ( 1. - expo_sub ) * (1. - qvap(i) / qsati(i)) * dz(i) & |
---|
[4833] | 166 | + ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_sub ) & |
---|
| 167 | ) )**( 1. / ( 1. - expo_sub ) ) - snowclr(i) |
---|
[4803] | 168 | |
---|
[4833] | 169 | !--Sublimation is limited by 0 |
---|
[4803] | 170 | ! TODO: change max when we will allow for vapor deposition in supersaturated regions |
---|
[4833] | 171 | dsnowsub = MIN(0., dsnowsub) |
---|
[4803] | 172 | |
---|
[4819] | 173 | !--Evaporation limit: we ensure that the layer's fraction below |
---|
| 174 | !--the clear sky does not reach saturation. In this case, we |
---|
| 175 | !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice |
---|
| 176 | !--Max evaporation is computed not to saturate the clear sky precip fraction |
---|
| 177 | !--(i.e., the fraction where evaporation occurs) |
---|
| 178 | !--It is expressed as a max flux dprecip_evasub_max |
---|
| 179 | |
---|
| 180 | dprecip_evasub_max = MIN(0., ( qvap(i) - qsat(i) ) * precipfracclr(i)) & |
---|
[4833] | 181 | * dhum_to_dflux(i) |
---|
[4819] | 182 | dprecip_evasub_tot = draineva + dsnowsub |
---|
[4803] | 183 | |
---|
[4819] | 184 | !--Barriers |
---|
| 185 | !--If activates if the total is LOWER than the max because |
---|
| 186 | !--everything is negative |
---|
| 187 | IF ( dprecip_evasub_tot .LT. dprecip_evasub_max ) THEN |
---|
| 188 | draineva = dprecip_evasub_max * draineva / dprecip_evasub_tot |
---|
| 189 | dsnowsub = dprecip_evasub_max * dsnowsub / dprecip_evasub_tot |
---|
[4803] | 190 | ENDIF |
---|
| 191 | |
---|
| 192 | |
---|
[4819] | 193 | !--New solid and liquid precipitation fluxes after evap and sublimation |
---|
[4833] | 194 | dqrevap = draineva / dhum_to_dflux(i) |
---|
| 195 | dqssubl = dsnowsub / dhum_to_dflux(i) |
---|
[4803] | 196 | |
---|
| 197 | |
---|
[4819] | 198 | !--Vapor is updated after evaporation/sublimation (it is increased) |
---|
| 199 | qvap(i) = qvap(i) - dqrevap - dqssubl |
---|
| 200 | !--qprecip is the total condensed water in the precip flux (it is decreased) |
---|
| 201 | qprecip(i) = qprecip(i) + dqrevap + dqssubl |
---|
| 202 | !--Air and precip temperature (i.e., gridbox temperature) |
---|
| 203 | !--is updated due to latent heat cooling |
---|
[4803] | 204 | temp(i) = temp(i) & |
---|
[4819] | 205 | + dqrevap * RLVTT / RCPD & |
---|
| 206 | / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) & |
---|
| 207 | + dqssubl * RLSTT / RCPD & |
---|
| 208 | / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) |
---|
[4803] | 209 | |
---|
[4819] | 210 | !--Add tendencies |
---|
[4833] | 211 | !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) |
---|
[4832] | 212 | rainclr(i) = MAX(0., rainclr(i) + draineva) |
---|
| 213 | snowclr(i) = MAX(0., snowclr(i) + dsnowsub) |
---|
| 214 | |
---|
[4819] | 215 | !--If there is no more precip fluxes, the precipitation fraction in clear |
---|
| 216 | !--sky is set to 0 |
---|
[4803] | 217 | IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0. |
---|
| 218 | |
---|
[4819] | 219 | !--Calculation of the total fluxes |
---|
[4803] | 220 | rain(i) = rainclr(i) + raincld(i) |
---|
| 221 | snow(i) = snowclr(i) + snowcld(i) |
---|
| 222 | |
---|
| 223 | ELSE |
---|
[4819] | 224 | !--If no precip, we reinitialize the cloud fraction used for the precip to 0 |
---|
[4803] | 225 | precipfraccld(i) = 0. |
---|
| 226 | precipfracclr(i) = 0. |
---|
| 227 | |
---|
| 228 | ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. ) |
---|
| 229 | |
---|
[4819] | 230 | !--Diagnostic tendencies |
---|
| 231 | dqssub(i) = dqssubl / dtime |
---|
| 232 | dqreva(i) = dqrevap / dtime |
---|
[4803] | 233 | |
---|
| 234 | ENDDO ! loop on klon |
---|
| 235 | |
---|
| 236 | |
---|
[4880] | 237 | END SUBROUTINE poprecip_precld |
---|
[4803] | 238 | |
---|
| 239 | !---------------------------------------------------------------- |
---|
| 240 | ! Computes the processes-oriented precipitation formulations for |
---|
| 241 | ! - autoconversion (auto) via a deposition process |
---|
| 242 | ! - aggregation (agg) |
---|
| 243 | ! - riming (rim) |
---|
[4830] | 244 | ! - collection (col) |
---|
[4803] | 245 | ! - melting (melt) |
---|
[4830] | 246 | ! - freezing (freez) |
---|
[4803] | 247 | ! |
---|
| 248 | SUBROUTINE poprecip_postcld( & |
---|
| 249 | klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, & |
---|
| 250 | temp, qvap, qliq, qice, icefrac, cldfra, & |
---|
| 251 | precipfracclr, precipfraccld, & |
---|
| 252 | rain, rainclr, raincld, snow, snowclr, snowcld, & |
---|
[4830] | 253 | qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, & |
---|
[4819] | 254 | dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez) |
---|
[4803] | 255 | |
---|
| 256 | USE lmdz_lscp_ini, ONLY : prt_level, lunout |
---|
[4818] | 257 | USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI |
---|
[4803] | 258 | USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf |
---|
| 259 | |
---|
| 260 | USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb, & |
---|
| 261 | cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, & |
---|
| 262 | thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, & |
---|
[4832] | 263 | rho_rain, rho_snow, r_rain, r_snow, rho_ice, r_ice, & |
---|
[4818] | 264 | tau_auto_snow_min, tau_auto_snow_max, & |
---|
[4882] | 265 | thresh_precip_frac, eps, & |
---|
[4830] | 266 | coef_ventil, alpha_freez, beta_freez, temp_nowater, & |
---|
| 267 | iflag_cloudth_vert, iflag_rain_incloud_vol, & |
---|
| 268 | cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez, & |
---|
| 269 | rain_fallspeed_clr, rain_fallspeed_cld, & |
---|
| 270 | snow_fallspeed_clr, snow_fallspeed_cld |
---|
[4803] | 271 | |
---|
[4818] | 272 | |
---|
[4803] | 273 | IMPLICIT NONE |
---|
| 274 | |
---|
| 275 | INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-] |
---|
| 276 | REAL, INTENT(IN) :: dtime !--time step [s] |
---|
| 277 | |
---|
| 278 | REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa] |
---|
| 279 | REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa] |
---|
| 280 | REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa] |
---|
| 281 | |
---|
[4819] | 282 | REAL, INTENT(IN), DIMENSION(klon) :: ctot_vol !--volumic cloud fraction [-] |
---|
| 283 | LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv !--true if we are in a convective point |
---|
[4803] | 284 | |
---|
| 285 | REAL, INTENT(INOUT), DIMENSION(klon) :: temp !--current temperature [K] |
---|
| 286 | REAL, INTENT(INOUT), DIMENSION(klon) :: qvap !--current water vapor specific humidity [kg/kg] |
---|
| 287 | REAL, INTENT(INOUT), DIMENSION(klon) :: qliq !--current liquid water specific humidity [kg/kg] |
---|
| 288 | REAL, INTENT(INOUT), DIMENSION(klon) :: qice !--current ice water specific humidity [kg/kg] |
---|
[4819] | 289 | REAL, INTENT(IN), DIMENSION(klon) :: icefrac !--ice fraction [-] |
---|
| 290 | REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction [-] |
---|
[4803] | 291 | |
---|
| 292 | REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-] |
---|
| 293 | REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-] |
---|
| 294 | !--NB. at the end of the routine, becomes the fraction of precip |
---|
| 295 | !--in the current layer |
---|
| 296 | |
---|
[4809] | 297 | REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2] |
---|
| 298 | REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2] |
---|
| 299 | REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2] |
---|
| 300 | REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2] |
---|
| 301 | REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2] |
---|
| 302 | REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2] |
---|
[4803] | 303 | |
---|
[4830] | 304 | REAL, INTENT(OUT), DIMENSION(klon) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg] |
---|
| 305 | REAL, INTENT(OUT), DIMENSION(klon) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg] |
---|
[4819] | 306 | REAL, INTENT(OUT), DIMENSION(klon) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s] |
---|
| 307 | REAL, INTENT(OUT), DIMENSION(klon) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s] |
---|
| 308 | REAL, INTENT(OUT), DIMENSION(klon) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s] |
---|
| 309 | REAL, INTENT(OUT), DIMENSION(klon) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s] |
---|
| 310 | REAL, INTENT(OUT), DIMENSION(klon) :: dqsrim !--snow tendency due to riming [kg/kg/s] |
---|
| 311 | REAL, INTENT(OUT), DIMENSION(klon) :: dqsmelt !--snow tendency due to melting [kg/kg/s] |
---|
| 312 | REAL, INTENT(OUT), DIMENSION(klon) :: dqrmelt !--rain tendency due to melting [kg/kg/s] |
---|
| 313 | REAL, INTENT(OUT), DIMENSION(klon) :: dqsfreez !--snow tendency due to freezing [kg/kg/s] |
---|
| 314 | REAL, INTENT(OUT), DIMENSION(klon) :: dqrfreez !--rain tendency due to freezing [kg/kg/s] |
---|
[4803] | 315 | |
---|
| 316 | |
---|
| 317 | |
---|
| 318 | !--Local variables |
---|
| 319 | |
---|
| 320 | INTEGER :: i |
---|
[4833] | 321 | REAL, DIMENSION(klon) :: dhum_to_dflux |
---|
[4818] | 322 | REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip |
---|
[4803] | 323 | |
---|
[4830] | 324 | !--Partition of the fluxes |
---|
[4803] | 325 | REAL :: dcldfra |
---|
| 326 | REAL :: precipfractot |
---|
| 327 | REAL :: dprecipfracclr, dprecipfraccld |
---|
| 328 | REAL :: drainclr, dsnowclr |
---|
| 329 | REAL :: draincld, dsnowcld |
---|
[4818] | 330 | |
---|
[4830] | 331 | !--Collection, aggregation and riming |
---|
[4803] | 332 | REAL :: eff_cldfra |
---|
[4819] | 333 | REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp |
---|
[4818] | 334 | REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq |
---|
[4830] | 335 | REAL :: dqlcol !--loss of liquid cloud content due to collection by rain [kg/kg/s] |
---|
| 336 | REAL :: dqiagg !--loss of ice cloud content due to collection by aggregation [kg/kg/s] |
---|
| 337 | REAL :: dqlrim !--loss of liquid cloud content due to riming on snow [kg/kg/s] |
---|
[4818] | 338 | |
---|
[4830] | 339 | !--Autoconversion |
---|
[4803] | 340 | REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain |
---|
| 341 | REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow |
---|
[4830] | 342 | REAL :: dqlauto !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s] |
---|
| 343 | REAL :: dqiauto !--loss of ice cloud content due to autoconversion to snow [kg/kg/s] |
---|
[4803] | 344 | |
---|
[4830] | 345 | !--Melting |
---|
[4882] | 346 | REAL :: dqsmelt_max, air_thermal_conduct |
---|
[4818] | 347 | REAL :: nb_snowflake_clr, nb_snowflake_cld |
---|
| 348 | REAL :: capa_snowflake, temp_wetbulb |
---|
| 349 | REAL :: dqsclrmelt, dqscldmelt, dqstotmelt |
---|
| 350 | REAL, DIMENSION(klon) :: qzero, qsat, dqsat |
---|
[4803] | 351 | |
---|
[4830] | 352 | !--Freezing |
---|
[4818] | 353 | REAL :: dqrfreez_max |
---|
| 354 | REAL :: tau_freez |
---|
[4832] | 355 | REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2 |
---|
[4830] | 356 | REAL :: coef_freez |
---|
| 357 | REAL :: dqifreez !--loss of ice cloud content due to collection of ice from rain [kg/kg/s] |
---|
| 358 | REAL :: Eff_rain_ice |
---|
[4818] | 359 | |
---|
| 360 | |
---|
[4803] | 361 | !--Initialisation of variables |
---|
| 362 | |
---|
[4888] | 363 | |
---|
[4818] | 364 | qzero(:) = 0. |
---|
| 365 | |
---|
| 366 | dqrcol(:) = 0. |
---|
| 367 | dqsagg(:) = 0. |
---|
| 368 | dqsauto(:) = 0. |
---|
| 369 | dqrauto(:) = 0. |
---|
| 370 | dqsrim(:) = 0. |
---|
| 371 | dqrmelt(:) = 0. |
---|
| 372 | dqsmelt(:) = 0. |
---|
[4803] | 373 | dqrfreez(:) = 0. |
---|
| 374 | dqsfreez(:) = 0. |
---|
| 375 | |
---|
| 376 | |
---|
[4809] | 377 | DO i = 1, klon |
---|
[4803] | 378 | |
---|
[4830] | 379 | !--Variables initialisation |
---|
[4818] | 380 | dqlcol = 0. |
---|
| 381 | dqiagg = 0. |
---|
| 382 | dqiauto = 0. |
---|
| 383 | dqlauto = 0. |
---|
[4830] | 384 | dqlrim = 0. |
---|
[4803] | 385 | |
---|
[4833] | 386 | !--dhum_to_dflux: coef to convert a delta in specific quantity to a flux |
---|
| 387 | !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt |
---|
| 388 | dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime |
---|
[4818] | 389 | qtot(i) = qvap(i) + qliq(i) + qice(i) & |
---|
[4833] | 390 | + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i) |
---|
[4818] | 391 | |
---|
[4809] | 392 | !------------------------------------------------------------ |
---|
| 393 | !-- PRECIPITATION FRACTIONS UPDATE |
---|
| 394 | !------------------------------------------------------------ |
---|
| 395 | !--The goal of this routine is to reattribute precipitation fractions |
---|
| 396 | !--and fluxes to clear or cloudy air, depending on the variation of |
---|
| 397 | !--the cloud fraction on the vertical dimension. We assume a |
---|
| 398 | !--maximum-random overlap of the cloud cover (see Jakob and Klein, 2000, |
---|
| 399 | !--and LTP thesis, 2021) |
---|
| 400 | !--NB. in fact, we assume a maximum-random overlap of the total precip. frac |
---|
[4803] | 401 | |
---|
| 402 | !--Initialisation |
---|
| 403 | precipfractot = precipfracclr(i) + precipfraccld(i) |
---|
| 404 | |
---|
[4809] | 405 | !--Instead of using the cloud cover which was use in LTP thesis, we use the |
---|
| 406 | !--total precip. fraction to compute the maximum-random overlap. This is |
---|
| 407 | !--because all the information of the cloud cover is embedded into |
---|
| 408 | !--precipfractot, and this allows for taking into account the potential |
---|
| 409 | !--reduction of the precipitation fraction because either the flux is too |
---|
| 410 | !--small (see barrier at the end of poprecip_postcld) or the flux is completely |
---|
| 411 | !--evaporated (see barrier at the end of poprecip_precld) |
---|
| 412 | !--NB. precipfraccld(i) is here the cloud fraction of the layer above |
---|
[4832] | 413 | !precipfractot = 1. - ( 1. - precipfractot ) * & |
---|
| 414 | ! ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) & |
---|
| 415 | ! / ( 1. - MIN( precipfraccld(i), 1. - eps ) ) |
---|
| 416 | |
---|
| 417 | |
---|
| 418 | IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN |
---|
| 419 | precipfractot = 1. |
---|
| 420 | ELSE |
---|
| 421 | precipfractot = 1. - ( 1. - precipfractot ) * & |
---|
[4803] | 422 | ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) & |
---|
[4832] | 423 | / ( 1. - precipfraccld(i) ) |
---|
| 424 | ENDIF |
---|
[4803] | 425 | |
---|
[4809] | 426 | !--precipfraccld(i) is here the cloud fraction of the layer above |
---|
[4803] | 427 | dcldfra = cldfra(i) - precipfraccld(i) |
---|
[4809] | 428 | !--Tendency of the clear-sky precipitation fraction. We add a MAX on the |
---|
| 429 | !--calculation of the current CS precip. frac. |
---|
[4830] | 430 | !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i) |
---|
| 431 | !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated |
---|
| 432 | !--if precipfractot < cldfra) |
---|
| 433 | dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i) |
---|
[4809] | 434 | !--Tendency of the cloudy precipitation fraction. We add a MAX on the |
---|
| 435 | !--calculation of the current CS precip. frac. |
---|
| 436 | !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) ) |
---|
[4830] | 437 | !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated |
---|
[4809] | 438 | !--if cldfra < 0) |
---|
| 439 | dprecipfraccld = dcldfra |
---|
[4803] | 440 | |
---|
| 441 | |
---|
[4809] | 442 | !--If the cloud extends |
---|
[4803] | 443 | IF ( dprecipfraccld .GT. 0. ) THEN |
---|
[4809] | 444 | !--If there is no CS precip, nothing happens. |
---|
| 445 | !--If there is, we reattribute some of the CS precip flux |
---|
| 446 | !--to the cloud precip flux, proportionnally to the |
---|
| 447 | !--decrease of the CS precip fraction |
---|
| 448 | IF ( precipfracclr(i) .LE. 0. ) THEN |
---|
| 449 | drainclr = 0. |
---|
| 450 | dsnowclr = 0. |
---|
| 451 | ELSE |
---|
[4803] | 452 | drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i) |
---|
| 453 | dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i) |
---|
| 454 | ENDIF |
---|
[4809] | 455 | !--If the cloud narrows |
---|
| 456 | ELSEIF ( dprecipfraccld .LT. 0. ) THEN |
---|
| 457 | !--We reattribute some of the cloudy precip flux |
---|
| 458 | !--to the CS precip flux, proportionnally to the |
---|
| 459 | !--decrease of the cloud precip fraction |
---|
| 460 | draincld = dprecipfraccld / precipfraccld(i) * raincld(i) |
---|
| 461 | dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i) |
---|
| 462 | drainclr = - draincld |
---|
| 463 | dsnowclr = - dsnowcld |
---|
| 464 | !--If the cloud stays the same or if there is no cloud above and |
---|
| 465 | !--in the current layer, nothing happens |
---|
[4803] | 466 | ELSE |
---|
[4809] | 467 | drainclr = 0. |
---|
| 468 | dsnowclr = 0. |
---|
[4803] | 469 | ENDIF |
---|
| 470 | |
---|
[4809] | 471 | !--We add the tendencies |
---|
[4833] | 472 | !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) |
---|
[4803] | 473 | precipfraccld(i) = precipfraccld(i) + dprecipfraccld |
---|
| 474 | precipfracclr(i) = precipfracclr(i) + dprecipfracclr |
---|
[4833] | 475 | rainclr(i) = MAX(0., rainclr(i) + drainclr) |
---|
| 476 | snowclr(i) = MAX(0., snowclr(i) + dsnowclr) |
---|
| 477 | raincld(i) = MAX(0., raincld(i) - drainclr) |
---|
| 478 | snowcld(i) = MAX(0., snowcld(i) - dsnowclr) |
---|
[4803] | 479 | |
---|
[4830] | 480 | !--If vertical heterogeneity is taken into account, we use |
---|
| 481 | !--the "true" volume fraction instead of a modified |
---|
| 482 | !--surface fraction (which is larger and artificially |
---|
| 483 | !--reduces the in-cloud water). |
---|
[4803] | 484 | IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN |
---|
| 485 | eff_cldfra = ctot_vol(i) |
---|
| 486 | ELSE |
---|
| 487 | eff_cldfra = cldfra(i) |
---|
| 488 | ENDIF |
---|
| 489 | |
---|
[4809] | 490 | |
---|
[4830] | 491 | !--Start precipitation growth processes |
---|
[4809] | 492 | |
---|
| 493 | !--If the cloud is big enough, the precipitation processes activate |
---|
[4803] | 494 | IF ( cldfra(i) .GE. seuil_neb ) THEN |
---|
| 495 | |
---|
| 496 | !--------------------------------------------------------- |
---|
| 497 | !-- COLLECTION AND AGGREGATION |
---|
| 498 | !--------------------------------------------------------- |
---|
[4809] | 499 | !--Collection: processus through which rain collects small liquid droplets |
---|
| 500 | !--in suspension, and add it to the rain flux |
---|
| 501 | !--Aggregation: same for snow (precip flux) and ice crystals (in suspension) |
---|
| 502 | !--Those processes are treated before autoconversion because we do not |
---|
| 503 | !--want to collect/aggregate the newly formed fluxes, which already |
---|
| 504 | !--"saw" the cloud as they come from it |
---|
[4819] | 505 | !--The formulas come from Muench and Lohmann 2020 |
---|
| 506 | |
---|
[4809] | 507 | !--gamma_col: tuning coefficient [-] |
---|
| 508 | !--rho_rain: volumic mass of rain [kg/m3] |
---|
| 509 | !--r_rain: size of the rain droplets [m] |
---|
| 510 | !--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1) |
---|
| 511 | !--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are |
---|
| 512 | !--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to |
---|
| 513 | !--get in-cloud mean quantities. The two divisions by eff_cldfra are |
---|
| 514 | !--then simplified. |
---|
[4830] | 515 | |
---|
[4833] | 516 | !--The collection efficiency is perfect. |
---|
[4818] | 517 | Eff_rain_liq = 1. |
---|
[4809] | 518 | coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq |
---|
[4833] | 519 | IF ( raincld(i) .GT. 0. ) THEN |
---|
[4819] | 520 | !-- ATTENTION Double implicit version |
---|
[4830] | 521 | !-- BEWARE the formule below is FALSE (because raincld is a flux, not a delta flux) |
---|
[4833] | 522 | !qrain_tmp = raincld(i) / dhum_to_dflux(i) |
---|
[4819] | 523 | !coef_tmp = coef_col * dtime * ( qrain_tmp / precipfraccld(i) + qliq(i) / eff_cldfra ) |
---|
[4818] | 524 | !dqlcol = qliq(i) * ( 1. / ( 1. + 0.5 * ( coef_tmp - 1. + SQRT( & |
---|
[4819] | 525 | ! ( 1. - coef_tmp )**2. + 4. * coef_col * dtime * qrain_tmp / precipfraccld(i) ) & |
---|
[4818] | 526 | ! ) ) - 1. ) |
---|
| 527 | !--Barriers so that the processes do not consume more liquid/ice than |
---|
| 528 | !--available. |
---|
| 529 | !dqlcol = MAX( - qliq(i), dqlcol ) |
---|
[4830] | 530 | !--Exact explicit version, which does not need a barrier because of |
---|
[4819] | 531 | !--the exponential decrease |
---|
[4818] | 532 | dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. ) |
---|
| 533 | |
---|
| 534 | !--Add tendencies |
---|
| 535 | qliq(i) = qliq(i) + dqlcol |
---|
[4833] | 536 | raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i) |
---|
[4818] | 537 | |
---|
[4819] | 538 | !--Diagnostic tendencies |
---|
[4818] | 539 | dqrcol(i) = - dqlcol / dtime |
---|
[4809] | 540 | ENDIF |
---|
[4803] | 541 | |
---|
[4809] | 542 | !--Same as for aggregation |
---|
[4819] | 543 | !--Eff_snow_liq formula: following Milbrandt and Yau 2005, |
---|
| 544 | !--it s a product of a collection efficiency and a sticking efficiency |
---|
[4818] | 545 | Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) ) |
---|
| 546 | coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice |
---|
[4833] | 547 | IF ( snowcld(i) .GT. 0. ) THEN |
---|
[4819] | 548 | !-- ATTENTION Double implicit version? |
---|
[4818] | 549 | !--Barriers so that the processes do not consume more liquid/ice than |
---|
| 550 | !--available. |
---|
| 551 | !dqiagg = MAX( - qice(i), dqiagg ) |
---|
[4830] | 552 | !--Exact explicit version, which does not need a barrier because of |
---|
[4819] | 553 | !--the exponential decrease |
---|
[4818] | 554 | dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. ) |
---|
| 555 | |
---|
| 556 | !--Add tendencies |
---|
| 557 | qice(i) = qice(i) + dqiagg |
---|
[4833] | 558 | snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i) |
---|
[4818] | 559 | |
---|
[4819] | 560 | !--Diagnostic tendencies |
---|
[4818] | 561 | dqsagg(i) = - dqiagg / dtime |
---|
[4809] | 562 | ENDIF |
---|
[4803] | 563 | |
---|
| 564 | |
---|
| 565 | !--------------------------------------------------------- |
---|
| 566 | !-- AUTOCONVERSION |
---|
| 567 | !--------------------------------------------------------- |
---|
[4819] | 568 | !--Autoconversion converts liquid droplets/ice crystals into |
---|
| 569 | !--rain drops/snowflakes. It relies on the formulations by |
---|
| 570 | !--Sundqvist 1978. |
---|
[4803] | 571 | |
---|
[4819] | 572 | !--If we are in a convective point, we have different parameters |
---|
| 573 | !--for the autoconversion |
---|
| 574 | IF ( ptconv(i) ) THEN |
---|
[4803] | 575 | qthresh_auto_rain = cld_lc_con |
---|
[4830] | 576 | qthresh_auto_snow = cld_lc_con_snow |
---|
[4803] | 577 | |
---|
| 578 | tau_auto_rain = cld_tau_con |
---|
[4819] | 579 | !--tau for snow depends on the ice fraction in mixed-phase clouds |
---|
[4803] | 580 | tau_auto_snow = tau_auto_snow_max & |
---|
| 581 | + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) ) |
---|
| 582 | |
---|
| 583 | expo_auto_rain = cld_expo_con |
---|
| 584 | expo_auto_snow = cld_expo_con |
---|
| 585 | ELSE |
---|
| 586 | qthresh_auto_rain = cld_lc_lsc |
---|
[4830] | 587 | qthresh_auto_snow = cld_lc_lsc_snow |
---|
[4803] | 588 | |
---|
| 589 | tau_auto_rain = cld_tau_lsc |
---|
[4819] | 590 | !--tau for snow depends on the ice fraction in mixed-phase clouds |
---|
[4803] | 591 | tau_auto_snow = tau_auto_snow_max & |
---|
| 592 | + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) ) |
---|
| 593 | |
---|
| 594 | expo_auto_rain = cld_expo_lsc |
---|
| 595 | expo_auto_snow = cld_expo_lsc |
---|
| 596 | ENDIF |
---|
| 597 | |
---|
| 598 | |
---|
[4809] | 599 | ! Liquid water quantity to remove according to (Sundqvist, 1978) |
---|
[4830] | 600 | ! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) ) |
---|
| 601 | ! |
---|
| 602 | !--And same formula for ice |
---|
| 603 | ! |
---|
| 604 | !--We first treat the second term (with exponential) in an explicit way |
---|
| 605 | !--and then treat the first term (-q/tau) in an exact way |
---|
[4803] | 606 | |
---|
| 607 | dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( & |
---|
| 608 | - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) ) |
---|
| 609 | |
---|
| 610 | dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( & |
---|
| 611 | - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) ) |
---|
| 612 | |
---|
| 613 | |
---|
[4819] | 614 | !--Barriers so that we don t create more rain/snow |
---|
| 615 | !--than there is liquid/ice |
---|
[4803] | 616 | dqlauto = MAX( - qliq(i), dqlauto ) |
---|
| 617 | dqiauto = MAX( - qice(i), dqiauto ) |
---|
| 618 | |
---|
[4819] | 619 | !--Add tendencies |
---|
[4803] | 620 | qliq(i) = qliq(i) + dqlauto |
---|
| 621 | qice(i) = qice(i) + dqiauto |
---|
[4833] | 622 | raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i) |
---|
| 623 | snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i) |
---|
[4803] | 624 | |
---|
[4819] | 625 | !--Diagnostic tendencies |
---|
[4818] | 626 | dqsauto(i) = - dqiauto / dtime |
---|
| 627 | dqrauto(i) = - dqlauto / dtime |
---|
[4803] | 628 | |
---|
[4818] | 629 | |
---|
[4803] | 630 | !--------------------------------------------------------- |
---|
| 631 | !-- RIMING |
---|
| 632 | !--------------------------------------------------------- |
---|
[4819] | 633 | !--Process which converts liquid droplets in suspension into |
---|
| 634 | !--snow (graupel in fact) because of the collision between |
---|
| 635 | !--those and falling snowflakes. |
---|
[4830] | 636 | !--The formula comes from Muench and Lohmann 2020 |
---|
[4819] | 637 | !--NB.: this process needs a temperature adjustment |
---|
[4803] | 638 | |
---|
[4819] | 639 | !--Eff_snow_liq formula: following Seifert and Beheng 2006, |
---|
| 640 | !--assuming a cloud droplet diameter of 20 microns. |
---|
[4818] | 641 | Eff_snow_liq = 0.2 |
---|
| 642 | coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq |
---|
[4833] | 643 | IF ( snowcld(i) .GT. 0. ) THEN |
---|
[4819] | 644 | !-- ATTENTION Double implicit version? |
---|
[4818] | 645 | !--Barriers so that the processes do not consume more liquid than |
---|
| 646 | !--available. |
---|
| 647 | !dqlrim = MAX( - qliq(i), dqlrim ) |
---|
[4819] | 648 | !--Exact version, which does not need a barrier because of |
---|
| 649 | !--the exponential decrease |
---|
[4879] | 650 | dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. ) |
---|
[4809] | 651 | |
---|
[4819] | 652 | !--Add tendencies |
---|
[4833] | 653 | !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) |
---|
[4818] | 654 | qliq(i) = qliq(i) + dqlrim |
---|
[4833] | 655 | snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i) |
---|
[4818] | 656 | |
---|
[4819] | 657 | !--Temperature adjustment with the release of latent |
---|
| 658 | !--heat because of solid condensation |
---|
[4818] | 659 | temp(i) = temp(i) - dqlrim * RLMLT / RCPD & |
---|
| 660 | / ( 1. + RVTMP2 * qtot(i) ) |
---|
| 661 | |
---|
[4819] | 662 | !--Diagnostic tendencies |
---|
[4818] | 663 | dqsrim(i) = - dqlrim / dtime |
---|
[4809] | 664 | ENDIF |
---|
[4803] | 665 | |
---|
[4819] | 666 | ENDIF ! cldfra .GE. seuil_neb |
---|
[4803] | 667 | |
---|
[4819] | 668 | ENDDO ! loop on klon |
---|
[4803] | 669 | |
---|
[4819] | 670 | |
---|
| 671 | !--Re-calculation of saturation specific humidity |
---|
| 672 | !--because riming changed temperature |
---|
[4818] | 673 | CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat) |
---|
[4803] | 674 | |
---|
[4818] | 675 | DO i = 1, klon |
---|
[4803] | 676 | |
---|
| 677 | !--------------------------------------------------------- |
---|
[4818] | 678 | !-- MELTING |
---|
[4803] | 679 | !--------------------------------------------------------- |
---|
[4819] | 680 | !--Process through which snow melts into rain. |
---|
| 681 | !--The formula is homemade. |
---|
| 682 | !--NB.: this process needs a temperature adjustment |
---|
[4803] | 683 | |
---|
[4819] | 684 | !--dqsmelt_max: maximum snow melting so that temperature |
---|
| 685 | !-- stays higher than 273 K [kg/kg] |
---|
| 686 | !--capa_snowflake: capacitance of a snowflake, equal to |
---|
| 687 | !-- the radius if the snowflake is a sphere [m] |
---|
| 688 | !--temp_wetbulb: wet-bulb temperature [K] |
---|
[4830] | 689 | !--snow_fallspeed: snow fall velocity (in clear/cloudy sky) [m/s] |
---|
[4819] | 690 | !--air_thermal_conduct: thermal conductivity of the air [J/m/K/s] |
---|
| 691 | !--coef_ventil: ventilation coefficient [-] |
---|
| 692 | !--nb_snowflake: number of snowflakes (in clear/cloudy air) [-] |
---|
| 693 | |
---|
[4818] | 694 | IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN |
---|
[4819] | 695 | !--Computed according to |
---|
| 696 | !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q |
---|
[4818] | 697 | dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD & |
---|
| 698 | * ( 1. + RVTMP2 * qtot(i) )) |
---|
[4819] | 699 | |
---|
| 700 | !--Initialisation |
---|
[4818] | 701 | dqsclrmelt = 0. |
---|
| 702 | dqscldmelt = 0. |
---|
[4803] | 703 | |
---|
[4818] | 704 | capa_snowflake = r_snow |
---|
| 705 | ! ATTENTION POUR L'INSTANT C'EST UN WBULB SELON FORMULE ECMWF |
---|
| 706 | ! ATTENTION EST-CE QVAP ????? |
---|
| 707 | temp_wetbulb = temp(i) - ( qsat(i) - qvap(i) ) & |
---|
| 708 | * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) & |
---|
| 709 | - 40.637 * ( temp(i) - 275. ) ) |
---|
[4884] | 710 | air_thermal_conduct=(5.69+0.017*(temp(i)-RTT))*1.e-3*4.184 ! thermal conductivity of the air, SI |
---|
[4818] | 711 | |
---|
[4888] | 712 | |
---|
[4819] | 713 | !--In clear air |
---|
[4888] | 714 | IF ( ( snowclr(i) .GT. 0. ) .AND. ( precipfracclr(i) .GT. 0. ) ) THEN |
---|
[4819] | 715 | !--Calculated according to |
---|
| 716 | !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow |
---|
[4830] | 717 | nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr & |
---|
[4818] | 718 | / ( 4. / 3. * RPI * r_snow**3. * rho_snow ) |
---|
| 719 | dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct & |
---|
| 720 | * capa_snowflake / RLMLT * coef_ventil & |
---|
| 721 | * MAX(0., temp_wetbulb - RTT) * dtime |
---|
[4888] | 722 | |
---|
[4830] | 723 | !--Barrier to limit the melting flux to the clr snow flux in the mesh |
---|
[4833] | 724 | dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i)) |
---|
[4832] | 725 | ENDIF |
---|
[4818] | 726 | |
---|
[4888] | 727 | |
---|
[4819] | 728 | !--In cloudy air |
---|
[4888] | 729 | IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN |
---|
[4819] | 730 | !--Calculated according to |
---|
| 731 | !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow |
---|
[4830] | 732 | nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld & |
---|
[4818] | 733 | / ( 4. / 3. * RPI * r_snow**3. * rho_snow ) |
---|
| 734 | dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct & |
---|
| 735 | * capa_snowflake / RLMLT * coef_ventil & |
---|
| 736 | * MAX(0., temp_wetbulb - RTT) * dtime |
---|
[4830] | 737 | |
---|
| 738 | !--Barrier to limit the melting flux to the cld snow flux in the mesh |
---|
[4833] | 739 | dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i)) |
---|
[4888] | 740 | ENDIF |
---|
[4818] | 741 | |
---|
[4888] | 742 | |
---|
[4830] | 743 | !--Barrier on temperature. If the total melting flux leads to a |
---|
| 744 | !--positive temperature, it is limited to keep temperature above 0 degC. |
---|
[4819] | 745 | !--It is activated if the total is LOWER than the max |
---|
| 746 | !--because everything is negative |
---|
[4818] | 747 | dqstotmelt = dqsclrmelt + dqscldmelt |
---|
| 748 | IF ( dqstotmelt .LT. dqsmelt_max ) THEN |
---|
[4819] | 749 | !--We redistribute the max melted snow keeping |
---|
| 750 | !--the clear/cloud partition of the melted snow |
---|
[4818] | 751 | dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt |
---|
| 752 | dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt |
---|
| 753 | dqstotmelt = dqsmelt_max |
---|
[4832] | 754 | |
---|
[4818] | 755 | ENDIF |
---|
| 756 | |
---|
[4819] | 757 | !--Add tendencies |
---|
[4833] | 758 | !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) |
---|
| 759 | rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i)) |
---|
| 760 | raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i)) |
---|
| 761 | snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i)) |
---|
| 762 | snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i)) |
---|
[4818] | 763 | |
---|
[4819] | 764 | !--Temperature adjustment with the release of latent |
---|
| 765 | !--heat because of melting |
---|
[4818] | 766 | temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD & |
---|
| 767 | / ( 1. + RVTMP2 * qtot(i) ) |
---|
| 768 | |
---|
[4819] | 769 | !--Diagnostic tendencies |
---|
[4818] | 770 | dqrmelt(i) = - dqstotmelt / dtime |
---|
| 771 | dqsmelt(i) = dqstotmelt / dtime |
---|
| 772 | |
---|
| 773 | ENDIF |
---|
| 774 | |
---|
| 775 | |
---|
[4803] | 776 | !--------------------------------------------------------- |
---|
[4818] | 777 | !-- FREEZING |
---|
[4803] | 778 | !--------------------------------------------------------- |
---|
[4832] | 779 | !--Process through which rain freezes into snow. |
---|
| 780 | !-- We parameterize it as a 2 step process: |
---|
| 781 | !-- first: freezing following collision with ice crystals |
---|
| 782 | !-- second: immersion freezing following (inspired by Bigg 1953) |
---|
| 783 | !-- the latter is parameterized as an exponential decrease of the rain |
---|
| 784 | !-- water content with a homemade formulya |
---|
[4819] | 785 | !--This is based on a caracteritic time of freezing, which |
---|
| 786 | !--exponentially depends on temperature so that it is |
---|
| 787 | !--equal to 1 for temp_nowater (see below) and is close to |
---|
| 788 | !--0 for RTT (=273.15 K). |
---|
| 789 | !--NB.: this process needs a temperature adjustment |
---|
| 790 | !--dqrfreez_max: maximum rain freezing so that temperature |
---|
| 791 | !-- stays lower than 273 K [kg/kg] |
---|
| 792 | !--tau_freez: caracteristic time of freezing [s] |
---|
| 793 | !--gamma_freez: tuning parameter [s-1] |
---|
| 794 | !--alpha_freez: tuning parameter for the shape of the exponential curve [-] |
---|
| 795 | !--temp_nowater: temperature below which no liquid water exists [K] (about -40 degC) |
---|
| 796 | |
---|
[4818] | 797 | IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN |
---|
[4803] | 798 | |
---|
[4832] | 799 | |
---|
| 800 | !--1st step: freezing following collision with ice crystals |
---|
| 801 | !--Sub-process of freezing which quantifies the collision between |
---|
| 802 | !--ice crystals in suspension and falling rain droplets. |
---|
| 803 | !--The rain droplets freeze, becoming graupel, and carrying |
---|
| 804 | !--the ice crystal (which acted as an ice nucleating particle). |
---|
| 805 | !--The formula is adapted from the riming formula. |
---|
| 806 | !-- it works only in the cloudy part |
---|
| 807 | |
---|
| 808 | dqifreez = 0. |
---|
| 809 | dqrtotfreez_step1 = 0. |
---|
| 810 | |
---|
[4888] | 811 | IF ( ( qice(i) .GT. 0. ) .AND. ( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN |
---|
[4832] | 812 | dqrclrfreez = 0. |
---|
| 813 | dqrcldfreez = 0. |
---|
| 814 | |
---|
| 815 | !--Computed according to |
---|
| 816 | !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q |
---|
| 817 | dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD & |
---|
| 818 | * ( 1. + RVTMP2 * qtot(i) )) |
---|
| 819 | |
---|
| 820 | !--Sub-process of freezing which quantifies the collision between |
---|
| 821 | !--ice crystals in suspension and falling rain droplets. |
---|
| 822 | !--The rain droplets freeze, becoming graupel, and carrying |
---|
| 823 | !--the ice crystal (which acted as an ice nucleating particle). |
---|
| 824 | !--The formula is adapted from the riming formula. |
---|
| 825 | |
---|
| 826 | !--The collision efficiency is assumed unity |
---|
| 827 | Eff_rain_ice = 1. |
---|
| 828 | coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice |
---|
| 829 | !--Exact version, which does not need a barrier because of |
---|
| 830 | !--the exponential decrease. |
---|
| 831 | dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. ) |
---|
| 832 | |
---|
| 833 | !--We add the part of rain water that freezes, limited by a temperature barrier |
---|
[4833] | 834 | !--this quantity is calculated assuming that the number of drop that freeze correspond to the number |
---|
| 835 | !--of crystals collected (and assuming uniform distributions of ice crystals and rain drops) |
---|
| 836 | !--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3 |
---|
| 837 | !--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3 |
---|
| 838 | !--The assumption above corresponds to dNi = dNr, i.e., |
---|
| 839 | !-- dqr = dqi * (4/3 PI rho_rain * r_rain**3) / (4/3 PI rho_ice * r_ice**3) |
---|
| 840 | dqrcldfreez = dqifreez * rho_rain * r_rain**3. / ( rho_ice * r_ice**3. ) |
---|
| 841 | dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i)) |
---|
[4832] | 842 | dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max) |
---|
| 843 | dqrtotfreez_step1 = dqrcldfreez |
---|
| 844 | |
---|
| 845 | !--Add tendencies |
---|
[4833] | 846 | !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) |
---|
[4832] | 847 | qice(i) = qice(i) + dqifreez |
---|
[4833] | 848 | raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i)) |
---|
| 849 | snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i)) |
---|
[4832] | 850 | temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD & |
---|
[4833] | 851 | / ( 1. + RVTMP2 * qtot(i) ) |
---|
[4832] | 852 | |
---|
| 853 | ENDIF |
---|
| 854 | |
---|
| 855 | !-- Second step immersion freezing of rain drops |
---|
| 856 | !-- with a homemade timeconstant depending on temperature |
---|
| 857 | |
---|
| 858 | dqrclrfreez = 0. |
---|
| 859 | dqrcldfreez = 0. |
---|
| 860 | dqrtotfreez_step2 = 0. |
---|
[4819] | 861 | !--Computed according to |
---|
| 862 | !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q |
---|
[4832] | 863 | |
---|
[4818] | 864 | dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD & |
---|
| 865 | * ( 1. + RVTMP2 * qtot(i) )) |
---|
[4832] | 866 | |
---|
| 867 | |
---|
[4830] | 868 | tau_freez = 1. / ( beta_freez & |
---|
[4818] | 869 | * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) ) |
---|
[4803] | 870 | |
---|
[4830] | 871 | |
---|
[4823] | 872 | !--In clear air |
---|
| 873 | IF ( rainclr(i) .GT. 0. ) THEN |
---|
[4830] | 874 | !--Exact solution of dqrain/dt = -qrain/tau_freez |
---|
[4833] | 875 | dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. ) |
---|
[4823] | 876 | ENDIF |
---|
[4803] | 877 | |
---|
[4823] | 878 | !--In cloudy air |
---|
| 879 | IF ( raincld(i) .GT. 0. ) THEN |
---|
[4830] | 880 | !--Exact solution of dqrain/dt = -qrain/tau_freez |
---|
[4833] | 881 | dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. ) |
---|
[4823] | 882 | ENDIF |
---|
[4803] | 883 | |
---|
[4832] | 884 | !--temperature barrie step 2 |
---|
[4823] | 885 | !--It is activated if the total is LOWER than the max |
---|
| 886 | !--because everything is negative |
---|
[4832] | 887 | dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez |
---|
| 888 | |
---|
| 889 | IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN |
---|
[4823] | 890 | !--We redistribute the max freezed rain keeping |
---|
| 891 | !--the clear/cloud partition of the freezing rain |
---|
[4832] | 892 | dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2 |
---|
| 893 | dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2 |
---|
| 894 | dqrtotfreez_step2 = dqrfreez_max |
---|
[4823] | 895 | ENDIF |
---|
| 896 | |
---|
| 897 | |
---|
[4819] | 898 | !--Add tendencies |
---|
[4833] | 899 | !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) |
---|
| 900 | rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i)) |
---|
| 901 | raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i)) |
---|
| 902 | snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i)) |
---|
| 903 | snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i)) |
---|
[4803] | 904 | |
---|
[4832] | 905 | |
---|
[4819] | 906 | !--Temperature adjustment with the uptake of latent |
---|
| 907 | !--heat because of freezing |
---|
[4832] | 908 | temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD & |
---|
[4818] | 909 | / ( 1. + RVTMP2 * qtot(i) ) |
---|
| 910 | |
---|
[4819] | 911 | !--Diagnostic tendencies |
---|
[4833] | 912 | dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2 |
---|
[4818] | 913 | dqrfreez(i) = dqrtotfreez / dtime |
---|
[4832] | 914 | dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime |
---|
[4818] | 915 | |
---|
| 916 | ENDIF |
---|
| 917 | |
---|
[4888] | 918 | |
---|
| 919 | |
---|
[4819] | 920 | !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min, |
---|
| 921 | !--we reduce the precipiration fraction in the clear/cloudy air so that the new |
---|
| 922 | !--local flux of rain+snow is equal to rain_int_min. |
---|
| 923 | !--Here, rain+snow is the gridbox-mean flux of precip. |
---|
| 924 | !--Therefore, (rain+snow)/precipfrac is the local flux of precip. |
---|
| 925 | !--If the local flux of precip is lower than rain_int_min, i.e., |
---|
| 926 | !-- (rain+snow)/precipfrac < rain_int_min , i.e., |
---|
| 927 | !-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce |
---|
| 928 | !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min. |
---|
| 929 | !--Note that this is physically different than what is proposed in LTP thesis. |
---|
[4803] | 930 | precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min ) |
---|
| 931 | precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min ) |
---|
| 932 | |
---|
[4819] | 933 | !--Calculate outputs |
---|
[4803] | 934 | rain(i) = rainclr(i) + raincld(i) |
---|
| 935 | snow(i) = snowclr(i) + snowcld(i) |
---|
| 936 | |
---|
[4819] | 937 | !--Diagnostics |
---|
[4830] | 938 | !--BEWARE this is indeed a diagnostic: this is an estimation from |
---|
| 939 | !--the value of the flux at the bottom interface of the mesh and |
---|
| 940 | !--and assuming an upstream numerical calculation |
---|
| 941 | !--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is |
---|
| 942 | !--used for computing the total ice water content in the mesh |
---|
| 943 | !--for radiation only |
---|
| 944 | qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr & |
---|
| 945 | + raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld ) |
---|
| 946 | qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr & |
---|
| 947 | + snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld ) |
---|
[4803] | 948 | |
---|
[4888] | 949 | |
---|
[4819] | 950 | ENDDO ! loop on klon |
---|
[4803] | 951 | |
---|
[4888] | 952 | |
---|
[4803] | 953 | END SUBROUTINE poprecip_postcld |
---|
| 954 | |
---|
| 955 | END MODULE lmdz_lscp_poprecip |
---|