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

Last change on this file since 4818 was 4818, checked in by evignon, 5 months ago

modifications suite à l'atelier nuages du jour
le nouveau schema de precip est full operationnel!
(+ petite commentaire d'un print dans wake)
File size: 29.6 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
66! saturation values
67REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
68! fluxes tendencies because of evaporation
69REAL :: flevapmax, flevapl, flevapi, flevaptot
70! specific humidity tendencies because of evaporation
71REAL :: dqevapl, dqevapi
72! specific heat constant
73REAL :: cpair, cpw
74
75qzero(:)  = 0.0
76dqreva(:) = 0.0
77dqssub(:) = 0.0
78dqevapl=0.0
79dqevapi=0.0
80
81! Calculation of saturation specific humidity
82! depending on temperature:
83CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:))
84! wrt liquid water
85CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
86! wrt ice
87CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
88
89
90
91! First step consists in "thermalizing" the layer:
92! as the flux of precip from layer above "advects" some heat (as the precip is at the temperature
93! of the overlying layer) we recalculate a mean temperature that both the air and the precip in the
94! layer have.
95
96IF (iftop) THEN
97
98   DO i = 1, klon
99      qprecip(i) = 0.
100   ENDDO
101
102ELSE
103
104   DO i = 1, klon
105       ! no condensed water so cp=cp(vapor+dry air)
106       ! RVTMP2=rcpv/rcpd-1
107       cpair=RCPD*(1.0+RVTMP2*qvap(i))
108       cpw=RCPD*RVTMP2
109       ! qprecip has to be thermalized with
110       ! layer's air so that precipitation at the ground has the
111       ! same temperature as the lowermost layer
112       ! we convert the flux into a specific quantity qprecip
113       qprecip(i) = (rain(i)+snow(i))*dtime/((paprsdn(i)-paprsup(i))/RG)
114       ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
115       temp(i) = ( (tempupnew(i))*qprecip(i)*cpw + cpair*temp(i) ) &
116             / (cpair + qprecip(i)*cpw)
117   ENDDO
118
119ENDIF
120
121
122DO i = 1, klon
123
124  ! if precipitation from the layer above
125  IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
126
127    ! Evaporation of liquid precipitation coming from above
128    ! dP/dz=beta*(1-q/qsat)*(P**expo_eva) (lines 1-2)
129    ! multiplying by dz = - dP / g / rho (line 3-4)
130    ! formula from Sundqvist 1988, Klemp & Wilhemson 1978
131    ! LTP: evaporation only in the clear sky part
132
133    flevapl = precipfracclr(i) * coef_eva * (1.0 - qvap(i) / qsatl(i)) &
134            * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva &
135            * temp(i) * RD / pplay(i) &
136            * ( paprsdn(i) - paprsup(i) ) / RG
137
138    ! evaporation is limited by 0 and by the total water amount in
139    ! the precipitation
140    flevapl = MAX(0.0, MIN(flevapl, rainclr(i)))
141
142
143    ! sublimation of the solid precipitation coming from above
144    ! (same formula as for liquid precip)
145    flevapi = precipfracclr(i) * coef_eva_i * (1.0 - qvap(i) / qsati(i)) &
146            * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva_i &
147            * temp(i) * RD / pplay(i) &
148            * ( paprsdn(i) - paprsup(i) ) / RG
149
150    ! sublimation is limited by 0 and by the total water amount in
151    ! the precipitation
152    ! TODO: change max when we will allow for vapor deposition in supersaturated regions
153    flevapi = MAX(0.0, MIN(flevapi, snowclr(i)))
154
155    ! Evaporation limit: we ensure that the layer's fraction below
156    ! the clear sky does not reach saturation. In this case, we
157    ! redistribute the maximum flux flevapmax conserving the ratio liquid/ice
158    ! Max evaporation is computed not to saturate the clear sky precip fraction
159    ! (i.e., the fraction where evaporation occurs)
160    ! It is expressed as a max flux flevapmax
161    !
162    flevapmax = MAX(0.0, ( qsat(i) - qvap(i) ) * precipfracclr(i)) &
163                * ( paprsdn(i) - paprsup(i) ) / RG / dtime
164    flevaptot = flevapl + flevapi
165
166    IF ( flevaptot .GT. flevapmax ) THEN
167        flevapl = flevapmax * flevapl / flevaptot
168        flevapi = flevapmax * flevapi / flevaptot
169    ENDIF
170
171
172    ! New solid and liquid precipitation fluxes after evap and sublimation
173    dqevapl = flevapl / ( paprsdn(i) - paprsup(i) ) * RG * dtime
174    dqevapi = flevapi / ( paprsdn(i) - paprsup(i) ) * RG * dtime
175
176
177    ! vapor is updated after evaporation/sublimation (it is increased)
178    qvap(i) = qvap(i) + dqevapl + dqevapi
179    ! qprecip is the total condensed water in the precip flux (it is decreased)
180    qprecip(i) = qprecip(i) - dqevapl - dqevapi
181    ! air and precip temperature (i.e., gridbox temperature)
182    ! is updated due to latent heat cooling
183    temp(i) = temp(i) &
184            - dqevapl * RLVTT / RCPD &
185            / ( 1.0 + RVTMP2 * ( qvap(i) + qprecip(i) ) ) &
186            - dqevapi * RLSTT / RCPD &
187            / ( 1.0 + RVTMP2 * ( qvap(i) + qprecip(i) ) )
188
189    ! New values of liquid and solid precipitation
190    rainclr(i) = rainclr(i) - flevapl
191    snowclr(i) = snowclr(i) - flevapi
192
193    ! if there is no more precip fluxes, the precipitation fraction in clear
194    ! sky is set to 0
195    IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0.
196
197    ! calculation of the total fluxes
198    rain(i) = rainclr(i) + raincld(i)
199    snow(i) = snowclr(i) + snowcld(i)
200
201  ELSE
202    ! if no precip, we reinitialize the cloud fraction used for the precip to 0
203    precipfraccld(i) = 0.
204    precipfracclr(i) = 0.
205
206  ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. )
207
208
209
210  ! write output tendencies for rain and snow
211
212  dqssub(i)  = -dqevapi/dtime
213  dqreva(i)  = -dqevapl/dtime
214
215ENDDO ! loop on klon
216
217
218END SUBROUTINE poprecip_evapsub
219
220!----------------------------------------------------------------
221! Computes the processes-oriented precipitation formulations for
222! - autoconversion (auto) via a deposition process
223! - aggregation (agg)
224! - riming (rim)
225! - collection (coll)
226! - melting (melt)
227! - freezing (free)
228!
229SUBROUTINE poprecip_postcld( &
230           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
231           temp, qvap, qliq, qice, icefrac, cldfra, &
232           precipfracclr, precipfraccld, &
233           rain, rainclr, raincld, snow, snowclr, snowcld, &
234           dqrauto, dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, &
235           dqsrim, dqsmelt, dqsfreez)
236
237USE lmdz_lscp_ini, ONLY : prt_level, lunout
238USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
239USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
240
241USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb,    &
242                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, &
243                          thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
244                          rho_rain, rho_snow, r_rain, r_snow,                  &
245                          tau_auto_snow_min, tau_auto_snow_max,                &
246                          thresh_precip_frac, eps, air_thermal_conduct,        &
247                          coef_ventil, alpha_freez, gamma_freez, temp_nowater, &
248                          iflag_cloudth_vert, iflag_rain_incloud_vol
249
250
251IMPLICIT NONE
252
253INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
254REAL,    INTENT(IN)                     :: dtime          !--time step [s]
255
256REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
257REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
258REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
259
260REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--
261LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--
262
263REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
264REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity [kg/kg]
265REAL,    INTENT(INOUT), DIMENSION(klon) :: qliq           !--current liquid water specific humidity [kg/kg]
266REAL,    INTENT(INOUT), DIMENSION(klon) :: qice           !--current ice water specific humidity [kg/kg]
267REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--
268REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--
269
270REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
271REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
272                                                          !--NB. at the end of the routine, becomes the fraction of precip
273                                                          !--in the current layer
274
275REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
276REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
277REAL,    INTENT(INOUT), DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
278REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
279REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
280REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
281
282REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrcol         !-- rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
283REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsagg         !-- snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
284REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !-- rain tendency due to autoconversion of cloud liquid [kg/kg/s]
285REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !-- snow tendency due to autoconversion of cloud ice [kg/kg/s]
286REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsrim         !-- snow tendency due to riming [kg/kg/s]
287REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsmelt        !-- snow tendency due to melting [kg/kg/s]
288REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrmelt        !-- rain tendency due to melting [kg/kg/s]
289REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !-- snow tendency due to freezing [kg/kg/s]
290REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !-- rain tendency due to freezing [kg/kg/s]
291
292
293
294!--Local variables
295
296INTEGER :: i
297REAL, DIMENSION(klon) :: hum_to_flux
298REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip
299
300! partition of the fluxes
301REAL :: dcldfra
302REAL :: precipfractot
303REAL :: dprecipfracclr, dprecipfraccld
304REAL :: drainclr, dsnowclr
305REAL :: draincld, dsnowcld
306
307! collection, aggregation and riming
308REAL :: eff_cldfra
309REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain
310REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
311REAL :: dqlcol           ! loss of liquid cloud content due to collection by rain [kg/kg/s]
312REAL :: dqiagg           ! loss of ice cloud content due to collection by aggregation [kg/kg/s]
313REAL :: dqlrim           ! loss of liquid cloud content due to riming on snow[kg/kg/s]
314
315! autoconversion
316REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
317REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
318REAL :: dqlauto          ! loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
319REAL :: dqiauto          ! loss of ice cloud content due to autoconversion to snow [kg/kg/s]
320
321! melting
322REAL :: dqsmelt_max
323REAL :: fallice_clr, fallice_cld
324REAL :: nb_snowflake_clr, nb_snowflake_cld
325REAL :: capa_snowflake, temp_wetbulb
326REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
327REAL, DIMENSION(klon) :: qzero, qsat, dqsat
328
329! freezing
330REAL :: dqrfreez_max
331REAL :: tau_freez
332REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez
333
334
335!--Initialisation of variables
336
337qzero(:)    = 0.
338
339dqrcol(:)   = 0.
340dqsagg(:)   = 0.
341dqsauto(:)  = 0.
342dqrauto(:)  = 0.
343dqsrim(:)   = 0.
344dqrmelt(:)  = 0.
345dqsmelt(:)  = 0.
346dqrfreez(:) = 0.
347dqsfreez(:) = 0.
348
349
350DO i = 1, klon
351
352  ! variables initialisation
353  dqlrim  = 0.
354  dqlcol  = 0.
355  dqiagg  = 0.
356  dqiauto = 0.
357  dqlauto = 0.
358
359  !--hum_to_flux: coef to convert a specific quantity to a flux
360  !-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt
361  hum_to_flux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
362  qtot(i) = qvap(i) + qliq(i) + qice(i) &
363          + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / hum_to_flux(i)
364
365  !------------------------------------------------------------
366  !--     PRECIPITATION FRACTIONS UPDATE
367  !------------------------------------------------------------
368  !--The goal of this routine is to reattribute precipitation fractions
369  !--and fluxes to clear or cloudy air, depending on the variation of
370  !--the cloud fraction on the vertical dimension. We assume a
371  !--maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
372  !--and LTP thesis, 2021)
373  !--NB. in fact, we assume a maximum-random overlap of the total precip. frac
374
375  !--Initialisation
376  precipfractot = precipfracclr(i) + precipfraccld(i)
377
378  !--Instead of using the cloud cover which was use in LTP thesis, we use the
379  !--total precip. fraction to compute the maximum-random overlap. This is
380  !--because all the information of the cloud cover is embedded into
381  !--precipfractot, and this allows for taking into account the potential
382  !--reduction of the precipitation fraction because either the flux is too
383  !--small (see barrier at the end of poprecip_postcld) or the flux is completely
384  !--evaporated (see barrier at the end of poprecip_precld)
385  !--NB. precipfraccld(i) is here the cloud fraction of the layer above
386  precipfractot = 1. - ( 1. - precipfractot ) * &
387                 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
388               / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
389
390
391  !--precipfraccld(i) is here the cloud fraction of the layer above
392  dcldfra = cldfra(i) - precipfraccld(i)
393  !--Tendency of the clear-sky precipitation fraction. We add a MAX on the
394  !--calculation of the current CS precip. frac.
395  dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
396  !--Tendency of the cloudy precipitation fraction. We add a MAX on the
397  !--calculation of the current CS precip. frac.
398  !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) )
399  !--We remove it, because cldfra is guaranteed to be > O (the MAX is activated
400  !--if cldfra < 0)
401  dprecipfraccld = dcldfra
402
403
404  !--If the cloud extends
405  IF ( dprecipfraccld .GT. 0. ) THEN
406    !--If there is no CS precip, nothing happens.
407    !--If there is, we reattribute some of the CS precip flux
408    !--to the cloud precip flux, proportionnally to the
409    !--decrease of the CS precip fraction
410    IF ( precipfracclr(i) .LE. 0. ) THEN
411      drainclr = 0.
412      dsnowclr = 0.
413    ELSE
414      drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i)
415      dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i)
416    ENDIF
417  !--If the cloud narrows
418  ELSEIF ( dprecipfraccld .LT. 0. ) THEN
419    !--We reattribute some of the cloudy precip flux
420    !--to the CS precip flux, proportionnally to the
421    !--decrease of the cloud precip fraction
422    draincld = dprecipfraccld / precipfraccld(i) * raincld(i)
423    dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i)
424    drainclr = - draincld
425    dsnowclr = - dsnowcld
426  !--If the cloud stays the same or if there is no cloud above and
427  !--in the current layer, nothing happens
428  ELSE
429    drainclr = 0.
430    dsnowclr = 0.
431  ENDIF
432
433  !--We add the tendencies
434  precipfraccld(i) = precipfraccld(i) + dprecipfraccld
435  precipfracclr(i) = precipfracclr(i) + dprecipfracclr
436  rainclr(i) = rainclr(i) + drainclr
437  snowclr(i) = snowclr(i) + dsnowclr
438  raincld(i) = raincld(i) - drainclr
439  snowcld(i) = snowcld(i) - dsnowclr
440 
441
442  ! if vertical heterogeneity is taken into account, we use
443  ! the "true" volume fraction instead of a modified
444  ! surface fraction (which is larger and artificially
445  ! reduces the in-cloud water).
446  IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
447    eff_cldfra = ctot_vol(i)
448  ELSE
449    eff_cldfra = cldfra(i)
450  ENDIF
451
452
453  ! Start precipitation growth processes
454
455  !--If the cloud is big enough, the precipitation processes activate
456  IF ( cldfra(i) .GE. seuil_neb ) THEN
457
458    !---------------------------------------------------------
459    !--           COLLECTION AND AGGREGATION
460    !---------------------------------------------------------
461    !--Collection: processus through which rain collects small liquid droplets
462    !--in suspension, and add it to the rain flux
463    !--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
464    !--Those processes are treated before autoconversion because we do not
465    !--want to collect/aggregate the newly formed fluxes, which already
466    !--"saw" the cloud as they come from it
467    !--gamma_col: tuning coefficient [-]
468    !--rho_rain: volumic mass of rain [kg/m3]
469    !--r_rain: size of the rain droplets [m]
470    !--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
471    !--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
472    !--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
473    !--get in-cloud mean quantities. The two divisions by eff_cldfra are
474    !--then simplified.
475    Eff_rain_liq = 1.
476    coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
477    IF ((raincld(i) .GT. 0.) .AND. (coef_col .GT. 0.)) THEN
478      !--Explicit version
479      !dqlcol = - coef_col * qliq(i) * raincld(i) / precipfraccld(i) *dtime
480      !--Semi-implicit version
481      !dqlcol = qliq(i) * ( 1. / ( 1. + coef_col * raincld(i) / precipfraccld(i)*dtime ) - 1. )
482      !--Implicit version
483      !qrain = raincld(i) / hum_to_flux(i)
484      !coef_tmp = coef_col * dtime * ( qrain / precipfraccld(i) + qliq(i) / eff_cldfra )
485      !dqlcol = qliq(i) * ( 1. / ( 1. + 0.5 * ( coef_tmp - 1. + SQRT( &
486      !    ( 1. - coef_tmp )**2. + 4. * coef_col * dtime * qrain / precipfraccld(i) ) &
487      !                   ) ) - 1. )
488      !--Barriers so that the processes do not consume more liquid/ice than
489      !--available.
490      !dqlcol = MAX( - qliq(i), dqlcol )
491      !--Exact version
492      dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
493
494      !--Add tendencies
495      qliq(i) = qliq(i) + dqlcol
496      raincld(i) = raincld(i) - dqlcol * hum_to_flux(i)
497
498      !--Outputs
499      dqrcol(i)  = - dqlcol  / dtime
500    ENDIF
501
502    !--Same as for aggregation
503    !--Following Milbrandt and Yau 2005, it s a product of a collection
504    !--efficiency and a sticking efficiency
505    Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
506    coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
507    IF ((snowcld(i) .GT. 0.) .AND. (coef_agg .GT. 0.)) THEN
508      !--Explicit version
509      !dqiagg = - coef_agg &
510      !      * qice(i) * snowcld(i) / precipfraccld(i) * dtime
511      !--Barriers so that the processes do not consume more liquid/ice than
512      !--available.
513      !dqiagg = MAX( - qice(i), dqiagg )
514      !--Exact version
515      dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
516
517      !--Add tendencies
518      qice(i) = qice(i) + dqiagg
519      snowcld(i) = snowcld(i) - dqiagg * hum_to_flux(i)
520
521      !--Outputs
522      dqsagg(i)  = - dqiagg  / dtime
523    ENDIF
524
525
526    !---------------------------------------------------------
527    !--                  AUTOCONVERSION
528    !---------------------------------------------------------
529
530    ! TODO
531    IF ( ptconv(i) ) THEN ! if convective point
532        qthresh_auto_rain = cld_lc_con
533        qthresh_auto_snow = cld_lc_con
534
535        tau_auto_rain = cld_tau_con
536        tau_auto_snow = tau_auto_snow_max &
537                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
538
539        expo_auto_rain = cld_expo_con
540        expo_auto_snow = cld_expo_con
541    ELSE
542        qthresh_auto_rain = cld_lc_lsc
543        qthresh_auto_snow = cld_lc_lsc
544
545        tau_auto_rain = cld_tau_lsc
546        tau_auto_snow = tau_auto_snow_max &
547                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
548
549        expo_auto_rain = cld_expo_lsc
550        expo_auto_snow = cld_expo_lsc
551    ENDIF
552
553
554    ! Liquid water quantity to remove according to (Sundqvist, 1978)
555    ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
556    !.........................................................
557    ! we first treat the second term (with exponential) in an explicit way
558    ! and then treat the first term (-q/tau) in an exact way
559
560    dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
561               - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
562
563    dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
564               - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
565
566
567    dqlauto = MAX( - qliq(i), dqlauto )
568    dqiauto = MAX( - qice(i), dqiauto )
569
570    qliq(i) = qliq(i) + dqlauto
571    qice(i) = qice(i) + dqiauto
572
573    raincld(i) = raincld(i) - dqlauto * hum_to_flux(i)
574    snowcld(i) = snowcld(i) - dqiauto * hum_to_flux(i)
575
576    !--Outputs
577    dqsauto(i) = - dqiauto / dtime
578    dqrauto(i) = - dqlauto / dtime
579
580
581    ! FOLLOWING PROCESSES IMPLY A PHASE CHANGE SO A TEMPERATURE
582    ! ADJUSTMENT
583
584    !---------------------------------------------------------
585    !--                  RIMING
586    !---------------------------------------------------------
587
588    !--Following Seifert and Beheng 2006, assuming a cloud droplet diameter
589    !--of 20 microns.
590    Eff_snow_liq = 0.2
591    coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
592    IF ((snowcld(i) .GT. 0.) .AND. (coef_rim .GT. 0.)) THEN
593      !--Explicit version
594      !dqlrim = - gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq &
595      !       * qliq(i) * snowcld(i) /  precipfraccld(i) * dtime
596      !--Barriers so that the processes do not consume more liquid than
597      !--available.
598      !dqlrim = MAX( - qliq(i), dqlrim )
599      !--Exact version
600      dqlrim = qliq(i) * ( EXP( - dtime * coef_col * snowcld(i) / precipfraccld(i) ) - 1. )
601
602      qliq(i) = qliq(i) + dqlrim
603      snowcld(i) = snowcld(i) - dqlrim * hum_to_flux(i)
604
605      ! Latent heat of melting with precipitation thermalization
606      temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
607                        / ( 1. + RVTMP2 * qtot(i) )
608
609      !--Outputs
610      dqsrim(i)  = - dqlrim  / dtime
611    ENDIF
612
613  ENDIF ! rneb .GE. seuil_neb
614
615ENDDO ! iteration on klon
616
617
618! Calculation of saturation specific humidity
619! depending on temperature:
620CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
621
622DO i = 1, klon
623
624  !---------------------------------------------------------
625  !--                  MELTING
626  !---------------------------------------------------------
627
628  IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
629    dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
630                        * ( 1. + RVTMP2 * qtot(i) ))
631   
632    dqsclrmelt = 0.
633    dqscldmelt = 0.
634
635    capa_snowflake = r_snow
636    ! ATTENTION POUR L'INSTANT C'EST UN WBULB SELON FORMULE ECMWF
637    ! ATTENTION EST-CE QVAP ?????
638    temp_wetbulb = temp(i) - ( qsat(i) - qvap(i) ) &
639                 * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
640                 - 40.637 * ( temp(i) - 275. ) )
641
642    IF ( snowclr(i) .GT. 0. ) THEN
643      ! ATTENTION ATTENTION ATTENTION
644      fallice_clr = 1.
645      nb_snowflake_clr = snowclr(i) / precipfracclr(i) / fallice_clr &
646                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
647      dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
648                   * capa_snowflake / RLMLT * coef_ventil &
649                   * MAX(0., temp_wetbulb - RTT) * dtime
650    ENDIF
651
652    IF ( snowcld(i) .GT. 0. ) THEN
653      ! ATTENTION ATTENTION ATTENTION
654      fallice_cld = 1.
655      nb_snowflake_cld = snowcld(i) / precipfraccld(i) / fallice_cld &
656                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
657      dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
658                   * capa_snowflake / RLMLT * coef_ventil &
659                   * MAX(0., temp_wetbulb - RTT) * dtime
660    ENDIF
661   
662    ! barrier
663    ! lower than bec. negative values
664    dqstotmelt = dqsclrmelt + dqscldmelt
665    IF ( dqstotmelt .LT. dqsmelt_max ) THEN
666      dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
667      dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
668      dqstotmelt = dqsmelt_max
669    ENDIF
670
671    ! update of rainfall and snowfall due to melting
672    rainclr(i) = rainclr(i) - dqsclrmelt * hum_to_flux(i)
673    raincld(i) = raincld(i) - dqscldmelt * hum_to_flux(i)
674    snowclr(i) = snowclr(i) + dqsclrmelt * hum_to_flux(i)
675    snowcld(i) = snowcld(i) + dqscldmelt * hum_to_flux(i)
676
677    ! Latent heat of melting with precipitation thermalization
678    temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
679                      / ( 1. + RVTMP2 * qtot(i) )
680
681    dqrmelt(i) = - dqstotmelt / dtime
682    dqsmelt(i) = dqstotmelt / dtime
683
684  ENDIF
685
686
687  !---------------------------------------------------------
688  !--                  FREEZING
689  !---------------------------------------------------------
690
691  IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
692
693    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
694                         * ( 1. + RVTMP2 * qtot(i) ))
695   
696    tau_freez = 1. / ( gamma_freez &
697              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
698    dqrtotfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
699
700    ! barrier
701    ! max bec. those are negative values
702    dqrtotfreez = MAX(dqrtotfreez, dqrfreez_max)
703
704    ! partition between clear and cloudy air
705    ! proportionnal to the rain fluxes in clear / cloud
706    dqrclrfreez = dqrtotfreez * rainclr(i) / ( rainclr(i) + raincld(i) )
707    dqrcldfreez = dqrtotfreez - dqrclrfreez
708
709    ! update of rainfall and snowfall due to melting
710    rainclr(i) = rainclr(i) + dqrclrfreez * hum_to_flux(i)
711    raincld(i) = raincld(i) + dqrcldfreez * hum_to_flux(i)
712    snowclr(i) = snowclr(i) - dqrclrfreez * hum_to_flux(i)
713    snowcld(i) = snowcld(i) - dqrcldfreez * hum_to_flux(i)
714
715    ! Latent heat of melting with precipitation thermalization
716    temp(i) = temp(i) - dqrtotfreez * RLMLT / RCPD &
717                      / ( 1. + RVTMP2 * qtot(i) )
718
719    dqrfreez(i) = dqrtotfreez / dtime
720    dqsfreez(i) = - dqrtotfreez / dtime
721
722  ENDIF
723
724
725  !! MISE A JOUR DES FRACTIONS PRECIP CLD et CS
726  ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
727
728  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
729  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
730
731  rain(i) = rainclr(i) + raincld(i)
732  snow(i) = snowclr(i) + snowcld(i)
733
734! write output tendencies for rain and snow
735
736ENDDO
737
738
739
740END SUBROUTINE poprecip_postcld
741
742END MODULE lmdz_lscp_poprecip
Note: See TracBrowser for help on using the repository browser.