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

Last change on this file since 4832 was 4832, checked in by evignon, 3 months ago

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