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

Last change on this file since 4898 was 4898, checked in by evignon, 5 weeks ago

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