source: LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90 @ 5408

Last change on this file since 5408 was 5408, checked in by evignon, 7 hours ago

debut de l'externalisation des processus de precipitation de lmdz_lscp avec ici, les processus precedents
la formation des nuages.
A Borella, E Vignon
File size: 53.5 KB
Line 
1MODULE lmdz_lscp_precip
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! historical (till CMIP6) version of the pre-cloud formation
16! precipitation scheme containing precip evaporation and melting
17
18SUBROUTINE histprecip_precld( &
19           klon, dtime, iftop, paprsdn, paprsup, pplay, zt, ztupnew, zq, &
20           zmqc, zneb, znebprecipclr, &
21           zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, dqreva, dqssub &
22           )
23
24USE lmdz_lscp_ini, ONLY : iflag_evap_prec
25USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, ztfondue
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) :: zt             !--current temperature [K]
42REAL,    INTENT(IN),    DIMENSION(klon) :: ztupnew        !--updated temperature of the overlying layer [K]
43
44REAL,    INTENT(INOUT), DIMENSION(klon) :: zq             !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
45REAL,    INTENT(INOUT), DIMENSION(klon) :: zmqc           !--specific humidity in the precipitation falling from the upper layer [kg/kg]
46
47REAL,    INTENT(IN),    DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
48REAL,    INTENT(INOUT), DIMENSION(klon) :: znebprecipclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
49
50REAL,    INTENT(INOUT), DIMENSION(klon) :: zrfl           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
51REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
52REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflcld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
53REAL,    INTENT(INOUT), DIMENSION(klon) :: zifl           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
54REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
55REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--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
61REAL :: zmair
62REAL :: zcpair, zcpeau
63REAL, DIMENSION(klon) :: znebprecip
64REAL, DIMENSION(klon) :: qzero
65REAL, DIMENSION(klon) :: zqs, zdqs
66REAL, DIMENSION(klon) :: qsl, qsi
67REAL, DIMENSION(klon) :: dqsl, dqsi
68
69REAL :: zqev0
70REAL :: zqev, zqevt
71REAL :: zqevi, zqevti
72
73REAL, DIMENSION(klon) :: zrfln, zifln
74
75REAL :: zmelt
76INTEGER :: i
77
78qzero(:) = 0.
79
80
81! --------------------------------------------------------------------
82! P1> Thermalization of precipitation falling from the overlying layer
83! --------------------------------------------------------------------
84! Computes air temperature variation due to enthalpy transported by
85! precipitation. Precipitation is then thermalized with the air in the
86! layer.
87! The precipitation should remain thermalized throughout the different
88! thermodynamical transformations.
89! The corresponding water mass should
90! be added when calculating the layer's enthalpy change with
91! temperature
92! ---------------------------------------------------------------------
93
94IF (iftop) THEN
95
96  DO i = 1, klon
97    zmqc(i) = 0.
98  ENDDO
99
100ELSE
101
102  DO i = 1, klon
103
104    zmair=(paprsdn(i)-paprsup(i))/RG
105    ! no condensed water so cp=cp(vapor+dry air)
106    ! RVTMP2=rcpv/rcpd-1
107    zcpair=RCPD*(1.0+RVTMP2*zq(i))
108    zcpeau=RCPD*RVTMP2
109
110    ! zmqc: precipitation mass that has to be thermalized with
111    ! layer's air so that precipitation at the ground has the
112    ! same temperature as the lowermost layer
113    zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
114    ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
115    zt(i) = ( ztupnew(i)*zmqc(i)*zcpeau + zcpair*zt(i) ) &
116          / (zcpair + zmqc(i)*zcpeau)
117
118  ENDDO
119
120ENDIF
121
122! --------------------------------------------------------------------
123! P2> Precipitation evaporation/sublimation/melting
124! --------------------------------------------------------------------
125! A part of the precipitation coming from above is evaporated/sublimated/melted.
126! --------------------------------------------------------------------
127
128IF (iflag_evap_prec.GE.1) THEN
129
130  ! Calculation of saturation specific humidity
131  ! depending on temperature:
132  CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,0,.false.,zqs,zdqs)
133  ! wrt liquid water
134  CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,1,.false.,qsl,dqsl)
135  ! wrt ice
136  CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,2,.false.,qsi,dqsi)
137
138  DO i = 1, klon
139
140    ! if precipitation
141    IF (zrfl(i)+zifl(i).GT.0.) THEN
142
143    ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
144    ! c_iso: likely important to distinguish cs from neb isotope precipitation
145
146      IF (iflag_evap_prec.GE.4) THEN
147        zrfl(i) = zrflclr(i)
148        zifl(i) = ziflclr(i)
149      ENDIF
150
151      IF (iflag_evap_prec.EQ.1) THEN
152        znebprecip(i)=zneb(i)
153      ELSE
154        znebprecip(i)=MAX(zneb(i),znebprecip(i))
155      ENDIF
156
157      IF (iflag_evap_prec.GT.4) THEN
158      ! Max evaporation not to saturate the clear sky precip fraction
159      ! i.e. the fraction where evaporation occurs
160        zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
161      ELSEIF (iflag_evap_prec .EQ. 4) THEN
162      ! Max evaporation not to saturate the whole mesh
163      ! Pay attention -> lead to unrealistic and excessive evaporation
164        zqev0 = MAX(0.0, zqs(i)-zq(i))
165      ELSE
166      ! Max evap not to saturate the fraction below the cloud
167        zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
168      ENDIF
169
170      ! Evaporation of liquid precipitation coming from above
171      ! dP/dz=beta*(1-q/qsat)*sqrt(P)
172      ! formula from Sundquist 1988, Klemp & Wilhemson 1978
173      ! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
174
175      IF (iflag_evap_prec.EQ.3) THEN
176        zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
177          *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
178          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
179      ELSE IF (iflag_evap_prec.GE.4) THEN
180        zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
181          *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
182          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
183      ELSE
184        zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
185          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
186      ENDIF
187
188      zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * RG*dtime/(paprsdn(i)-paprsup(i))
189
190      ! sublimation of the solid precipitation coming from above
191      IF (iflag_evap_prec.EQ.3) THEN
192        zqevti = znebprecip(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
193          *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
194          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
195      ELSE IF (iflag_evap_prec.GE.4) THEN
196        zqevti = znebprecipclr(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
197          *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
198          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
199      ELSE
200        zqevti = 1.*coef_sub*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
201          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
202      ENDIF
203
204      zqevti = MAX(0.0,MIN(zqevti,zifl(i))) * RG*dtime/(paprsdn(i)-paprsup(i))
205
206      ! A. JAM
207      ! Evaporation limit: we ensure that the layer's fraction below
208      ! the cloud or the whole mesh (depending on iflag_evap_prec)
209      ! does not reach saturation. In this case, we
210      ! redistribute zqev0 conserving the ratio liquid/ice
211
212      IF (zqevt+zqevti.GT.zqev0) THEN
213        zqev=zqev0*zqevt/(zqevt+zqevti)
214        zqevi=zqev0*zqevti/(zqevt+zqevti)
215      ELSE
216        zqev=zqevt
217        zqevi=zqevti
218      ENDIF
219
220      !--Diagnostics
221      dqreva(i) = - zqev / dtime
222      dqssub(i) = - zqevti / dtime
223
224      ! New solid and liquid precipitation fluxes after evap and sublimation
225      zrfln(i) = Max(0.,zrfl(i) - zqev*(paprsdn(i)-paprsup(i))/RG/dtime)
226      zifln(i) = Max(0.,zifl(i) - zqevi*(paprsdn(i)-paprsup(i))/RG/dtime)
227     
228
229      ! vapor, temperature, precip fluxes update
230      ! vapor is updated after evaporation/sublimation (it is increased)
231      zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
232        * (RG/(paprsdn(i)-paprsup(i)))*dtime
233      ! zmqc is the total condensed water in the precip flux (it is decreased)
234      zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
235        * (RG/(paprsdn(i)-paprsup(i)))*dtime
236      ! air and precip temperature (i.e., gridbox temperature)
237      ! is updated due to latent heat cooling
238      zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
239        * (RG/(paprsdn(i)-paprsup(i)))*dtime    &
240        * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
241        + (zifln(i)-zifl(i))                      &
242        * (RG/(paprsdn(i)-paprsup(i)))*dtime    &
243        * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
244
245      ! New values of liquid and solid precipitation
246      zrfl(i) = zrfln(i)
247      zifl(i) = zifln(i)
248
249      ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
250      !                         due to evap + sublim                                     
251                                 
252
253      IF (iflag_evap_prec.GE.4) THEN
254        zrflclr(i) = zrfl(i)
255        ziflclr(i) = zifl(i)
256        IF(zrflclr(i) + ziflclr(i).LE.0) THEN
257          znebprecipclr(i) = 0.0
258        ENDIF
259        zrfl(i) = zrflclr(i) + zrflcld(i)
260        zifl(i) = ziflclr(i) + ziflcld(i)
261      ENDIF
262
263      ! c_iso duplicate for isotopes or loop on isotopes
264
265      ! Melting:
266      zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
267      ! precip fraction that is melted
268      zmelt = MIN(MAX(zmelt,0.),1.)
269
270      ! update of rainfall and snowfall due to melting
271      IF (iflag_evap_prec.GE.4) THEN
272        zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
273        zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
274        zrfl(i)=zrflclr(i)+zrflcld(i)
275      ELSE
276        zrfl(i)=zrfl(i)+zmelt*zifl(i)
277      ENDIF
278
279
280      ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
281
282      ! Latent heat of melting because of precipitation melting
283      ! NB: the air + precip temperature is simultaneously updated
284      zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprsdn(i)-paprsup(i)) &
285        *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
286     
287      IF (iflag_evap_prec.GE.4) THEN
288        ziflclr(i)=ziflclr(i)*(1.-zmelt)
289        ziflcld(i)=ziflcld(i)*(1.-zmelt)
290        zifl(i)=ziflclr(i)+ziflcld(i)
291      ELSE
292        zifl(i)=zifl(i)*(1.-zmelt)
293      ENDIF
294
295    ELSE
296      ! if no precip, we reinitialize the cloud fraction used for the precip to 0
297      ! AB note that this assignment is useless, as znebprecip is not re-used
298      znebprecip(i)=0.
299
300    ENDIF ! (zrfl(i)+zifl(i).GT.0.)
301
302  ENDDO ! loop on klon
303
304ENDIF ! (iflag_evap_prec>=1)
305
306END SUBROUTINE histprecip_precld
307
308
309!----------------------------------------------------------------
310! Computes the processes-oriented precipitation formulations for
311! evaporation and sublimation
312!
313SUBROUTINE poprecip_precld( &
314           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
315           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
316           rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub &
317           )
318
319USE lmdz_lscp_ini, ONLY : prt_level, lunout
320USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
321USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
322USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub
323USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
324
325IMPLICIT NONE
326
327
328INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
329REAL,    INTENT(IN)                     :: dtime          !--time step [s]
330LOGICAL, INTENT(IN)                     :: iftop          !--if top of the column
331
332
333REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
334REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
335REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
336
337REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
338REAL,    INTENT(IN),    DIMENSION(klon) :: tempupnew      !--updated temperature of the overlying layer [K]
339
340REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
341REAL,    INTENT(INOUT), DIMENSION(klon) :: qprecip        !--specific humidity in the precipitation falling from the upper layer [kg/kg]
342
343REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
344REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
345
346REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
347REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
348
349REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
350REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
351REAL,    INTENT(IN),    DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
352REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
353REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
354REAL,    INTENT(IN),    DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
355
356REAL,    INTENT(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
357REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
358
359
360!--Integer for interating over klon
361INTEGER :: i
362!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
363REAL, DIMENSION(klon) :: dhum_to_dflux
364!--
365REAL, DIMENSION(klon) :: rho, dz
366
367!--Saturation values
368REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
369!--Vapor in the clear sky
370REAL :: qvapclr
371!--Fluxes tendencies because of evaporation and sublimation
372REAL :: dprecip_evasub_max, draineva, dsnowsub, dprecip_evasub_tot
373!--Specific humidity tendencies because of evaporation and sublimation
374REAL :: dqrevap, dqssubl
375!--Specific heat constant
376REAL :: cpair, cpw
377
378!--Initialisation
379qzero(:)  = 0.
380dqreva(:) = 0.
381dqssub(:) = 0.
382dqrevap   = 0.
383dqssubl   = 0.
384
385!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
386dhum_to_dflux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime
387rho(:) = pplay(:) / temp(:) / RD
388dz(:) = ( paprsdn(:) - paprsup(:) ) / RG / rho(:)
389
390!--Calculation of saturation specific humidity
391!--depending on temperature:
392CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:))
393!--wrt liquid water
394CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
395!--wrt ice
396CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
397
398
399
400!--First step consists in "thermalizing" the layer:
401!--as the flux of precip from layer above "advects" some heat (as the precip is at the temperature
402!--of the overlying layer) we recalculate a mean temperature that both the air and the precip in the
403!--layer have.
404
405IF (iftop) THEN
406
407  DO i = 1, klon
408    qprecip(i) = 0.
409  ENDDO
410
411ELSE
412
413  DO i = 1, klon
414    !--No condensed water so cp=cp(vapor+dry air)
415    !-- RVTMP2=rcpv/rcpd-1
416    cpair = RCPD * ( 1. + RVTMP2 * qvap(i) )
417    cpw = RCPD * RVTMP2
418    !--qprecip has to be thermalized with
419    !--layer's air so that precipitation at the ground has the
420    !--same temperature as the lowermost layer
421    !--we convert the flux into a specific quantity qprecip
422    qprecip(i) = ( rain(i) + snow(i) ) / dhum_to_dflux(i)
423    !-- t(i,k+1) + d_t(i,k+1): new temperature of the overlying layer
424    temp(i) = ( tempupnew(i) * qprecip(i) * cpw + cpair * temp(i) ) &
425            / ( cpair + qprecip(i) * cpw )
426  ENDDO
427
428ENDIF
429
430! TODO Probleme : on utilise qvap total dans la maille pour l'evap / sub
431! alors qu'on n'evap / sub que dans le ciel clair
432! deux options pour cette routine :
433! - soit on diagnostique le nuage AVANT l'evap / sub et on estime donc
434! la fraction precipitante ciel clair dans la maille, ce qui permet de travailler
435! avec des fractions, des fluxs et surtout un qvap dans le ciel clair
436! - soit on pousse la param de Ludo au bout, et on prend un qvap de k+1
437! dans le ciel clair, avec un truc comme :
438!   qvapclr(k) = qvapclr(k+1)/qtot(k+1) * qtot(k)
439! UPDATE : on code la seconde version. A voir si on veut mettre la premiere version.
440
441
442DO i = 1, klon
443
444  !--If there is precipitation from the layer above
445  ! NOTE TODO here we could replace the condition on precipfracclr(i) by a condition
446  ! such as eps or thresh_precip_frac, to remove the senseless barrier in the formulas
447  ! of evap / sublim
448  IF ( ( ( rain(i) + snow(i) ) .GT. 0. ) .AND. ( precipfracclr(i) .GT. 0. ) ) THEN
449
450    IF ( ok_corr_vap_evasub ) THEN
451      !--Corrected version - we use the same water ratio between
452      !--the clear and the cloudy sky as in the layer above. This
453      !--extends the assumption that the cloud fraction is the same
454      !--as the layer above. This is assumed only for the evap / subl
455      !--process
456      !--Note that qvap(i) is the total water in the gridbox, and
457      !--precipfraccld(i) is the cloud fraction in the layer above
458      qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
459    ELSE
460      !--Legacy version from Ludo - we use the total specific humidity
461      !--for the evap / subl process
462      qvapclr = qvap(i)
463    ENDIF
464
465    !--Evaporation of liquid precipitation coming from above
466    !--in the clear sky only
467    !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva)
468    !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
469    !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly)
470    !--which does not need a barrier on rainclr, because included in the formula
471    draineva = precipfracclr(i) * ( MAX(0., &
472             - coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) &
473             + ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_eva ) &
474               ) )**( 1. / ( 1. - expo_eva ) ) - rainclr(i)
475             
476    !--Evaporation is limited by 0
477    draineva = MIN(0., draineva)
478
479
480    !--Sublimation of the solid precipitation coming from above
481    !--(same formula as for liquid precip)
482    !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly)
483    !--which does not need a barrier on snowclr, because included in the formula
484    dsnowsub = precipfracclr(i) * ( MAX(0., &
485             - coef_sub * ( 1. - expo_sub ) * (1. - qvapclr / qsati(i)) * dz(i) &
486             + ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_sub ) &
487             ) )**( 1. / ( 1. - expo_sub ) ) - snowclr(i)
488
489    !--Sublimation is limited by 0
490    ! TODO: change max when we will allow for vapor deposition in supersaturated regions
491    dsnowsub = MIN(0., dsnowsub)
492
493    !--Evaporation limit: we ensure that the layer's fraction below
494    !--the clear sky does not reach saturation. In this case, we
495    !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
496    !--Max evaporation is computed not to saturate the clear sky precip fraction
497    !--(i.e., the fraction where evaporation occurs)
498    !--It is expressed as a max flux dprecip_evasub_max
499   
500    dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr(i)) &
501                     * dhum_to_dflux(i)
502    dprecip_evasub_tot = draineva + dsnowsub
503
504    !--Barriers
505    !--If activates if the total is LOWER than the max because
506    !--everything is negative
507    IF ( dprecip_evasub_tot .LT. dprecip_evasub_max ) THEN
508      draineva = dprecip_evasub_max * draineva / dprecip_evasub_tot
509      dsnowsub = dprecip_evasub_max * dsnowsub / dprecip_evasub_tot
510    ENDIF
511
512
513    !--New solid and liquid precipitation fluxes after evap and sublimation
514    dqrevap = draineva / dhum_to_dflux(i)
515    dqssubl = dsnowsub / dhum_to_dflux(i)
516
517
518    !--Vapor is updated after evaporation/sublimation (it is increased)
519    qvap(i) = qvap(i) - dqrevap - dqssubl
520    !--qprecip is the total condensed water in the precip flux (it is decreased)
521    qprecip(i) = qprecip(i) + dqrevap + dqssubl
522    !--Air and precip temperature (i.e., gridbox temperature)
523    !--is updated due to latent heat cooling
524    temp(i) = temp(i) &
525            + dqrevap * RLVTT / RCPD &
526            / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) &
527            + dqssubl * RLSTT / RCPD &
528            / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
529
530    !--Add tendencies
531    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
532    rainclr(i) = MAX(0., rainclr(i) + draineva)
533    snowclr(i) = MAX(0., snowclr(i) + dsnowsub)
534
535    !--If there is no more precip fluxes, the precipitation fraction in clear
536    !--sky is set to 0
537    IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0.
538
539    !--Calculation of the total fluxes
540    rain(i) = rainclr(i) + raincld(i)
541    snow(i) = snowclr(i) + snowcld(i)
542
543  ELSEIF ( ( rain(i) + snow(i) ) .LE. 0. ) THEN
544    !--If no precip, we reinitialize the cloud fraction used for the precip to 0
545    precipfraccld(i) = 0.
546    precipfracclr(i) = 0.
547
548  ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. )
549
550  !--Diagnostic tendencies
551  dqssub(i) = dqssubl / dtime
552  dqreva(i) = dqrevap / dtime
553
554ENDDO ! loop on klon
555
556END SUBROUTINE poprecip_precld
557
558
559!----------------------------------------------------------------
560! Computes the processes-oriented precipitation formulations for
561! - autoconversion (auto) via a deposition process
562! - aggregation (agg)
563! - riming (rim)
564! - collection (col)
565! - melting (melt)
566! - freezing (freez)
567!
568SUBROUTINE poprecip_postcld( &
569           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
570           temp, qvap, qliq, qice, icefrac, cldfra, &
571           precipfracclr, precipfraccld, &
572           rain, rainclr, raincld, snow, snowclr, snowcld, &
573           qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
574           dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
575
576USE lmdz_lscp_ini, ONLY : prt_level, lunout
577USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
578USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
579
580USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb,    &
581                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, &
582                          thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
583                          rho_rain, r_rain, r_snow, rho_ice,                   &
584                          tau_auto_snow_min, tau_auto_snow_max,                &
585                          thresh_precip_frac, eps,                             &
586                          gamma_melt, alpha_freez, beta_freez, temp_nowater,   &
587                          iflag_cloudth_vert, iflag_rain_incloud_vol,          &
588                          cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez,       &
589                          rain_fallspeed_clr, rain_fallspeed_cld,              &
590                          snow_fallspeed_clr, snow_fallspeed_cld
591
592
593IMPLICIT NONE
594
595INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
596REAL,    INTENT(IN)                     :: dtime          !--time step [s]
597
598REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
599REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
600REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
601
602REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
603LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
604
605REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
606REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity [kg/kg]
607REAL,    INTENT(INOUT), DIMENSION(klon) :: qliq           !--current liquid water specific humidity [kg/kg]
608REAL,    INTENT(INOUT), DIMENSION(klon) :: qice           !--current ice water specific humidity [kg/kg]
609REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--ice fraction [-]
610REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction [-]
611
612REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
613REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
614                                                          !--NB. at the end of the routine, becomes the fraction of precip
615                                                          !--in the current layer
616
617REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
618REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
619REAL,    INTENT(INOUT), DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
620REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
621REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
622REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
623
624REAL,    INTENT(OUT),   DIMENSION(klon) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
625REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
626REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
627REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
628REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
629REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
630REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsrim         !--snow tendency due to riming [kg/kg/s]
631REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
632REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
633REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
634REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
635
636
637
638!--Local variables
639
640INTEGER :: i
641!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
642REAL, DIMENSION(klon) :: dhum_to_dflux
643REAL, DIMENSION(klon) :: qtot                             !--includes vap, liq, ice and precip
644
645!--Partition of the fluxes
646REAL :: dcldfra
647REAL :: precipfractot
648REAL :: dprecipfracclr, dprecipfraccld
649REAL :: drainclr, dsnowclr
650REAL :: draincld, dsnowcld
651
652!--Collection, aggregation and riming
653REAL :: eff_cldfra
654REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
655REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
656REAL :: rho_snow
657REAL :: dqlcol           !--loss of liquid cloud content due to collection by rain [kg/kg/s]
658REAL :: dqiagg           !--loss of ice cloud content due to collection by aggregation [kg/kg/s]
659REAL :: dqlrim           !--loss of liquid cloud content due to riming on snow [kg/kg/s]
660
661!--Autoconversion
662REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
663REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
664REAL :: dqlauto          !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
665REAL :: dqiauto          !--loss of ice cloud content due to autoconversion to snow [kg/kg/s]
666
667!--Melting
668REAL :: dqsmelt_max, air_thermal_conduct
669REAL :: nb_snowflake_clr, nb_snowflake_cld
670REAL :: capa_snowflake, temp_wetbulb
671REAL :: rho, r_ice
672REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
673REAL, DIMENSION(klon) :: qzero, qsat, dqsat
674
675!--Freezing
676REAL :: dqrfreez_max
677REAL :: tau_freez
678REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2
679REAL :: coef_freez
680REAL :: dqifreez        !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
681REAL :: Eff_rain_ice
682
683
684!--Initialisation of variables
685
686
687qzero(:)    = 0.
688
689dqrcol(:)   = 0.
690dqsagg(:)   = 0.
691dqsauto(:)  = 0.
692dqrauto(:)  = 0.
693dqsrim(:)   = 0.
694dqrmelt(:)  = 0.
695dqsmelt(:)  = 0.
696dqrfreez(:) = 0.
697dqsfreez(:) = 0.
698
699
700DO i = 1, klon
701
702  !--Variables initialisation
703  dqlcol  = 0.
704  dqiagg  = 0.
705  dqiauto = 0.
706  dqlauto = 0.
707  dqlrim  = 0.
708
709  !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
710  dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
711  qtot(i) = qvap(i) + qliq(i) + qice(i) &
712          + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i)
713
714  !------------------------------------------------------------
715  !--     PRECIPITATION FRACTIONS UPDATE
716  !------------------------------------------------------------
717  !--The goal of this routine is to reattribute precipitation fractions
718  !--and fluxes to clear or cloudy air, depending on the variation of
719  !--the cloud fraction on the vertical dimension. We assume a
720  !--maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
721  !--and LTP thesis, 2021)
722  !--NB. in fact, we assume a maximum-random overlap of the total precip. frac
723
724  !--Initialisation
725  precipfractot = precipfracclr(i) + precipfraccld(i)
726
727  !--Instead of using the cloud cover which was use in LTP thesis, we use the
728  !--total precip. fraction to compute the maximum-random overlap. This is
729  !--because all the information of the cloud cover is embedded into
730  !--precipfractot, and this allows for taking into account the potential
731  !--reduction of the precipitation fraction because either the flux is too
732  !--small (see barrier at the end of poprecip_postcld) or the flux is completely
733  !--evaporated (see barrier at the end of poprecip_precld)
734  !--NB. precipfraccld(i) is here the cloud fraction of the layer above
735  !precipfractot = 1. - ( 1. - precipfractot ) * &
736  !               ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
737  !             / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
738
739
740  IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN
741    precipfractot = 1.
742  ELSE
743    precipfractot = 1. - ( 1. - precipfractot ) * &
744                 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
745               / ( 1. - precipfraccld(i) )
746  ENDIF
747
748  !--precipfraccld(i) is here the cloud fraction of the layer above
749  dcldfra = cldfra(i) - precipfraccld(i)
750  !--Tendency of the clear-sky precipitation fraction. We add a MAX on the
751  !--calculation of the current CS precip. frac.
752  !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
753  !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated
754  !--if precipfractot < cldfra)
755  dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i)
756  !--Tendency of the cloudy precipitation fraction. We add a MAX on the
757  !--calculation of the current CS precip. frac.
758  !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) )
759  !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated
760  !--if cldfra < 0)
761  dprecipfraccld = dcldfra
762
763
764  !--If the cloud extends
765  IF ( dprecipfraccld .GT. 0. ) THEN
766    !--If there is no CS precip, nothing happens.
767    !--If there is, we reattribute some of the CS precip flux
768    !--to the cloud precip flux, proportionnally to the
769    !--decrease of the CS precip fraction
770    IF ( precipfracclr(i) .LE. 0. ) THEN
771      drainclr = 0.
772      dsnowclr = 0.
773    ELSE
774      drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i)
775      dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i)
776    ENDIF
777  !--If the cloud narrows
778  ELSEIF ( dprecipfraccld .LT. 0. ) THEN
779    !--We reattribute some of the cloudy precip flux
780    !--to the CS precip flux, proportionnally to the
781    !--decrease of the cloud precip fraction
782    draincld = dprecipfraccld / precipfraccld(i) * raincld(i)
783    dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i)
784    drainclr = - draincld
785    dsnowclr = - dsnowcld
786  !--If the cloud stays the same or if there is no cloud above and
787  !--in the current layer, nothing happens
788  ELSE
789    drainclr = 0.
790    dsnowclr = 0.
791  ENDIF
792
793  !--We add the tendencies
794  !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
795  precipfraccld(i) = precipfraccld(i) + dprecipfraccld
796  precipfracclr(i) = precipfracclr(i) + dprecipfracclr
797  rainclr(i) = MAX(0., rainclr(i) + drainclr)
798  snowclr(i) = MAX(0., snowclr(i) + dsnowclr)
799  raincld(i) = MAX(0., raincld(i) - drainclr)
800  snowcld(i) = MAX(0., snowcld(i) - dsnowclr)
801 
802  !--If vertical heterogeneity is taken into account, we use
803  !--the "true" volume fraction instead of a modified
804  !--surface fraction (which is larger and artificially
805  !--reduces the in-cloud water).
806  IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
807    eff_cldfra = ctot_vol(i)
808  ELSE
809    eff_cldfra = cldfra(i)
810  ENDIF
811
812
813  !--Start precipitation growth processes
814
815  !--If the cloud is big enough, the precipitation processes activate
816  ! TODO met on seuil_neb ici ?
817  IF ( cldfra(i) .GE. seuil_neb ) THEN
818
819    !---------------------------------------------------------
820    !--           COLLECTION AND AGGREGATION
821    !---------------------------------------------------------
822    !--Collection: processus through which rain collects small liquid droplets
823    !--in suspension, and add it to the rain flux
824    !--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
825    !--Those processes are treated before autoconversion because we do not
826    !--want to collect/aggregate the newly formed fluxes, which already
827    !--"saw" the cloud as they come from it
828    !--The formulas come from Muench and Lohmann 2020
829
830    !--gamma_col: tuning coefficient [-]
831    !--rho_rain: volumic mass of rain [kg/m3]
832    !--r_rain: size of the rain droplets [m]
833    !--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
834    !--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
835    !--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
836    !--get in-cloud mean quantities. The two divisions by eff_cldfra are
837    !--then simplified.
838
839    !--The collection efficiency is perfect.
840    Eff_rain_liq = 1.
841    coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
842    IF ( raincld(i) .GT. 0. ) THEN
843      !--Exact explicit version, which does not need a barrier because of
844      !--the exponential decrease
845      dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
846
847      !--Add tendencies
848      qliq(i) = qliq(i) + dqlcol
849      raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i)
850
851      !--Diagnostic tendencies
852      dqrcol(i)  = - dqlcol  / dtime
853    ENDIF
854
855    !--Same as for aggregation
856    !--Eff_snow_liq formula:
857    !--it s a product of a collection efficiency and a sticking efficiency
858    ! Milbrandt and Yau formula that gives very low values:
859    ! Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
860    ! Lin 1983's formula
861    Eff_snow_ice = EXP( 0.025 * MIN( ( temp(i) - RTT ), 0.) )
862    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
863    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
864    coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
865    IF ( snowcld(i) .GT. 0. ) THEN
866      !--Exact explicit version, which does not need a barrier because of
867      !--the exponential decrease
868      dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
869
870      !--Add tendencies
871      qice(i) = qice(i) + dqiagg
872      snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i)
873
874      !--Diagnostic tendencies
875      dqsagg(i)  = - dqiagg  / dtime
876    ENDIF
877
878
879    !---------------------------------------------------------
880    !--                  AUTOCONVERSION
881    !---------------------------------------------------------
882    !--Autoconversion converts liquid droplets/ice crystals into
883    !--rain drops/snowflakes. It relies on the formulations by
884    !--Sundqvist 1978.
885
886    !--If we are in a convective point, we have different parameters
887    !--for the autoconversion
888    IF ( ptconv(i) ) THEN
889        qthresh_auto_rain = cld_lc_con
890        qthresh_auto_snow = cld_lc_con_snow
891
892        tau_auto_rain = cld_tau_con
893        !--tau for snow depends on the ice fraction in mixed-phase clouds     
894        tau_auto_snow = tau_auto_snow_max &
895                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
896
897        expo_auto_rain = cld_expo_con
898        expo_auto_snow = cld_expo_con
899    ELSE
900        qthresh_auto_rain = cld_lc_lsc
901        qthresh_auto_snow = cld_lc_lsc_snow
902
903        tau_auto_rain = cld_tau_lsc
904        !--tau for snow depends on the ice fraction in mixed-phase clouds     
905        tau_auto_snow = tau_auto_snow_max &
906                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
907
908        expo_auto_rain = cld_expo_lsc
909        expo_auto_snow = cld_expo_lsc
910    ENDIF
911
912
913    ! Liquid water quantity to remove according to (Sundqvist, 1978)
914    ! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) )
915    !
916    !--And same formula for ice
917    !
918    !--We first treat the second term (with exponential) in an explicit way
919    !--and then treat the first term (-q/tau) in an exact way
920
921    dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
922               - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
923
924    dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
925               - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
926
927
928    !--Barriers so that we don't create more rain/snow
929    !--than there is liquid/ice
930    dqlauto = MAX( - qliq(i), dqlauto )
931    dqiauto = MAX( - qice(i), dqiauto )
932
933    !--Add tendencies
934    qliq(i) = qliq(i) + dqlauto
935    qice(i) = qice(i) + dqiauto
936    raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
937    snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
938
939    !--Diagnostic tendencies
940    dqsauto(i) = - dqiauto / dtime
941    dqrauto(i) = - dqlauto / dtime
942
943
944    !---------------------------------------------------------
945    !--                  RIMING
946    !---------------------------------------------------------
947    !--Process which converts liquid droplets in suspension into
948    !--snow because of the collision between
949    !--those and falling snowflakes.
950    !--The formula comes from Muench and Lohmann 2020
951    !--NB.: this process needs a temperature adjustment
952
953    !--Eff_snow_liq formula: following Ferrier 1994,
954    !--assuming 1
955    Eff_snow_liq = 1.0
956    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
957    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
958    coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
959    IF ( snowcld(i) .GT. 0. ) THEN
960      !--Exact version, which does not need a barrier because of
961      !--the exponential decrease
962      dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. )
963
964      !--Add tendencies
965      !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
966      qliq(i) = qliq(i) + dqlrim
967      snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i)
968
969      !--Temperature adjustment with the release of latent
970      !--heat because of solid condensation
971      temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
972                        / ( 1. + RVTMP2 * qtot(i) )
973
974      !--Diagnostic tendencies
975      dqsrim(i)  = - dqlrim  / dtime
976    ENDIF
977
978  ENDIF ! cldfra .GE. seuil_neb
979
980ENDDO ! loop on klon
981
982
983!--Re-calculation of saturation specific humidity
984!--because riming changed temperature
985CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
986
987DO i = 1, klon
988
989  !---------------------------------------------------------
990  !--                  MELTING
991  !---------------------------------------------------------
992  !--Process through which snow melts into rain.
993  !--The formula is homemade.
994  !--NB.: this process needs a temperature adjustment
995
996  !--dqsmelt_max         : maximum snow melting so that temperature
997  !--                      stays higher than 273 K [kg/kg]
998  !--capa_snowflake      : capacitance of a snowflake, equal to
999  !--                      the radius if the snowflake is a sphere [m]
1000  !--temp_wetbulb        : wet-bulb temperature [K]
1001  !--snow_fallspeed      : snow fall velocity (in clear/cloudy sky) [m/s]
1002  !--air_thermal_conduct : thermal conductivity of the air [J/m/K/s]
1003  !--gamma_melt          : tuning parameter for melting [-]
1004  !--nb_snowflake        : number of snowflakes (in clear/cloudy air) [-]
1005
1006  IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
1007    !--Computed according to
1008    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1009    dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
1010                        * ( 1. + RVTMP2 * qtot(i) ))
1011   
1012    !--Initialisation
1013    dqsclrmelt = 0.
1014    dqscldmelt = 0.
1015
1016    !--We assume that the snowflakes are spherical
1017    capa_snowflake = r_snow
1018    !--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
1019    air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184   
1020    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
1021    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
1022
1023    !--In clear air
1024    IF ( ( snowclr(i) .GT. 0. ) .AND.  ( precipfracclr(i) .GT. 0. ) ) THEN
1025      !--Formula for the wet-bulb temperature from ECMWF (IFS)
1026      !--The vapor used is the vapor in the clear sky
1027      temp_wetbulb = temp(i) &
1028                   - ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) &
1029                   * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
1030                   - 40.637 * ( temp(i) - 275. ) )
1031      !--Calculated according to
1032      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
1033      nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr &
1034                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
1035      dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
1036                   * capa_snowflake / RLMLT * gamma_melt &
1037                   * MAX(0., temp_wetbulb - RTT) * dtime
1038     
1039      !--Barrier to limit the melting flux to the clr snow flux in the mesh
1040      dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i))
1041    ENDIF
1042
1043
1044    !--In cloudy air
1045    IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
1046      !--As the air is saturated, the wet-bulb temperature is equal to the
1047      !--temperature
1048      temp_wetbulb = temp(i)
1049      !--Calculated according to
1050      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
1051      nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld &
1052                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
1053      dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
1054                   * capa_snowflake / RLMLT * gamma_melt &
1055                   * MAX(0., temp_wetbulb - RTT) * dtime
1056
1057      !--Barrier to limit the melting flux to the cld snow flux in the mesh
1058      dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i))
1059   ENDIF
1060   
1061
1062    !--Barrier on temperature. If the total melting flux leads to a
1063    !--positive temperature, it is limited to keep temperature above 0 degC.
1064    !--It is activated if the total is LOWER than the max
1065    !--because everything is negative
1066    dqstotmelt = dqsclrmelt + dqscldmelt
1067    IF ( dqstotmelt .LT. dqsmelt_max ) THEN
1068      !--We redistribute the max melted snow keeping
1069      !--the clear/cloud partition of the melted snow
1070      dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
1071      dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
1072      dqstotmelt = dqsmelt_max
1073
1074    ENDIF
1075
1076    !--Add tendencies
1077    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1078    rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i))
1079    raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i))
1080    snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i))
1081    snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i))
1082
1083    !--Temperature adjustment with the release of latent
1084    !--heat because of melting
1085    temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
1086                      / ( 1. + RVTMP2 * qtot(i) )
1087
1088    !--Diagnostic tendencies
1089    dqrmelt(i) = - dqstotmelt / dtime
1090    dqsmelt(i) = dqstotmelt / dtime
1091
1092  ENDIF
1093
1094
1095  !---------------------------------------------------------
1096  !--                  FREEZING
1097  !---------------------------------------------------------
1098  !--Process through which rain freezes into snow.
1099  !-- We parameterize it as a 2 step process:
1100  !--first: freezing following collision with ice crystals
1101  !--second: immersion freezing following (inspired by Bigg 1953)
1102  !--the latter is parameterized as an exponential decrease of the rain
1103  !--water content with a homemade formulya
1104  !--This is based on a caracteritic time of freezing, which
1105  !--exponentially depends on temperature so that it is
1106  !--equal to 1 for temp_nowater (see below) and is close to
1107  !--0 for RTT (=273.15 K).
1108  !--NB.: this process needs a temperature adjustment
1109  !--dqrfreez_max : maximum rain freezing so that temperature
1110  !--              stays lower than 273 K [kg/kg]
1111  !--tau_freez    : caracteristic time of freezing [s]
1112  !--gamma_freez  : tuning parameter [s-1]
1113  !--alpha_freez  : tuning parameter for the shape of the exponential curve [-]
1114  !--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC)
1115
1116  IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
1117
1118 
1119    !--1st step: freezing following collision with ice crystals
1120    !--Sub-process of freezing which quantifies the collision between
1121    !--ice crystals in suspension and falling rain droplets.
1122    !--The rain droplets freeze, becoming graupel, and carrying
1123    !--the ice crystal (which acted as an ice nucleating particle).
1124    !--The formula is adapted from the riming formula.
1125    !--it works only in the cloudy part
1126
1127    dqifreez = 0.
1128    dqrtotfreez_step1 = 0.
1129
1130    IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. &
1131         ( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
1132      dqrclrfreez = 0.
1133      dqrcldfreez = 0.
1134
1135      !--Computed according to
1136      !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1137      dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
1138                         * ( 1. + RVTMP2 * qtot(i) ))
1139 
1140
1141      !--The collision efficiency is assumed unity
1142      Eff_rain_ice = 1.
1143      coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
1144      !--Exact version, which does not need a barrier because of
1145      !--the exponential decrease.
1146      dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
1147
1148      !--We add the part of rain water that freezes, limited by a temperature barrier
1149      !--This quantity is calculated assuming that the number of drop that freeze correspond to the number
1150      !--of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
1151      !--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3
1152      !--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3
1153      !--The assumption above corresponds to dNi = dNr, i.e.,
1154      !-- dqr = dqi * (4/3 PI rho_rain * r_rain**3) / (4/3 PI rho_ice * r_ice**3)
1155      !--Dry density [kg/m3]
1156      rho = pplay(i) / temp(i) / RD
1157      !--r_ice formula from Sun and Rikus (1999)
1158      r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 &
1159            + 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
1160      dqrcldfreez = dqifreez * rho_rain * r_rain**3. / ( rho_ice * r_ice**3. )
1161      dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i))
1162      dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max)
1163      dqrtotfreez_step1 = dqrcldfreez
1164
1165      !--Add tendencies
1166      !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1167      qice(i) = qice(i) + dqifreez
1168      raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
1169      snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i))
1170      temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD &
1171                        / ( 1. + RVTMP2 * qtot(i) )
1172
1173    ENDIF
1174   
1175    !-- Second step immersion freezing of rain drops
1176    !-- with a homemade timeconstant depending on temperature
1177
1178    dqrclrfreez = 0.
1179    dqrcldfreez = 0.
1180    dqrtotfreez_step2 = 0.
1181    !--Computed according to
1182    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1183
1184    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
1185                         * ( 1. + RVTMP2 * qtot(i) ))
1186
1187   
1188    tau_freez = 1. / ( beta_freez &
1189              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
1190
1191
1192    !--In clear air
1193    IF ( rainclr(i) .GT. 0. ) THEN
1194      !--Exact solution of dqrain/dt = -qrain/tau_freez
1195      dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
1196    ENDIF
1197
1198    !--In cloudy air
1199    IF ( raincld(i) .GT. 0. ) THEN
1200      !--Exact solution of dqrain/dt = -qrain/tau_freez
1201      dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
1202    ENDIF
1203
1204    !--temperature barrier step 2
1205    !--It is activated if the total is LOWER than the max
1206    !--because everything is negative
1207    dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
1208 
1209    IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
1210      !--We redistribute the max freezed rain keeping
1211      !--the clear/cloud partition of the freezing rain
1212      dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
1213      dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
1214      dqrtotfreez_step2 = dqrfreez_max
1215    ENDIF
1216
1217
1218    !--Add tendencies
1219    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1220    rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i))
1221    raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
1222    snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i))
1223    snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i))
1224
1225
1226    !--Temperature adjustment with the uptake of latent
1227    !--heat because of freezing
1228    temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
1229                      / ( 1. + RVTMP2 * qtot(i) )
1230
1231    !--Diagnostic tendencies
1232    dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2         
1233    dqrfreez(i) = dqrtotfreez / dtime
1234    dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
1235
1236  ENDIF
1237
1238
1239
1240  !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
1241  !--we reduce the precipiration fraction in the clear/cloudy air so that the new
1242  !--local flux of rain+snow is equal to rain_int_min.
1243  !--Here, rain+snow is the gridbox-mean flux of precip.
1244  !--Therefore, (rain+snow)/precipfrac is the local flux of precip.
1245  !--If the local flux of precip is lower than rain_int_min, i.e.,
1246  !-- (rain+snow)/precipfrac < rain_int_min , i.e.,
1247  !-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
1248  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
1249  !--Note that this is physically different than what is proposed in LTP thesis.
1250  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
1251  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
1252
1253  !--Calculate outputs
1254  rain(i) = rainclr(i) + raincld(i)
1255  snow(i) = snowclr(i) + snowcld(i)
1256
1257  !--Diagnostics
1258  !--BEWARE this is indeed a diagnostic: this is an estimation from
1259  !--the value of the flux at the bottom interface of the mesh and
1260  !--and assuming an upstream numerical calculation
1261  !--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is
1262  !--used for computing the total ice water content in the mesh
1263  !--for radiation only
1264  qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr &
1265                 + raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld )
1266  qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr &
1267                 + snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld )
1268
1269
1270ENDDO ! loop on klon
1271
1272END SUBROUTINE poprecip_postcld
1273
1274END MODULE lmdz_lscp_precip
Note: See TracBrowser for help on using the repository browser.