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

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

modifications du commit precedent a la suite de l'atelier nuages
File size: 34.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_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_eva_i, expo_eva, expo_eva_i, 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
139    draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) &
140             * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva &
141             * temp(i) * RD / pplay(i) &
142             * ( paprsdn(i) - paprsup(i) ) / RG
143
144    !--Evaporation is limited by 0 and by the total water amount in
145    !--the precipitation
146    draineva = MIN(0., MAX(draineva, -rainclr(i)))
147
148
149    !--Sublimation of the solid precipitation coming from above
150    !--(same formula as for liquid precip)
151    dsnowsub = - precipfracclr(i) * coef_eva_i * (1. - qvap(i) / qsati(i)) &
152             * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva_i &
153             * temp(i) * RD / pplay(i) &
154             * ( paprsdn(i) - paprsup(i) ) / RG
155
156    !--Sublimation is limited by 0 and by the total water amount in
157    !--the precipitation
158    ! TODO: change max when we will allow for vapor deposition in supersaturated regions
159    dsnowsub = MIN(0., MAX(dsnowsub, -snowclr(i)))
160
161    !--Evaporation limit: we ensure that the layer's fraction below
162    !--the clear sky does not reach saturation. In this case, we
163    !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
164    !--Max evaporation is computed not to saturate the clear sky precip fraction
165    !--(i.e., the fraction where evaporation occurs)
166    !--It is expressed as a max flux dprecip_evasub_max
167   
168    dprecip_evasub_max = MIN(0., ( qvap(i) - qsat(i) ) * precipfracclr(i)) &
169                     * hum_to_flux(i)
170    dprecip_evasub_tot = draineva + dsnowsub
171
172    !--Barriers
173    !--If activates if the total is LOWER than the max because
174    !--everything is negative
175    IF ( dprecip_evasub_tot .LT. dprecip_evasub_max ) THEN
176      draineva = dprecip_evasub_max * draineva / dprecip_evasub_tot
177      dsnowsub = dprecip_evasub_max * dsnowsub / dprecip_evasub_tot
178    ENDIF
179
180
181    !--New solid and liquid precipitation fluxes after evap and sublimation
182    dqrevap = draineva / hum_to_flux(i)
183    dqssubl = dsnowsub / hum_to_flux(i)
184
185
186    !--Vapor is updated after evaporation/sublimation (it is increased)
187    qvap(i) = qvap(i) - dqrevap - dqssubl
188    !--qprecip is the total condensed water in the precip flux (it is decreased)
189    qprecip(i) = qprecip(i) + dqrevap + dqssubl
190    !--Air and precip temperature (i.e., gridbox temperature)
191    !--is updated due to latent heat cooling
192    temp(i) = temp(i) &
193            + dqrevap * RLVTT / RCPD &
194            / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) &
195            + dqssubl * RLSTT / RCPD &
196            / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
197
198    !--Add tendencies
199    rainclr(i) = rainclr(i) + draineva
200    snowclr(i) = 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 (coll)
232! - melting (melt)
233! - freezing (free)
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           qrain, qsnow, 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,                  &
251                          tau_auto_snow_min, tau_auto_snow_max,                &
252                          thresh_precip_frac, eps, air_thermal_conduct,        &
253                          coef_ventil, alpha_freez, gamma_freez, temp_nowater, &
254                          iflag_cloudth_vert, iflag_rain_incloud_vol
255
256
257IMPLICIT NONE
258
259INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
260REAL,    INTENT(IN)                     :: dtime          !--time step [s]
261
262REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
263REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
264REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
265
266REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
267LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
268
269REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
270REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity [kg/kg]
271REAL,    INTENT(INOUT), DIMENSION(klon) :: qliq           !--current liquid water specific humidity [kg/kg]
272REAL,    INTENT(INOUT), DIMENSION(klon) :: qice           !--current ice water specific humidity [kg/kg]
273REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--ice fraction [-]
274REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction [-]
275
276REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
277REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
278                                                          !--NB. at the end of the routine, becomes the fraction of precip
279                                                          !--in the current layer
280
281REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
282REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
283REAL,    INTENT(INOUT), DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
284REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
285REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
286REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
287
288REAL,    INTENT(OUT),   DIMENSION(klon) :: qrain          !--specific rain content [kg/kg]
289REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnow          !--specific snow content [kg/kg]
290REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
291REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
292REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
293REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
294REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsrim         !--snow tendency due to riming [kg/kg/s]
295REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
296REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
297REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
298REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
299
300
301
302!--Local variables
303
304INTEGER :: i
305REAL, DIMENSION(klon) :: hum_to_flux
306REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip
307
308! partition of the fluxes
309REAL :: dcldfra
310REAL :: precipfractot
311REAL :: dprecipfracclr, dprecipfraccld
312REAL :: drainclr, dsnowclr
313REAL :: draincld, dsnowcld
314
315! collection, aggregation and riming
316REAL :: eff_cldfra
317REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
318REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
319REAL :: dqlcol           ! loss of liquid cloud content due to collection by rain [kg/kg/s]
320REAL :: dqiagg           ! loss of ice cloud content due to collection by aggregation [kg/kg/s]
321REAL :: dqlrim           ! loss of liquid cloud content due to riming on snow[kg/kg/s]
322
323! autoconversion
324REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
325REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
326REAL :: dqlauto          ! loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
327REAL :: dqiauto          ! loss of ice cloud content due to autoconversion to snow [kg/kg/s]
328
329! melting
330REAL :: dqsmelt_max
331REAL :: fallice_clr, fallice_cld
332REAL :: nb_snowflake_clr, nb_snowflake_cld
333REAL :: capa_snowflake, temp_wetbulb
334REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
335REAL, DIMENSION(klon) :: qzero, qsat, dqsat
336
337! freezing
338REAL :: dqrfreez_max
339REAL :: tau_freez
340REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez
341
342
343!--Initialisation of variables
344
345qzero(:)    = 0.
346
347dqrcol(:)   = 0.
348dqsagg(:)   = 0.
349dqsauto(:)  = 0.
350dqrauto(:)  = 0.
351dqsrim(:)   = 0.
352dqrmelt(:)  = 0.
353dqsmelt(:)  = 0.
354dqrfreez(:) = 0.
355dqsfreez(:) = 0.
356
357
358DO i = 1, klon
359
360  ! variables initialisation
361  dqlrim  = 0.
362  dqlcol  = 0.
363  dqiagg  = 0.
364  dqiauto = 0.
365  dqlauto = 0.
366
367  !--hum_to_flux: coef to convert a specific quantity to a flux
368  !-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt
369  hum_to_flux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
370  qtot(i) = qvap(i) + qliq(i) + qice(i) &
371          + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / hum_to_flux(i)
372
373  !------------------------------------------------------------
374  !--     PRECIPITATION FRACTIONS UPDATE
375  !------------------------------------------------------------
376  !--The goal of this routine is to reattribute precipitation fractions
377  !--and fluxes to clear or cloudy air, depending on the variation of
378  !--the cloud fraction on the vertical dimension. We assume a
379  !--maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
380  !--and LTP thesis, 2021)
381  !--NB. in fact, we assume a maximum-random overlap of the total precip. frac
382
383  !--Initialisation
384  precipfractot = precipfracclr(i) + precipfraccld(i)
385
386  !--Instead of using the cloud cover which was use in LTP thesis, we use the
387  !--total precip. fraction to compute the maximum-random overlap. This is
388  !--because all the information of the cloud cover is embedded into
389  !--precipfractot, and this allows for taking into account the potential
390  !--reduction of the precipitation fraction because either the flux is too
391  !--small (see barrier at the end of poprecip_postcld) or the flux is completely
392  !--evaporated (see barrier at the end of poprecip_precld)
393  !--NB. precipfraccld(i) is here the cloud fraction of the layer above
394  precipfractot = 1. - ( 1. - precipfractot ) * &
395                 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
396               / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
397
398
399  !--precipfraccld(i) is here the cloud fraction of the layer above
400  dcldfra = cldfra(i) - precipfraccld(i)
401  !--Tendency of the clear-sky precipitation fraction. We add a MAX on the
402  !--calculation of the current CS precip. frac.
403  dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
404  !--Tendency of the cloudy precipitation fraction. We add a MAX on the
405  !--calculation of the current CS precip. frac.
406  !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) )
407  !--We remove it, because cldfra is guaranteed to be > O (the MAX is activated
408  !--if cldfra < 0)
409  dprecipfraccld = dcldfra
410
411
412  !--If the cloud extends
413  IF ( dprecipfraccld .GT. 0. ) THEN
414    !--If there is no CS precip, nothing happens.
415    !--If there is, we reattribute some of the CS precip flux
416    !--to the cloud precip flux, proportionnally to the
417    !--decrease of the CS precip fraction
418    IF ( precipfracclr(i) .LE. 0. ) THEN
419      drainclr = 0.
420      dsnowclr = 0.
421    ELSE
422      drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i)
423      dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i)
424    ENDIF
425  !--If the cloud narrows
426  ELSEIF ( dprecipfraccld .LT. 0. ) THEN
427    !--We reattribute some of the cloudy precip flux
428    !--to the CS precip flux, proportionnally to the
429    !--decrease of the cloud precip fraction
430    draincld = dprecipfraccld / precipfraccld(i) * raincld(i)
431    dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i)
432    drainclr = - draincld
433    dsnowclr = - dsnowcld
434  !--If the cloud stays the same or if there is no cloud above and
435  !--in the current layer, nothing happens
436  ELSE
437    drainclr = 0.
438    dsnowclr = 0.
439  ENDIF
440
441  !--We add the tendencies
442  precipfraccld(i) = precipfraccld(i) + dprecipfraccld
443  precipfracclr(i) = precipfracclr(i) + dprecipfracclr
444  rainclr(i) = rainclr(i) + drainclr
445  snowclr(i) = snowclr(i) + dsnowclr
446  raincld(i) = raincld(i) - drainclr
447  snowcld(i) = snowcld(i) - dsnowclr
448 
449
450  ! if vertical heterogeneity is taken into account, we use
451  ! the "true" volume fraction instead of a modified
452  ! surface fraction (which is larger and artificially
453  ! reduces the in-cloud water).
454  IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
455    eff_cldfra = ctot_vol(i)
456  ELSE
457    eff_cldfra = cldfra(i)
458  ENDIF
459
460
461  ! Start precipitation growth processes
462
463  !--If the cloud is big enough, the precipitation processes activate
464  IF ( cldfra(i) .GE. seuil_neb ) THEN
465
466    !---------------------------------------------------------
467    !--           COLLECTION AND AGGREGATION
468    !---------------------------------------------------------
469    !--Collection: processus through which rain collects small liquid droplets
470    !--in suspension, and add it to the rain flux
471    !--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
472    !--Those processes are treated before autoconversion because we do not
473    !--want to collect/aggregate the newly formed fluxes, which already
474    !--"saw" the cloud as they come from it
475    !--The formulas come from Muench and Lohmann 2020
476
477    !--gamma_col: tuning coefficient [-]
478    !--rho_rain: volumic mass of rain [kg/m3]
479    !--r_rain: size of the rain droplets [m]
480    !--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
481    !--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
482    !--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
483    !--get in-cloud mean quantities. The two divisions by eff_cldfra are
484    !--then simplified.
485    Eff_rain_liq = 1.
486    coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
487    IF ((raincld(i) .GT. 0.) .AND. (coef_col .GT. 0.)) THEN
488      !-- ATTENTION Double implicit version
489      !qrain_tmp = raincld(i) / hum_to_flux(i)
490      !coef_tmp = coef_col * dtime * ( qrain_tmp / precipfraccld(i) + qliq(i) / eff_cldfra )
491      !dqlcol = qliq(i) * ( 1. / ( 1. + 0.5 * ( coef_tmp - 1. + SQRT( &
492      !    ( 1. - coef_tmp )**2. + 4. * coef_col * dtime * qrain_tmp / precipfraccld(i) ) &
493      !                   ) ) - 1. )
494      !--Barriers so that the processes do not consume more liquid/ice than
495      !--available.
496      !dqlcol = MAX( - qliq(i), dqlcol )
497      !--Exact version, which does not need a barrier because of
498      !--the exponential decrease
499      dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
500
501      !--Add tendencies
502      qliq(i) = qliq(i) + dqlcol
503      raincld(i) = raincld(i) - dqlcol * hum_to_flux(i)
504
505      !--Diagnostic tendencies
506      dqrcol(i)  = - dqlcol  / dtime
507    ENDIF
508
509    !--Same as for aggregation
510    !--Eff_snow_liq formula: following Milbrandt and Yau 2005,
511    !--it s a product of a collection efficiency and a sticking efficiency
512    Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
513    coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
514    IF ((snowcld(i) .GT. 0.) .AND. (coef_agg .GT. 0.)) THEN
515      !-- ATTENTION Double implicit version?
516      !--Barriers so that the processes do not consume more liquid/ice than
517      !--available.
518      !dqiagg = MAX( - qice(i), dqiagg )
519      !--Exact version, which does not need a barrier because of
520      !--the exponential decrease
521      dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
522
523      !--Add tendencies
524      qice(i) = qice(i) + dqiagg
525      snowcld(i) = snowcld(i) - dqiagg * hum_to_flux(i)
526
527      !--Diagnostic tendencies
528      dqsagg(i)  = - dqiagg  / dtime
529    ENDIF
530
531
532    !---------------------------------------------------------
533    !--                  AUTOCONVERSION
534    !---------------------------------------------------------
535    !--Autoconversion converts liquid droplets/ice crystals into
536    !--rain drops/snowflakes. It relies on the formulations by
537    !--Sundqvist 1978.
538    ! TODO what else?
539
540    !--If we are in a convective point, we have different parameters
541    !--for the autoconversion
542    IF ( ptconv(i) ) THEN
543        ! ATTENTION voir les constantes ici qui ne devraient pas etre identiques
544        qthresh_auto_rain = cld_lc_con
545        qthresh_auto_snow = cld_lc_con
546
547        tau_auto_rain = cld_tau_con
548        !--tau for snow depends on the ice fraction in mixed-phase clouds     
549        tau_auto_snow = tau_auto_snow_max &
550                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
551
552        expo_auto_rain = cld_expo_con
553        expo_auto_snow = cld_expo_con
554    ELSE
555        ! ATTENTION voir les constantes ici qui ne devraient pas etre identiques
556        qthresh_auto_rain = cld_lc_lsc
557        qthresh_auto_snow = cld_lc_lsc
558
559        tau_auto_rain = cld_tau_lsc
560        !--tau for snow depends on the ice fraction in mixed-phase clouds     
561        tau_auto_snow = tau_auto_snow_max &
562                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
563
564        expo_auto_rain = cld_expo_lsc
565        expo_auto_snow = cld_expo_lsc
566    ENDIF
567
568
569    ! Liquid water quantity to remove according to (Sundqvist, 1978)
570    ! dqliq/dt = -qliq/tau * ( 1-exp(-(qcin/clw)**2) )
571    !.........................................................
572    ! we first treat the second term (with exponential) in an explicit way
573    ! and then treat the first term (-q/tau) in an exact way
574
575    dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
576               - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
577
578    dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
579               - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
580
581
582    !--Barriers so that we don t create more rain/snow
583    !--than there is liquid/ice
584    dqlauto = MAX( - qliq(i), dqlauto )
585    dqiauto = MAX( - qice(i), dqiauto )
586
587    !--Add tendencies
588    qliq(i) = qliq(i) + dqlauto
589    qice(i) = qice(i) + dqiauto
590    raincld(i) = raincld(i) - dqlauto * hum_to_flux(i)
591    snowcld(i) = snowcld(i) - dqiauto * hum_to_flux(i)
592
593    !--Diagnostic tendencies
594    dqsauto(i) = - dqiauto / dtime
595    dqrauto(i) = - dqlauto / dtime
596
597
598    !---------------------------------------------------------
599    !--                  RIMING
600    !---------------------------------------------------------
601    !--Process which converts liquid droplets in suspension into
602    !--snow (graupel in fact) because of the collision between
603    !--those and falling snowflakes.
604    !--The formula come from Muench and Lohmann 2020
605    !--NB.: this process needs a temperature adjustment
606
607    !--Eff_snow_liq formula: following Seifert and Beheng 2006,
608    !--assuming a cloud droplet diameter of 20 microns.
609    Eff_snow_liq = 0.2
610    coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
611    IF ((snowcld(i) .GT. 0.) .AND. (coef_rim .GT. 0.)) THEN
612      !-- ATTENTION Double implicit version?
613      !--Barriers so that the processes do not consume more liquid than
614      !--available.
615      !dqlrim = MAX( - qliq(i), dqlrim )
616      !--Exact version, which does not need a barrier because of
617      !--the exponential decrease
618      dqlrim = qliq(i) * ( EXP( - dtime * coef_col * snowcld(i) / precipfraccld(i) ) - 1. )
619
620      !--Add tendencies
621      qliq(i) = qliq(i) + dqlrim
622      snowcld(i) = snowcld(i) - dqlrim * hum_to_flux(i)
623
624      !--Temperature adjustment with the release of latent
625      !--heat because of solid condensation
626      temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
627                        / ( 1. + RVTMP2 * qtot(i) )
628
629      !--Diagnostic tendencies
630      dqsrim(i)  = - dqlrim  / dtime
631    ENDIF
632
633  ! ATTENTION veut on faire un processus similaire et symetrique qui
634  ! convertirait la pluie en surfusion qui tombe sur des cristaux de glace
635  ! en flux de neige ? (si la temperature est negative)
636
637  ENDIF ! cldfra .GE. seuil_neb
638
639ENDDO ! loop on klon
640
641
642!--Re-calculation of saturation specific humidity
643!--because riming changed temperature
644CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
645
646DO i = 1, klon
647
648  !---------------------------------------------------------
649  !--                  MELTING
650  !---------------------------------------------------------
651  !--Process through which snow melts into rain.
652  !--The formula is homemade.
653  !--NB.: this process needs a temperature adjustment
654
655  !--dqsmelt_max: maximum snow melting so that temperature
656  !--             stays higher than 273 K [kg/kg]
657  !--capa_snowflake: capacitance of a snowflake, equal to
658  !--                the radius if the snowflake is a sphere [m]
659  !--temp_wetbulb: wet-bulb temperature [K]
660  !--fallice: snow fall velocity (in clear/cloudy sky) [m/s]
661  !--air_thermal_conduct: thermal conductivity of the air [J/m/K/s]
662  !--coef_ventil: ventilation coefficient [-]
663  !--nb_snowflake: number of snowflakes (in clear/cloudy air) [-]
664
665  IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
666    !--Computed according to
667    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
668    dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
669                        * ( 1. + RVTMP2 * qtot(i) ))
670   
671    !--Initialisation
672    dqsclrmelt = 0.
673    dqscldmelt = 0.
674
675    capa_snowflake = r_snow
676    ! ATTENTION POUR L'INSTANT C'EST UN WBULB SELON FORMULE ECMWF
677    ! ATTENTION EST-CE QVAP ?????
678    temp_wetbulb = temp(i) - ( qsat(i) - qvap(i) ) &
679                 * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
680                 - 40.637 * ( temp(i) - 275. ) )
681
682    !--In clear air
683    IF ( snowclr(i) .GT. 0. ) THEN
684      ! ATTENTION fallice a definir
685      fallice_clr = 1.
686      !--Calculated according to
687      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
688      nb_snowflake_clr = snowclr(i) / precipfracclr(i) / fallice_clr &
689                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
690      dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
691                   * capa_snowflake / RLMLT * coef_ventil &
692                   * MAX(0., temp_wetbulb - RTT) * dtime
693    ENDIF
694
695    !--In cloudy air
696    IF ( snowcld(i) .GT. 0. ) THEN
697      ! ATTENTION fallice a definir
698      fallice_cld = 1.
699      !--Calculated according to
700      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
701      nb_snowflake_cld = snowcld(i) / precipfraccld(i) / fallice_cld &
702                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
703      dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
704                   * capa_snowflake / RLMLT * coef_ventil &
705                   * MAX(0., temp_wetbulb - RTT) * dtime
706    ENDIF
707   
708    !--Barriers
709    !--It is activated if the total is LOWER than the max
710    !--because everything is negative
711    dqstotmelt = dqsclrmelt + dqscldmelt
712    IF ( dqstotmelt .LT. dqsmelt_max ) THEN
713      !--We redistribute the max melted snow keeping
714      !--the clear/cloud partition of the melted snow
715      dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
716      dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
717      dqstotmelt = dqsmelt_max
718    ENDIF
719
720    !--Add tendencies
721    rainclr(i) = rainclr(i) - dqsclrmelt * hum_to_flux(i)
722    raincld(i) = raincld(i) - dqscldmelt * hum_to_flux(i)
723    snowclr(i) = snowclr(i) + dqsclrmelt * hum_to_flux(i)
724    snowcld(i) = snowcld(i) + dqscldmelt * hum_to_flux(i)
725
726    !--Temperature adjustment with the release of latent
727    !--heat because of melting
728    temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
729                      / ( 1. + RVTMP2 * qtot(i) )
730
731    !--Diagnostic tendencies
732    dqrmelt(i) = - dqstotmelt / dtime
733    dqsmelt(i) = dqstotmelt / dtime
734
735  ENDIF
736
737
738  !---------------------------------------------------------
739  !--                  FREEZING
740  !---------------------------------------------------------
741  !--Process through which rain freezes into snow. This is
742  !--parameterized as an exponential decrease of the rain
743  !--water content.
744  !--The formula is homemade.
745  !--This is based on a caracteritic time of freezing, which
746  !--exponentially depends on temperature so that it is
747  !--equal to 1 for temp_nowater (see below) and is close to
748  !--0 for RTT (=273.15 K).
749  !--NB.: this process needs a temperature adjustment
750
751  !--dqrfreez_max: maximum rain freezing so that temperature
752  !--              stays lower than 273 K [kg/kg]
753  !--tau_freez: caracteristic time of freezing [s]
754  !--gamma_freez: tuning parameter [s-1]
755  !--alpha_freez: tuning parameter for the shape of the exponential curve [-]
756  !--temp_nowater: temperature below which no liquid water exists [K] (about -40 degC)
757
758  IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
759
760    !--Computed according to
761    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
762    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
763                         * ( 1. + RVTMP2 * qtot(i) ))
764 
765    tau_freez = 1. / ( gamma_freez &
766              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
767    !--Exact solution of dqrain/dt = -qrain/tau_freez
768    dqrtotfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
769
770    !--Barrier
771    !--It is a MAX because everything is negative
772    dqrtotfreez = MAX(dqrtotfreez, dqrfreez_max)
773
774    !--The partition between clear and cloudy air
775    !--is proportionnal to the rain fluxes in clear / cloudy air
776    dqrclrfreez = dqrtotfreez * rainclr(i) / ( rainclr(i) + raincld(i) )
777    dqrcldfreez = dqrtotfreez - dqrclrfreez
778
779    !--Add tendencies
780    rainclr(i) = rainclr(i) + dqrclrfreez * hum_to_flux(i)
781    raincld(i) = raincld(i) + dqrcldfreez * hum_to_flux(i)
782    snowclr(i) = snowclr(i) - dqrclrfreez * hum_to_flux(i)
783    snowcld(i) = snowcld(i) - dqrcldfreez * hum_to_flux(i)
784
785    !--Temperature adjustment with the uptake of latent
786    !--heat because of freezing
787    temp(i) = temp(i) - dqrtotfreez * RLMLT / RCPD &
788                      / ( 1. + RVTMP2 * qtot(i) )
789
790    !--Diagnostic tendencies
791    dqrfreez(i) = dqrtotfreez / dtime
792    dqsfreez(i) = - dqrtotfreez / dtime
793
794  ENDIF
795
796  !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
797  !--we reduce the precipiration fraction in the clear/cloudy air so that the new
798  !--local flux of rain+snow is equal to rain_int_min.
799  !--Here, rain+snow is the gridbox-mean flux of precip.
800  !--Therefore, (rain+snow)/precipfrac is the local flux of precip.
801  !--If the local flux of precip is lower than rain_int_min, i.e.,
802  !-- (rain+snow)/precipfrac < rain_int_min , i.e.,
803  !-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
804  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
805  !--Note that this is physically different than what is proposed in LTP thesis.
806  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
807  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
808
809  !--Calculate outputs
810  rain(i) = rainclr(i) + raincld(i)
811  snow(i) = snowclr(i) + snowcld(i)
812
813  !--Diagnostics
814  ! ATTENTION A REPRENDRE
815  qrain(i) = rain(i) / hum_to_flux(i)
816  qsnow(i) = snow(i) / hum_to_flux(i)
817
818ENDDO ! loop on klon
819
820END SUBROUTINE poprecip_postcld
821
822END MODULE lmdz_lscp_poprecip
Note: See TracBrowser for help on using the repository browser.