source: LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90 @ 4913

Last change on this file since 4913 was 4913, checked in by evignon, 4 weeks ago

amelioration de poprecip avec:

  • formule du rayon des cristaux
  • qvap du ciel clair pour l'évap des precip

Lea et Audran
File size: 42.2 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  ELSE
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: following Milbrandt and Yau 2005,
562    !--it s a product of a collection efficiency and a sticking efficiency
563    Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
564    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
565    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
566    coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
567    IF ( snowcld(i) .GT. 0. ) THEN
568      !--Exact explicit version, which does not need a barrier because of
569      !--the exponential decrease
570      dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
571
572      !--Add tendencies
573      qice(i) = qice(i) + dqiagg
574      snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i)
575
576      !--Diagnostic tendencies
577      dqsagg(i)  = - dqiagg  / dtime
578    ENDIF
579
580
581    !---------------------------------------------------------
582    !--                  AUTOCONVERSION
583    !---------------------------------------------------------
584    !--Autoconversion converts liquid droplets/ice crystals into
585    !--rain drops/snowflakes. It relies on the formulations by
586    !--Sundqvist 1978.
587
588    !--If we are in a convective point, we have different parameters
589    !--for the autoconversion
590    IF ( ptconv(i) ) THEN
591        qthresh_auto_rain = cld_lc_con
592        qthresh_auto_snow = cld_lc_con_snow
593
594        tau_auto_rain = cld_tau_con
595        !--tau for snow depends on the ice fraction in mixed-phase clouds     
596        tau_auto_snow = tau_auto_snow_max &
597                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
598
599        expo_auto_rain = cld_expo_con
600        expo_auto_snow = cld_expo_con
601    ELSE
602        qthresh_auto_rain = cld_lc_lsc
603        qthresh_auto_snow = cld_lc_lsc_snow
604
605        tau_auto_rain = cld_tau_lsc
606        !--tau for snow depends on the ice fraction in mixed-phase clouds     
607        tau_auto_snow = tau_auto_snow_max &
608                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
609
610        expo_auto_rain = cld_expo_lsc
611        expo_auto_snow = cld_expo_lsc
612    ENDIF
613
614
615    ! Liquid water quantity to remove according to (Sundqvist, 1978)
616    ! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) )
617    !
618    !--And same formula for ice
619    !
620    !--We first treat the second term (with exponential) in an explicit way
621    !--and then treat the first term (-q/tau) in an exact way
622
623    dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
624               - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
625
626    dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
627               - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
628
629
630    !--Barriers so that we don't create more rain/snow
631    !--than there is liquid/ice
632    dqlauto = MAX( - qliq(i), dqlauto )
633    dqiauto = MAX( - qice(i), dqiauto )
634
635    !--Add tendencies
636    qliq(i) = qliq(i) + dqlauto
637    qice(i) = qice(i) + dqiauto
638    raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
639    snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
640
641    !--Diagnostic tendencies
642    dqsauto(i) = - dqiauto / dtime
643    dqrauto(i) = - dqlauto / dtime
644
645
646    !---------------------------------------------------------
647    !--                  RIMING
648    !---------------------------------------------------------
649    !--Process which converts liquid droplets in suspension into
650    !--snow because of the collision between
651    !--those and falling snowflakes.
652    !--The formula comes from Muench and Lohmann 2020
653    !--NB.: this process needs a temperature adjustment
654
655    !--Eff_snow_liq formula: following Seifert and Beheng 2006,
656    !--assuming a cloud droplet diameter of 20 microns.
657    Eff_snow_liq = 0.2
658    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
659    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
660    coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
661    IF ( snowcld(i) .GT. 0. ) THEN
662      !--Exact version, which does not need a barrier because of
663      !--the exponential decrease
664      dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. )
665
666      !--Add tendencies
667      !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
668      qliq(i) = qliq(i) + dqlrim
669      snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i)
670
671      !--Temperature adjustment with the release of latent
672      !--heat because of solid condensation
673      temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
674                        / ( 1. + RVTMP2 * qtot(i) )
675
676      !--Diagnostic tendencies
677      dqsrim(i)  = - dqlrim  / dtime
678    ENDIF
679
680  ENDIF ! cldfra .GE. seuil_neb
681
682ENDDO ! loop on klon
683
684
685!--Re-calculation of saturation specific humidity
686!--because riming changed temperature
687CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
688
689DO i = 1, klon
690
691  !---------------------------------------------------------
692  !--                  MELTING
693  !---------------------------------------------------------
694  !--Process through which snow melts into rain.
695  !--The formula is homemade.
696  !--NB.: this process needs a temperature adjustment
697
698  !--dqsmelt_max         : maximum snow melting so that temperature
699  !--                      stays higher than 273 K [kg/kg]
700  !--capa_snowflake      : capacitance of a snowflake, equal to
701  !--                      the radius if the snowflake is a sphere [m]
702  !--temp_wetbulb        : wet-bulb temperature [K]
703  !--snow_fallspeed      : snow fall velocity (in clear/cloudy sky) [m/s]
704  !--air_thermal_conduct : thermal conductivity of the air [J/m/K/s]
705  !--gamma_melt          : tuning parameter for melting [-]
706  !--nb_snowflake        : number of snowflakes (in clear/cloudy air) [-]
707
708  IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
709    !--Computed according to
710    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
711    dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
712                        * ( 1. + RVTMP2 * qtot(i) ))
713   
714    !--Initialisation
715    dqsclrmelt = 0.
716    dqscldmelt = 0.
717
718    !--We assume that the snowflakes are spherical
719    capa_snowflake = r_snow
720    !--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
721    air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184   
722    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
723    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
724
725    !--In clear air
726    IF ( ( snowclr(i) .GT. 0. ) .AND.  ( precipfracclr(i) .GT. 0. ) ) THEN
727      !--Formula for the wet-bulb temperature from ECMWF (IFS)
728      !--The vapor used is the vapor in the clear sky
729      temp_wetbulb = temp(i) &
730                   - ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) &
731                   * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
732                   - 40.637 * ( temp(i) - 275. ) )
733      !--Calculated according to
734      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
735      nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr &
736                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
737      dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
738                   * capa_snowflake / RLMLT * gamma_melt &
739                   * MAX(0., temp_wetbulb - RTT) * dtime
740     
741      !--Barrier to limit the melting flux to the clr snow flux in the mesh
742      dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i))
743    ENDIF
744
745
746    !--In cloudy air
747    IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
748      !--As the air is saturated, the wet-bulb temperature is equal to the
749      !--temperature
750      temp_wetbulb = temp(i)
751      !--Calculated according to
752      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
753      nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld &
754                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
755      dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
756                   * capa_snowflake / RLMLT * gamma_melt &
757                   * MAX(0., temp_wetbulb - RTT) * dtime
758
759      !--Barrier to limit the melting flux to the cld snow flux in the mesh
760      dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i))
761   ENDIF
762   
763
764    !--Barrier on temperature. If the total melting flux leads to a
765    !--positive temperature, it is limited to keep temperature above 0 degC.
766    !--It is activated if the total is LOWER than the max
767    !--because everything is negative
768    dqstotmelt = dqsclrmelt + dqscldmelt
769    IF ( dqstotmelt .LT. dqsmelt_max ) THEN
770      !--We redistribute the max melted snow keeping
771      !--the clear/cloud partition of the melted snow
772      dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
773      dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
774      dqstotmelt = dqsmelt_max
775
776    ENDIF
777
778    !--Add tendencies
779    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
780    rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i))
781    raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i))
782    snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i))
783    snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i))
784
785    !--Temperature adjustment with the release of latent
786    !--heat because of melting
787    temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
788                      / ( 1. + RVTMP2 * qtot(i) )
789
790    !--Diagnostic tendencies
791    dqrmelt(i) = - dqstotmelt / dtime
792    dqsmelt(i) = dqstotmelt / dtime
793
794  ENDIF
795
796
797  !---------------------------------------------------------
798  !--                  FREEZING
799  !---------------------------------------------------------
800  !--Process through which rain freezes into snow.
801  !-- We parameterize it as a 2 step process:
802  !--first: freezing following collision with ice crystals
803  !--second: immersion freezing following (inspired by Bigg 1953)
804  !--the latter is parameterized as an exponential decrease of the rain
805  !--water content with a homemade formulya
806  !--This is based on a caracteritic time of freezing, which
807  !--exponentially depends on temperature so that it is
808  !--equal to 1 for temp_nowater (see below) and is close to
809  !--0 for RTT (=273.15 K).
810  !--NB.: this process needs a temperature adjustment
811  !--dqrfreez_max : maximum rain freezing so that temperature
812  !--              stays lower than 273 K [kg/kg]
813  !--tau_freez    : caracteristic time of freezing [s]
814  !--gamma_freez  : tuning parameter [s-1]
815  !--alpha_freez  : tuning parameter for the shape of the exponential curve [-]
816  !--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC)
817
818  IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
819
820 
821    !--1st step: freezing following collision with ice crystals
822    !--Sub-process of freezing which quantifies the collision between
823    !--ice crystals in suspension and falling rain droplets.
824    !--The rain droplets freeze, becoming graupel, and carrying
825    !--the ice crystal (which acted as an ice nucleating particle).
826    !--The formula is adapted from the riming formula.
827    !--it works only in the cloudy part
828
829    dqifreez = 0.
830    dqrtotfreez_step1 = 0.
831
832    IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. &
833         ( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
834      dqrclrfreez = 0.
835      dqrcldfreez = 0.
836
837      !--Computed according to
838      !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
839      dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
840                         * ( 1. + RVTMP2 * qtot(i) ))
841 
842
843      !--The collision efficiency is assumed unity
844      Eff_rain_ice = 1.
845      coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
846      !--Exact version, which does not need a barrier because of
847      !--the exponential decrease.
848      dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
849
850      !--We add the part of rain water that freezes, limited by a temperature barrier
851      !--This quantity is calculated assuming that the number of drop that freeze correspond to the number
852      !--of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
853      !--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3
854      !--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3
855      !--The assumption above corresponds to dNi = dNr, i.e.,
856      !-- dqr = dqi * (4/3 PI rho_rain * r_rain**3) / (4/3 PI rho_ice * r_ice**3)
857      !--Dry density [kg/m3]
858      rho = pplay(i) / temp(i) / RD
859      !--r_ice formula from Sun and Rikus (1999)
860      r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 &
861            + 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
862      dqrcldfreez = dqifreez * rho_rain * r_rain**3. / ( rho_ice * r_ice**3. )
863      dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i))
864      dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max)
865      dqrtotfreez_step1 = dqrcldfreez
866
867      !--Add tendencies
868      !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
869      qice(i) = qice(i) + dqifreez
870      raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
871      snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i))
872      temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD &
873                        / ( 1. + RVTMP2 * qtot(i) )
874
875    ENDIF
876   
877    !-- Second step immersion freezing of rain drops
878    !-- with a homemade timeconstant depending on temperature
879
880    dqrclrfreez = 0.
881    dqrcldfreez = 0.
882    dqrtotfreez_step2 = 0.
883    !--Computed according to
884    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
885
886    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
887                         * ( 1. + RVTMP2 * qtot(i) ))
888
889   
890    tau_freez = 1. / ( beta_freez &
891              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
892
893
894    !--In clear air
895    IF ( rainclr(i) .GT. 0. ) THEN
896      !--Exact solution of dqrain/dt = -qrain/tau_freez
897      dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
898    ENDIF
899
900    !--In cloudy air
901    IF ( raincld(i) .GT. 0. ) THEN
902      !--Exact solution of dqrain/dt = -qrain/tau_freez
903      dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
904    ENDIF
905
906    !--temperature barrier step 2
907    !--It is activated if the total is LOWER than the max
908    !--because everything is negative
909    dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
910 
911    IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
912      !--We redistribute the max freezed rain keeping
913      !--the clear/cloud partition of the freezing rain
914      dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
915      dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
916      dqrtotfreez_step2 = dqrfreez_max
917    ENDIF
918
919
920    !--Add tendencies
921    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
922    rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i))
923    raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
924    snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i))
925    snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i))
926
927
928    !--Temperature adjustment with the uptake of latent
929    !--heat because of freezing
930    temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
931                      / ( 1. + RVTMP2 * qtot(i) )
932
933    !--Diagnostic tendencies
934    dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2         
935    dqrfreez(i) = dqrtotfreez / dtime
936    dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
937
938  ENDIF
939
940
941
942  !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
943  !--we reduce the precipiration fraction in the clear/cloudy air so that the new
944  !--local flux of rain+snow is equal to rain_int_min.
945  !--Here, rain+snow is the gridbox-mean flux of precip.
946  !--Therefore, (rain+snow)/precipfrac is the local flux of precip.
947  !--If the local flux of precip is lower than rain_int_min, i.e.,
948  !-- (rain+snow)/precipfrac < rain_int_min , i.e.,
949  !-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
950  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
951  !--Note that this is physically different than what is proposed in LTP thesis.
952  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
953  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
954
955  !--Calculate outputs
956  rain(i) = rainclr(i) + raincld(i)
957  snow(i) = snowclr(i) + snowcld(i)
958
959  !--Diagnostics
960  !--BEWARE this is indeed a diagnostic: this is an estimation from
961  !--the value of the flux at the bottom interface of the mesh and
962  !--and assuming an upstream numerical calculation
963  !--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is
964  !--used for computing the total ice water content in the mesh
965  !--for radiation only
966  qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr &
967                 + raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld )
968  qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr &
969                 + snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld )
970
971
972ENDDO ! loop on klon
973
974
975END SUBROUTINE poprecip_postcld
976
977END MODULE lmdz_lscp_poprecip
Note: See TracBrowser for help on using the repository browser.