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

Last change on this file since 4816 was 4809, checked in by evignon, 22 months ago

mise a jour des routines proprecip suite à l'atelier nuage
File size: 25.4 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,dqsrim,dqsmelt,dqsfreez)
235
236USE lmdz_lscp_ini, ONLY : prt_level, lunout
237USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
238USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
239
240USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb,    &
241                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, &
242                          thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
243                          rho_rain, rho_snow, r_rain, r_snow, Eff_rain_liq,    &
244                          Eff_snow_ice, Eff_snow_liq, tau_auto_snow_min,       &
245                          tau_auto_snow_max, thresh_precip_frac, eps,          &
246                          iflag_cloudth_vert, iflag_rain_incloud_vol
247
248IMPLICIT NONE
249
250INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
251REAL,    INTENT(IN)                     :: dtime          !--time step [s]
252
253REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
254REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
255REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
256
257REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--
258LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--
259
260REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
261REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity [kg/kg]
262REAL,    INTENT(INOUT), DIMENSION(klon) :: qliq           !--current liquid water specific humidity [kg/kg]
263REAL,    INTENT(INOUT), DIMENSION(klon) :: qice           !--current ice water specific humidity [kg/kg]
264REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--
265REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--
266
267REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
268REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
269                                                          !--NB. at the end of the routine, becomes the fraction of precip
270                                                          !--in the current layer
271
272REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
273REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
274REAL,    INTENT(INOUT), DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
275REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
276REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
277REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
278
279REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrcol         !-- rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
280REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsagg         !-- snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
281REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !-- rain tendency due to autoconversion of cloud liquid [kg/kg/s]
282REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !-- snow tendency due to autoconversion of cloud ice [kg/kg/s]
283REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsrim         !-- snow tendency due to riming [kg/kg/s]
284REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsmelt        !-- snow tendency due to melting [kg/kg/s]
285REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrmelt        !-- rain tendency due to melting [kg/kg/s]
286REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !-- snow tendency due to freezing [kg/kg/s]
287REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !-- rain tendency due to freezing [kg/kg/s]
288
289
290
291!--Local variables
292
293INTEGER :: i
294
295REAL :: hum_to_flux
296REAL :: dcldfra
297REAL :: precipfractot
298REAL :: dprecipfracclr, dprecipfraccld
299REAL :: drainclr, dsnowclr
300REAL :: draincld, dsnowcld
301REAL :: eff_cldfra
302REAL :: coef_col, coef_agg, coef_tmp, qrain
303REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
304REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
305REAL ::  dqlcol           ! loss of liquid cloud content due to collection by rain [kg/kg/s]
306REAL ::  dqiagg           ! loss of ice cloud content due to collection by aggregation [kg/kg/s]
307REAL ::  dqlauto          ! loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
308REAL ::  dqiauto          ! loss of ice cloud content due to autoconversion to snow [kg/kg/s]
309REAL ::  dqlrim           ! loss of liquid cloud content due to riming on snow[kg/kg/s]
310
311
312!--Initialisation of variables
313
314dqrcol(:) = 0.
315dqsagg(:) = 0.
316dqsauto(:) = 0.
317dqrauto(:) = 0.
318dqsrim(:) = 0.
319dqrmelt(:) = 0.
320dqsmelt(:) = 0.
321dqrfreez(:) = 0.
322dqsfreez(:) = 0.
323
324
325DO i = 1, klon
326
327
328  ! variables initialisation
329  dqlrim = 0.0
330  dqlcol = 0.0
331  dqiagg = 0.0
332  dqiauto = 0.0
333  dqlauto = 0.0
334
335  !------------------------------------------------------------
336  !--     PRECIPITATION FRACTIONS UPDATE
337  !------------------------------------------------------------
338  !--The goal of this routine is to reattribute precipitation fractions
339  !--and fluxes to clear or cloudy air, depending on the variation of
340  !--the cloud fraction on the vertical dimension. We assume a
341  !--maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
342  !--and LTP thesis, 2021)
343  !--NB. in fact, we assume a maximum-random overlap of the total precip. frac
344
345  !--Initialisation
346  !--hum_to_flux: coef to convert a specific quantity to a flux
347  !-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt
348  hum_to_flux = ( paprsdn(i) - paprsup(i) ) / RG / dtime
349  precipfractot = precipfracclr(i) + precipfraccld(i)
350
351  !--Instead of using the cloud cover which was use in LTP thesis, we use the
352  !--total precip. fraction to compute the maximum-random overlap. This is
353  !--because all the information of the cloud cover is embedded into
354  !--precipfractot, and this allows for taking into account the potential
355  !--reduction of the precipitation fraction because either the flux is too
356  !--small (see barrier at the end of poprecip_postcld) or the flux is completely
357  !--evaporated (see barrier at the end of poprecip_precld)
358  !--NB. precipfraccld(i) is here the cloud fraction of the layer above
359  precipfractot = 1. - ( 1. - precipfractot ) * &
360                 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
361               / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
362
363
364  !--precipfraccld(i) is here the cloud fraction of the layer above
365  dcldfra = cldfra(i) - precipfraccld(i)
366  !--Tendency of the clear-sky precipitation fraction. We add a MAX on the
367  !--calculation of the current CS precip. frac.
368  dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
369  !--Tendency of the cloudy precipitation fraction. We add a MAX on the
370  !--calculation of the current CS precip. frac.
371  !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) )
372  !--We remove it, because cldfra is guaranteed to be > O (the MAX is activated
373  !--if cldfra < 0)
374  dprecipfraccld = dcldfra
375
376
377  !--If the cloud extends
378  IF ( dprecipfraccld .GT. 0. ) THEN
379    !--If there is no CS precip, nothing happens.
380    !--If there is, we reattribute some of the CS precip flux
381    !--to the cloud precip flux, proportionnally to the
382    !--decrease of the CS precip fraction
383    IF ( precipfracclr(i) .LE. 0. ) THEN
384      drainclr = 0.
385      dsnowclr = 0.
386    ELSE
387      drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i)
388      dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i)
389    ENDIF
390  !--If the cloud narrows
391  ELSEIF ( dprecipfraccld .LT. 0. ) THEN
392    !--We reattribute some of the cloudy precip flux
393    !--to the CS precip flux, proportionnally to the
394    !--decrease of the cloud precip fraction
395    draincld = dprecipfraccld / precipfraccld(i) * raincld(i)
396    dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i)
397    drainclr = - draincld
398    dsnowclr = - dsnowcld
399  !--If the cloud stays the same or if there is no cloud above and
400  !--in the current layer, nothing happens
401  ELSE
402    drainclr = 0.
403    dsnowclr = 0.
404  ENDIF
405
406  !--We add the tendencies
407  precipfraccld(i) = precipfraccld(i) + dprecipfraccld
408  precipfracclr(i) = precipfracclr(i) + dprecipfracclr
409  rainclr(i) = rainclr(i) + drainclr
410  snowclr(i) = snowclr(i) + dsnowclr
411  raincld(i) = raincld(i) - drainclr
412  snowcld(i) = snowcld(i) - dsnowclr
413 
414
415  ! if vertical heterogeneity is taken into account, we use
416  ! the "true" volume fraction instead of a modified
417  ! surface fraction (which is larger and artificially
418  ! reduces the in-cloud water).
419  IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
420    eff_cldfra = ctot_vol(i)
421  ELSE
422    eff_cldfra = cldfra(i)
423  ENDIF
424
425
426  ! Start precipitation growth processes
427
428  !--If the cloud is big enough, the precipitation processes activate
429  IF ( cldfra(i) .GE. seuil_neb ) THEN
430
431    !---------------------------------------------------------
432    !--           COLLECTION AND AGGREGATION
433    !---------------------------------------------------------
434    !--Collection: processus through which rain collects small liquid droplets
435    !--in suspension, and add it to the rain flux
436    !--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
437    !--Those processes are treated before autoconversion because we do not
438    !--want to collect/aggregate the newly formed fluxes, which already
439    !--"saw" the cloud as they come from it
440    !--gamma_col: tuning coefficient [-]
441    !--rho_rain: volumic mass of rain [kg/m3]
442    !--r_rain: size of the rain droplets [m]
443    !--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
444    !--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
445    !--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
446    !--get in-cloud mean quantities. The two divisions by eff_cldfra are
447    !--then simplified.
448
449    coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
450    IF ((raincld(i) .GT. 0.) .AND. (coef_col .GT. 0.)) THEN
451       !--Explicit version
452       !dqlcol = - coef_col * qliq(i) * raincld(i) / precipfraccld(i) *dtime
453       !--Semi-implicit version
454       !dqlcol = qliq(i) * ( 1. / ( 1. + coef_col * raincld(i) / precipfraccld(i)*dtime ) - 1. )
455       !--Implicit version
456       !qrain = raincld(i) / hum_to_flux
457       !coef_tmp = coef_col * dtime *( qrain / precipfraccld(i) + qliq(i) / eff_cldfra )
458       !dqlcol = qliq(i) * ( 1. / ( 1. + 0.5 * ( coef_tmp - 1. + SQRT( &
459       !    ( 1. - coef_tmp )**2. + 4. * coef_col * dtime *qrain / precipfraccld(i) ) &
460       !                   ) ) - 1. )
461       !dqlcol=max(dqlcol,-qliq(i))
462       ! Exact version
463       dqlcol=qliq(i)*(exp(-dtime * coef_col * raincld(i) / precipfraccld(i))-1.)
464    ENDIF
465
466    !--Same as for aggregation
467    coef_agg=gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
468    IF ((snowcld(i) .GT. 0.) .AND. (coef_agg .GT. 0.)) THEN
469        !--Explicit version     
470        !dqiagg = - coef_agg &
471        !      * qice(i) * snowcld(i) / precipfraccld(i) * dtime
472        ! Exact version
473        dqiagg=qice(i)*(exp(-dtime * coef_agg * snowcld(i) / precipfraccld(i))-1.)
474    ENDIF
475    !--Barriers so that the processes do not consume more liquid/ice than
476    !--available.
477    dqlcol = MAX( - qliq(i), dqlcol )
478    dqiagg = MAX( - qice(i), dqiagg )
479
480    !--Add tendencies
481    qliq(i) = qliq(i) + dqlcol
482    qice(i) = qice(i) + dqiagg
483    raincld(i) = raincld(i) - dqlcol * hum_to_flux
484    snowcld(i) = snowcld(i) - dqiagg * hum_to_flux
485
486
487    !---------------------------------------------------------
488    !--                  AUTOCONVERSION
489    !---------------------------------------------------------
490
491    ! TODO
492    IF ( ptconv(i) ) THEN ! if convective point
493        qthresh_auto_rain = cld_lc_con
494        qthresh_auto_snow = cld_lc_con
495
496        tau_auto_rain = cld_tau_con
497        tau_auto_snow = tau_auto_snow_max &
498                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
499
500        expo_auto_rain = cld_expo_con
501        expo_auto_snow = cld_expo_con
502    ELSE
503        qthresh_auto_rain = cld_lc_lsc
504        qthresh_auto_snow = cld_lc_lsc
505
506        tau_auto_rain = cld_tau_lsc
507        tau_auto_snow = tau_auto_snow_max &
508                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
509
510        expo_auto_rain = cld_expo_lsc
511        expo_auto_snow = cld_expo_lsc
512    ENDIF
513
514
515    ! Liquid water quantity to remove according to (Sundqvist, 1978)
516    ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
517    !.........................................................
518    ! we first treat the second term (with exponential) in an explicit way
519    ! and then treat the first term (-q/tau) in an exact way
520
521    dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
522               - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
523
524    dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
525               - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
526
527
528    dqlauto = MAX( - qliq(i), dqlauto )
529    dqiauto = MAX( - qice(i), dqiauto )
530
531    qliq(i) = qliq(i) + dqlauto
532    qice(i) = qice(i) + dqiauto
533
534    raincld(i) = raincld(i) - dqlauto * hum_to_flux
535    snowcld(i) = snowcld(i) - dqiauto * hum_to_flux
536
537
538    ! FOLLOWING PROCESSES IMPLY A PHASE CHANGE SO A TEMPERATURE
539    ! ADJUSTMENT
540
541    !---------------------------------------------------------
542    !--                  RIMING
543    !---------------------------------------------------------
544
545    dqlrim=0.0
546
547    ! remplacer la premiere ligne par "coef_rim" ?
548    IF (snowcld(i) .GT. 0.) THEN
549       dqlrim = - gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq &
550              * qliq(i) * snowcld(i) /  precipfraccld(i) * dtime
551    ENDIF
552    dqlrim = MAX( - qliq(i), dqlrim )
553
554    qliq(i) = qliq(i) + dqlrim
555
556    snowcld(i) = snowcld(i) - dqlrim * hum_to_flux
557
558
559  ENDIF ! rneb .GE. seuil_neb
560
561
562
563  !---------------------------------------------------------
564  !--                  FREEZING
565  !---------------------------------------------------------
566
567  !dqrfree_max = MINOUMAX(0., ( RTT - temp(i) ) / RLMLT * RCPD / ( 1. + RVTMP2 * ( qtot(i) + qprecip(i) ) ))
568
569  !---------------------------------------------------------
570  !--                  MELTING
571  !---------------------------------------------------------
572
573  !flux = velocity * N_0 * 4. / 3. * PI * r_snow**3. * rho_snow
574
575  !IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
576  !  dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD / ( 1. + RVTMP2 * ( qtot(i) + qprecip(i) ) ))
577  !  dsnowtotmelt_max = dqsmelt_max * hum_to_flux(i)
578  ! 
579  !  dsnowtotmelt = - nb_snowflake * 4. * PI * mol_diff_vap * snowflake_capa / RLMLT * coef_ventil &
580  !              * MAX(0., ticebulb - RTT) &
581  !              * ( paprsdn(i) - paprsup(i) ) / RG
582  !  ! max bec. negative values
583  !  dsnowtotmelt = MAX(dsnowtotmelt, dsnowmelt_max)
584  !  dsnowclrmelt = dsnowtotmelt * snowclr(i) / ( snowclr(i) + snowcld(i) )
585  !  dsnowcldmelt = dsnowtotmelt - dsnowclrmelt
586
587
588  !  ! update of rainfall and snowfall due to melting
589  !  rainclr(i) = rainclr(i) - dsnowclrmelt(i)
590  !  raincld(i) = raincld(i) - dsnowcldmelt(i)
591  !  snowclr(i) = snowclr(i) + dsnowclrmelt(i)
592  !  snowcld(i) = snowcld(i) + dsnowcldmelt(i)
593  !ENDIF
594  !
595  !! Latent heat of melting with precipitation thermalization
596  !zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
597  !*RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
598
599
600  !! MISE A JOUR DES FRACTIONS PRECIP CLD et CS
601  ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
602
603  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
604  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
605
606  rain(i) = rainclr(i) + raincld(i)
607  snow(i) = snowclr(i) + snowcld(i)
608
609! write output tendencies for rain and snow
610
611  dqsrim(i)  = -dqlrim/dtime
612  dqrcol(i)  = -dqlcol/dtime
613  dqsagg(i)  = -dqiagg/dtime
614  dqsauto(i) = -dqiauto/dtime
615  dqrauto(i) = -dqlauto/dtime
616
617
618
619ENDDO
620
621
622
623END SUBROUTINE poprecip_postcld
624
625END MODULE lmdz_lscp_poprecip
Note: See TracBrowser for help on using the repository browser.