source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_lscp_poprecip.F90 @ 5119

Last change on this file since 5119 was 5103, checked in by abarral, 4 months ago

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