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

Last change on this file since 5006 was 4960, checked in by evignon, 7 months ago

petite commission dans poprecip suite à l'atelier nuage du jour
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  ELSEIF ( ( rain(i) + snow(i) ) .LE. 0. ) THEN
249
250    !--If no precip, we reinitialize the cloud fraction used for the precip to 0
251    precipfraccld(i) = 0.
252    precipfracclr(i) = 0.
253
254  ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. )
255
256  !--Diagnostic tendencies
257  dqssub(i) = dqssubl / dtime
258  dqreva(i) = dqrevap / dtime
259
260ENDDO ! loop on klon
261
262
263END SUBROUTINE poprecip_precld
264
265!----------------------------------------------------------------
266! Computes the processes-oriented precipitation formulations for
267! - autoconversion (auto) via a deposition process
268! - aggregation (agg)
269! - riming (rim)
270! - collection (col)
271! - melting (melt)
272! - freezing (freez)
273!
274SUBROUTINE poprecip_postcld( &
275           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
276           temp, qvap, qliq, qice, icefrac, cldfra, &
277           precipfracclr, precipfraccld, &
278           rain, rainclr, raincld, snow, snowclr, snowcld, &
279           qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
280           dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
281
282USE lmdz_lscp_ini, ONLY : prt_level, lunout
283USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
284USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
285
286USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb,    &
287                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, &
288                          thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
289                          rho_rain, r_rain, r_snow, rho_ice,                   &
290                          tau_auto_snow_min, tau_auto_snow_max,                &
291                          thresh_precip_frac, eps,                             &
292                          gamma_melt, alpha_freez, beta_freez, temp_nowater,   &
293                          iflag_cloudth_vert, iflag_rain_incloud_vol,          &
294                          cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez,       &
295                          rain_fallspeed_clr, rain_fallspeed_cld,              &
296                          snow_fallspeed_clr, snow_fallspeed_cld
297
298
299IMPLICIT NONE
300
301INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
302REAL,    INTENT(IN)                     :: dtime          !--time step [s]
303
304REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
305REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
306REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
307
308REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
309LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
310
311REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
312REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity [kg/kg]
313REAL,    INTENT(INOUT), DIMENSION(klon) :: qliq           !--current liquid water specific humidity [kg/kg]
314REAL,    INTENT(INOUT), DIMENSION(klon) :: qice           !--current ice water specific humidity [kg/kg]
315REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--ice fraction [-]
316REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction [-]
317
318REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
319REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
320                                                          !--NB. at the end of the routine, becomes the fraction of precip
321                                                          !--in the current layer
322
323REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
324REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
325REAL,    INTENT(INOUT), DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
326REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
327REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
328REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
329
330REAL,    INTENT(OUT),   DIMENSION(klon) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
331REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
332REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
333REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
334REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
335REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
336REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsrim         !--snow tendency due to riming [kg/kg/s]
337REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
338REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
339REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
340REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
341
342
343
344!--Local variables
345
346INTEGER :: i
347!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
348REAL, DIMENSION(klon) :: dhum_to_dflux
349REAL, DIMENSION(klon) :: qtot                             !--includes vap, liq, ice and precip
350
351!--Partition of the fluxes
352REAL :: dcldfra
353REAL :: precipfractot
354REAL :: dprecipfracclr, dprecipfraccld
355REAL :: drainclr, dsnowclr
356REAL :: draincld, dsnowcld
357
358!--Collection, aggregation and riming
359REAL :: eff_cldfra
360REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
361REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
362REAL :: rho_snow
363REAL :: dqlcol           !--loss of liquid cloud content due to collection by rain [kg/kg/s]
364REAL :: dqiagg           !--loss of ice cloud content due to collection by aggregation [kg/kg/s]
365REAL :: dqlrim           !--loss of liquid cloud content due to riming on snow [kg/kg/s]
366
367!--Autoconversion
368REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
369REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
370REAL :: dqlauto          !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
371REAL :: dqiauto          !--loss of ice cloud content due to autoconversion to snow [kg/kg/s]
372
373!--Melting
374REAL :: dqsmelt_max, air_thermal_conduct
375REAL :: nb_snowflake_clr, nb_snowflake_cld
376REAL :: capa_snowflake, temp_wetbulb
377REAL :: rho, r_ice
378REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
379REAL, DIMENSION(klon) :: qzero, qsat, dqsat
380
381!--Freezing
382REAL :: dqrfreez_max
383REAL :: tau_freez
384REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2
385REAL :: coef_freez
386REAL :: dqifreez        !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
387REAL :: Eff_rain_ice
388
389
390!--Initialisation of variables
391
392
393qzero(:)    = 0.
394
395dqrcol(:)   = 0.
396dqsagg(:)   = 0.
397dqsauto(:)  = 0.
398dqrauto(:)  = 0.
399dqsrim(:)   = 0.
400dqrmelt(:)  = 0.
401dqsmelt(:)  = 0.
402dqrfreez(:) = 0.
403dqsfreez(:) = 0.
404
405
406DO i = 1, klon
407
408  !--Variables initialisation
409  dqlcol  = 0.
410  dqiagg  = 0.
411  dqiauto = 0.
412  dqlauto = 0.
413  dqlrim  = 0.
414
415  !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
416  dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
417  qtot(i) = qvap(i) + qliq(i) + qice(i) &
418          + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i)
419
420  !------------------------------------------------------------
421  !--     PRECIPITATION FRACTIONS UPDATE
422  !------------------------------------------------------------
423  !--The goal of this routine is to reattribute precipitation fractions
424  !--and fluxes to clear or cloudy air, depending on the variation of
425  !--the cloud fraction on the vertical dimension. We assume a
426  !--maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
427  !--and LTP thesis, 2021)
428  !--NB. in fact, we assume a maximum-random overlap of the total precip. frac
429
430  !--Initialisation
431  precipfractot = precipfracclr(i) + precipfraccld(i)
432
433  !--Instead of using the cloud cover which was use in LTP thesis, we use the
434  !--total precip. fraction to compute the maximum-random overlap. This is
435  !--because all the information of the cloud cover is embedded into
436  !--precipfractot, and this allows for taking into account the potential
437  !--reduction of the precipitation fraction because either the flux is too
438  !--small (see barrier at the end of poprecip_postcld) or the flux is completely
439  !--evaporated (see barrier at the end of poprecip_precld)
440  !--NB. precipfraccld(i) is here the cloud fraction of the layer above
441  !precipfractot = 1. - ( 1. - precipfractot ) * &
442  !               ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
443  !             / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
444
445
446  IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN
447    precipfractot = 1.
448  ELSE
449    precipfractot = 1. - ( 1. - precipfractot ) * &
450                 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
451               / ( 1. - precipfraccld(i) )
452  ENDIF
453
454  !--precipfraccld(i) is here the cloud fraction of the layer above
455  dcldfra = cldfra(i) - precipfraccld(i)
456  !--Tendency of the clear-sky precipitation fraction. We add a MAX on the
457  !--calculation of the current CS precip. frac.
458  !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
459  !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated
460  !--if precipfractot < cldfra)
461  dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i)
462  !--Tendency of the cloudy precipitation fraction. We add a MAX on the
463  !--calculation of the current CS precip. frac.
464  !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) )
465  !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated
466  !--if cldfra < 0)
467  dprecipfraccld = dcldfra
468
469
470  !--If the cloud extends
471  IF ( dprecipfraccld .GT. 0. ) THEN
472    !--If there is no CS precip, nothing happens.
473    !--If there is, we reattribute some of the CS precip flux
474    !--to the cloud precip flux, proportionnally to the
475    !--decrease of the CS precip fraction
476    IF ( precipfracclr(i) .LE. 0. ) THEN
477      drainclr = 0.
478      dsnowclr = 0.
479    ELSE
480      drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i)
481      dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i)
482    ENDIF
483  !--If the cloud narrows
484  ELSEIF ( dprecipfraccld .LT. 0. ) THEN
485    !--We reattribute some of the cloudy precip flux
486    !--to the CS precip flux, proportionnally to the
487    !--decrease of the cloud precip fraction
488    draincld = dprecipfraccld / precipfraccld(i) * raincld(i)
489    dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i)
490    drainclr = - draincld
491    dsnowclr = - dsnowcld
492  !--If the cloud stays the same or if there is no cloud above and
493  !--in the current layer, nothing happens
494  ELSE
495    drainclr = 0.
496    dsnowclr = 0.
497  ENDIF
498
499  !--We add the tendencies
500  !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
501  precipfraccld(i) = precipfraccld(i) + dprecipfraccld
502  precipfracclr(i) = precipfracclr(i) + dprecipfracclr
503  rainclr(i) = MAX(0., rainclr(i) + drainclr)
504  snowclr(i) = MAX(0., snowclr(i) + dsnowclr)
505  raincld(i) = MAX(0., raincld(i) - drainclr)
506  snowcld(i) = MAX(0., snowcld(i) - dsnowclr)
507 
508  !--If vertical heterogeneity is taken into account, we use
509  !--the "true" volume fraction instead of a modified
510  !--surface fraction (which is larger and artificially
511  !--reduces the in-cloud water).
512  IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
513    eff_cldfra = ctot_vol(i)
514  ELSE
515    eff_cldfra = cldfra(i)
516  ENDIF
517
518
519  !--Start precipitation growth processes
520
521  !--If the cloud is big enough, the precipitation processes activate
522  ! TODO met on seuil_neb ici ?
523  IF ( cldfra(i) .GE. seuil_neb ) THEN
524
525    !---------------------------------------------------------
526    !--           COLLECTION AND AGGREGATION
527    !---------------------------------------------------------
528    !--Collection: processus through which rain collects small liquid droplets
529    !--in suspension, and add it to the rain flux
530    !--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
531    !--Those processes are treated before autoconversion because we do not
532    !--want to collect/aggregate the newly formed fluxes, which already
533    !--"saw" the cloud as they come from it
534    !--The formulas come from Muench and Lohmann 2020
535
536    !--gamma_col: tuning coefficient [-]
537    !--rho_rain: volumic mass of rain [kg/m3]
538    !--r_rain: size of the rain droplets [m]
539    !--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
540    !--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
541    !--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
542    !--get in-cloud mean quantities. The two divisions by eff_cldfra are
543    !--then simplified.
544
545    !--The collection efficiency is perfect.
546    Eff_rain_liq = 1.
547    coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
548    IF ( raincld(i) .GT. 0. ) THEN
549      !--Exact explicit version, which does not need a barrier because of
550      !--the exponential decrease
551      dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
552
553      !--Add tendencies
554      qliq(i) = qliq(i) + dqlcol
555      raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i)
556
557      !--Diagnostic tendencies
558      dqrcol(i)  = - dqlcol  / dtime
559    ENDIF
560
561    !--Same as for aggregation
562    !--Eff_snow_liq formula: following Milbrandt and Yau 2005,
563    !--it s a product of a collection efficiency and a sticking efficiency
564    Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
565    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
566    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
567    coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
568    IF ( snowcld(i) .GT. 0. ) THEN
569      !--Exact explicit version, which does not need a barrier because of
570      !--the exponential decrease
571      dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
572
573      !--Add tendencies
574      qice(i) = qice(i) + dqiagg
575      snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i)
576
577      !--Diagnostic tendencies
578      dqsagg(i)  = - dqiagg  / dtime
579    ENDIF
580
581
582    !---------------------------------------------------------
583    !--                  AUTOCONVERSION
584    !---------------------------------------------------------
585    !--Autoconversion converts liquid droplets/ice crystals into
586    !--rain drops/snowflakes. It relies on the formulations by
587    !--Sundqvist 1978.
588
589    !--If we are in a convective point, we have different parameters
590    !--for the autoconversion
591    IF ( ptconv(i) ) THEN
592        qthresh_auto_rain = cld_lc_con
593        qthresh_auto_snow = cld_lc_con_snow
594
595        tau_auto_rain = cld_tau_con
596        !--tau for snow depends on the ice fraction in mixed-phase clouds     
597        tau_auto_snow = tau_auto_snow_max &
598                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
599
600        expo_auto_rain = cld_expo_con
601        expo_auto_snow = cld_expo_con
602    ELSE
603        qthresh_auto_rain = cld_lc_lsc
604        qthresh_auto_snow = cld_lc_lsc_snow
605
606        tau_auto_rain = cld_tau_lsc
607        !--tau for snow depends on the ice fraction in mixed-phase clouds     
608        tau_auto_snow = tau_auto_snow_max &
609                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
610
611        expo_auto_rain = cld_expo_lsc
612        expo_auto_snow = cld_expo_lsc
613    ENDIF
614
615
616    ! Liquid water quantity to remove according to (Sundqvist, 1978)
617    ! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) )
618    !
619    !--And same formula for ice
620    !
621    !--We first treat the second term (with exponential) in an explicit way
622    !--and then treat the first term (-q/tau) in an exact way
623
624    dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
625               - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
626
627    dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
628               - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
629
630
631    !--Barriers so that we don't create more rain/snow
632    !--than there is liquid/ice
633    dqlauto = MAX( - qliq(i), dqlauto )
634    dqiauto = MAX( - qice(i), dqiauto )
635
636    !--Add tendencies
637    qliq(i) = qliq(i) + dqlauto
638    qice(i) = qice(i) + dqiauto
639    raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
640    snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
641
642    !--Diagnostic tendencies
643    dqsauto(i) = - dqiauto / dtime
644    dqrauto(i) = - dqlauto / dtime
645
646
647    !---------------------------------------------------------
648    !--                  RIMING
649    !---------------------------------------------------------
650    !--Process which converts liquid droplets in suspension into
651    !--snow because of the collision between
652    !--those and falling snowflakes.
653    !--The formula comes from Muench and Lohmann 2020
654    !--NB.: this process needs a temperature adjustment
655
656    !--Eff_snow_liq formula: following Seifert and Beheng 2006,
657    !--assuming a cloud droplet diameter of 20 microns.
658    Eff_snow_liq = 0.2
659    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
660    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
661    coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
662    IF ( snowcld(i) .GT. 0. ) THEN
663      !--Exact version, which does not need a barrier because of
664      !--the exponential decrease
665      dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. )
666
667      !--Add tendencies
668      !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
669      qliq(i) = qliq(i) + dqlrim
670      snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i)
671
672      !--Temperature adjustment with the release of latent
673      !--heat because of solid condensation
674      temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
675                        / ( 1. + RVTMP2 * qtot(i) )
676
677      !--Diagnostic tendencies
678      dqsrim(i)  = - dqlrim  / dtime
679    ENDIF
680
681  ENDIF ! cldfra .GE. seuil_neb
682
683ENDDO ! loop on klon
684
685
686!--Re-calculation of saturation specific humidity
687!--because riming changed temperature
688CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
689
690DO i = 1, klon
691
692  !---------------------------------------------------------
693  !--                  MELTING
694  !---------------------------------------------------------
695  !--Process through which snow melts into rain.
696  !--The formula is homemade.
697  !--NB.: this process needs a temperature adjustment
698
699  !--dqsmelt_max         : maximum snow melting so that temperature
700  !--                      stays higher than 273 K [kg/kg]
701  !--capa_snowflake      : capacitance of a snowflake, equal to
702  !--                      the radius if the snowflake is a sphere [m]
703  !--temp_wetbulb        : wet-bulb temperature [K]
704  !--snow_fallspeed      : snow fall velocity (in clear/cloudy sky) [m/s]
705  !--air_thermal_conduct : thermal conductivity of the air [J/m/K/s]
706  !--gamma_melt          : tuning parameter for melting [-]
707  !--nb_snowflake        : number of snowflakes (in clear/cloudy air) [-]
708
709  IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
710    !--Computed according to
711    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
712    dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
713                        * ( 1. + RVTMP2 * qtot(i) ))
714   
715    !--Initialisation
716    dqsclrmelt = 0.
717    dqscldmelt = 0.
718
719    !--We assume that the snowflakes are spherical
720    capa_snowflake = r_snow
721    !--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
722    air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184   
723    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
724    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
725
726    !--In clear air
727    IF ( ( snowclr(i) .GT. 0. ) .AND.  ( precipfracclr(i) .GT. 0. ) ) THEN
728      !--Formula for the wet-bulb temperature from ECMWF (IFS)
729      !--The vapor used is the vapor in the clear sky
730      temp_wetbulb = temp(i) &
731                   - ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) &
732                   * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
733                   - 40.637 * ( temp(i) - 275. ) )
734      !--Calculated according to
735      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
736      nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr &
737                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
738      dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
739                   * capa_snowflake / RLMLT * gamma_melt &
740                   * MAX(0., temp_wetbulb - RTT) * dtime
741     
742      !--Barrier to limit the melting flux to the clr snow flux in the mesh
743      dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i))
744    ENDIF
745
746
747    !--In cloudy air
748    IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
749      !--As the air is saturated, the wet-bulb temperature is equal to the
750      !--temperature
751      temp_wetbulb = temp(i)
752      !--Calculated according to
753      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
754      nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld &
755                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
756      dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
757                   * capa_snowflake / RLMLT * gamma_melt &
758                   * MAX(0., temp_wetbulb - RTT) * dtime
759
760      !--Barrier to limit the melting flux to the cld snow flux in the mesh
761      dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i))
762   ENDIF
763   
764
765    !--Barrier on temperature. If the total melting flux leads to a
766    !--positive temperature, it is limited to keep temperature above 0 degC.
767    !--It is activated if the total is LOWER than the max
768    !--because everything is negative
769    dqstotmelt = dqsclrmelt + dqscldmelt
770    IF ( dqstotmelt .LT. dqsmelt_max ) THEN
771      !--We redistribute the max melted snow keeping
772      !--the clear/cloud partition of the melted snow
773      dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
774      dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
775      dqstotmelt = dqsmelt_max
776
777    ENDIF
778
779    !--Add tendencies
780    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
781    rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i))
782    raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i))
783    snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i))
784    snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i))
785
786    !--Temperature adjustment with the release of latent
787    !--heat because of melting
788    temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
789                      / ( 1. + RVTMP2 * qtot(i) )
790
791    !--Diagnostic tendencies
792    dqrmelt(i) = - dqstotmelt / dtime
793    dqsmelt(i) = dqstotmelt / dtime
794
795  ENDIF
796
797
798  !---------------------------------------------------------
799  !--                  FREEZING
800  !---------------------------------------------------------
801  !--Process through which rain freezes into snow.
802  !-- We parameterize it as a 2 step process:
803  !--first: freezing following collision with ice crystals
804  !--second: immersion freezing following (inspired by Bigg 1953)
805  !--the latter is parameterized as an exponential decrease of the rain
806  !--water content with a homemade formulya
807  !--This is based on a caracteritic time of freezing, which
808  !--exponentially depends on temperature so that it is
809  !--equal to 1 for temp_nowater (see below) and is close to
810  !--0 for RTT (=273.15 K).
811  !--NB.: this process needs a temperature adjustment
812  !--dqrfreez_max : maximum rain freezing so that temperature
813  !--              stays lower than 273 K [kg/kg]
814  !--tau_freez    : caracteristic time of freezing [s]
815  !--gamma_freez  : tuning parameter [s-1]
816  !--alpha_freez  : tuning parameter for the shape of the exponential curve [-]
817  !--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC)
818
819  IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
820
821 
822    !--1st step: freezing following collision with ice crystals
823    !--Sub-process of freezing which quantifies the collision between
824    !--ice crystals in suspension and falling rain droplets.
825    !--The rain droplets freeze, becoming graupel, and carrying
826    !--the ice crystal (which acted as an ice nucleating particle).
827    !--The formula is adapted from the riming formula.
828    !--it works only in the cloudy part
829
830    dqifreez = 0.
831    dqrtotfreez_step1 = 0.
832
833    IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. &
834         ( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
835      dqrclrfreez = 0.
836      dqrcldfreez = 0.
837
838      !--Computed according to
839      !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
840      dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
841                         * ( 1. + RVTMP2 * qtot(i) ))
842 
843
844      !--The collision efficiency is assumed unity
845      Eff_rain_ice = 1.
846      coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
847      !--Exact version, which does not need a barrier because of
848      !--the exponential decrease.
849      dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
850
851      !--We add the part of rain water that freezes, limited by a temperature barrier
852      !--This quantity is calculated assuming that the number of drop that freeze correspond to the number
853      !--of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
854      !--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3
855      !--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3
856      !--The assumption above corresponds to dNi = dNr, i.e.,
857      !-- dqr = dqi * (4/3 PI rho_rain * r_rain**3) / (4/3 PI rho_ice * r_ice**3)
858      !--Dry density [kg/m3]
859      rho = pplay(i) / temp(i) / RD
860      !--r_ice formula from Sun and Rikus (1999)
861      r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 &
862            + 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
863      dqrcldfreez = dqifreez * rho_rain * r_rain**3. / ( rho_ice * r_ice**3. )
864      dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i))
865      dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max)
866      dqrtotfreez_step1 = dqrcldfreez
867
868      !--Add tendencies
869      !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
870      qice(i) = qice(i) + dqifreez
871      raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
872      snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i))
873      temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD &
874                        / ( 1. + RVTMP2 * qtot(i) )
875
876    ENDIF
877   
878    !-- Second step immersion freezing of rain drops
879    !-- with a homemade timeconstant depending on temperature
880
881    dqrclrfreez = 0.
882    dqrcldfreez = 0.
883    dqrtotfreez_step2 = 0.
884    !--Computed according to
885    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
886
887    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
888                         * ( 1. + RVTMP2 * qtot(i) ))
889
890   
891    tau_freez = 1. / ( beta_freez &
892              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
893
894
895    !--In clear air
896    IF ( rainclr(i) .GT. 0. ) THEN
897      !--Exact solution of dqrain/dt = -qrain/tau_freez
898      dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
899    ENDIF
900
901    !--In cloudy air
902    IF ( raincld(i) .GT. 0. ) THEN
903      !--Exact solution of dqrain/dt = -qrain/tau_freez
904      dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
905    ENDIF
906
907    !--temperature barrier step 2
908    !--It is activated if the total is LOWER than the max
909    !--because everything is negative
910    dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
911 
912    IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
913      !--We redistribute the max freezed rain keeping
914      !--the clear/cloud partition of the freezing rain
915      dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
916      dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
917      dqrtotfreez_step2 = dqrfreez_max
918    ENDIF
919
920
921    !--Add tendencies
922    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
923    rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i))
924    raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
925    snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i))
926    snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i))
927
928
929    !--Temperature adjustment with the uptake of latent
930    !--heat because of freezing
931    temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
932                      / ( 1. + RVTMP2 * qtot(i) )
933
934    !--Diagnostic tendencies
935    dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2         
936    dqrfreez(i) = dqrtotfreez / dtime
937    dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
938
939  ENDIF
940
941
942
943  !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
944  !--we reduce the precipiration fraction in the clear/cloudy air so that the new
945  !--local flux of rain+snow is equal to rain_int_min.
946  !--Here, rain+snow is the gridbox-mean flux of precip.
947  !--Therefore, (rain+snow)/precipfrac is the local flux of precip.
948  !--If the local flux of precip is lower than rain_int_min, i.e.,
949  !-- (rain+snow)/precipfrac < rain_int_min , i.e.,
950  !-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
951  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
952  !--Note that this is physically different than what is proposed in LTP thesis.
953  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
954  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
955
956  !--Calculate outputs
957  rain(i) = rainclr(i) + raincld(i)
958  snow(i) = snowclr(i) + snowcld(i)
959
960  !--Diagnostics
961  !--BEWARE this is indeed a diagnostic: this is an estimation from
962  !--the value of the flux at the bottom interface of the mesh and
963  !--and assuming an upstream numerical calculation
964  !--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is
965  !--used for computing the total ice water content in the mesh
966  !--for radiation only
967  qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr &
968                 + raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld )
969  qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr &
970                 + snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld )
971
972
973ENDDO ! loop on klon
974
975
976END SUBROUTINE poprecip_postcld
977
978END MODULE lmdz_lscp_poprecip
Note: See TracBrowser for help on using the repository browser.