source: LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.f90 @ 5296

Last change on this file since 5296 was 5268, checked in by abarral, 8 days ago

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