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

Last change on this file since 5426 was 5413, checked in by evignon, 40 hours ago

fin de l'externalisation des precips dans lscp.

  1. Borella, E Vignon
File size: 69.3 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, znebprecip, 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) :: znebprecip     !--fraction of precipitation IN THE LAYER ABOVE [-]
49REAL,    INTENT(INOUT), DIMENSION(klon) :: znebprecipclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
50
51REAL,    INTENT(INOUT), DIMENSION(klon) :: zrfl           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
52REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
53REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflcld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
54REAL,    INTENT(INOUT), DIMENSION(klon) :: zifl           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
55REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
56REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
57
58REAL,    INTENT(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
59REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
60
61
62REAL :: zmair
63REAL :: zcpair, zcpeau
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      znebprecip(i)=0.
298
299    ENDIF ! (zrfl(i)+zifl(i).GT.0.)
300
301  ENDDO ! loop on klon
302
303ENDIF ! (iflag_evap_prec>=1)
304
305END SUBROUTINE histprecip_precld
306
307
308SUBROUTINE histprecip_postcld( &
309        klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, zdqsdT_raw, &
310        zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
311        rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, &
312        zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
313        zradocond, zradoice, dqrauto, dqsauto &
314        )
315
316USE lmdz_lscp_ini, ONLY : RD, RG, RCPD, RVTMP2, RLSTT, RLMLT, RTT
317USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, ffallv_con, &
318                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, ffallv_lsc, &
319                          seuil_neb, rain_int_min, cice_velo, dice_velo
320USE lmdz_lscp_ini, ONLY : iflag_evap_prec, iflag_cloudth_vert, iflag_rain_incloud_vol, &
321                          iflag_autoconversion, ok_radocond_snow, ok_bug_phase_lscp, &
322                          niter_lscp
323USE lmdz_lscp_tools, ONLY : fallice_velocity
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
332REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
333REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
334REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
335REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
336LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
337REAL,    INTENT(IN),    DIMENSION(klon) :: zdqsdT_raw     !--derivative of qsat w.r.t. temperature times L/Cp [SI]
338
339REAL,    INTENT(INOUT), DIMENSION(klon) :: zt             !--current temperature [K]
340REAL,    INTENT(INOUT), DIMENSION(klon) :: zq             !--current water vapor specific humidity [kg/kg]
341REAL,    INTENT(INOUT), DIMENSION(klon) :: zoliq          !--current liquid+ice water specific humidity [kg/kg]
342REAL,    INTENT(INOUT), DIMENSION(klon) :: zoliql         !--current liquid water specific humidity [kg/kg]
343REAL,    INTENT(INOUT), DIMENSION(klon) :: zoliqi         !--current ice water specific humidity [kg/kg]
344REAL,    INTENT(INOUT), DIMENSION(klon) :: zcond          !--liquid+ice water specific humidity AFTER CONDENSATION (no precip process) [kg/kg]
345REAL,    INTENT(IN),    DIMENSION(klon) :: zfice          !--ice fraction AFTER CONDENSATION [-]
346REAL,    INTENT(IN),    DIMENSION(klon) :: zmqc           !--specific humidity in the precipitation falling from the upper layer [kg/kg]
347
348REAL,    INTENT(IN),    DIMENSION(klon) :: rneb           !--cloud fraction [-]
349REAL,    INTENT(INOUT), DIMENSION(klon) :: znebprecipclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
350REAL,    INTENT(INOUT), DIMENSION(klon) :: znebprecipcld  !--fraction of precipitation in the cloud IN THE LAYER ABOVE [-]
351REAL,    INTENT(INOUT), DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
352REAL,    INTENT(INOUT), DIMENSION(klon) :: tot_zneb       !--total cloud cover above the current mesh [-]
353REAL,    INTENT(INOUT), DIMENSION(klon) :: zrho_up        !--air density IN THE LAYER ABOVE [kg/m3]
354REAL,    INTENT(INOUT), DIMENSION(klon) :: zvelo_up       !--ice fallspeed velocity IN THE LAYER ABOVE [m/s]
355
356REAL,    INTENT(INOUT), DIMENSION(klon) :: zrfl           !--flux of rain gridbox-mean [kg/s/m2]
357REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflclr        !--flux of rain gridbox-mean in clear sky [kg/s/m2]
358REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflcld        !--flux of rain gridbox-mean in cloudy air [kg/s/m2]
359REAL,    INTENT(INOUT), DIMENSION(klon) :: zifl           !--flux of snow gridbox-mean [kg/s/m2]
360REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky [kg/s/m2]
361REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
362
363REAL,    INTENT(OUT),   DIMENSION(klon) :: zradocond      !--condensed water used in the radiation scheme [kg/kg]
364REAL,    INTENT(OUT),   DIMENSION(klon) :: zradoice       !--condensed ice water used in the radiation scheme [kg/kg]
365REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
366REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
367
368
369! Local variables for precip fraction update
370REAL :: smallestreal
371REAL, DIMENSION(klon) :: tot_znebn, d_tot_zneb
372REAL, DIMENSION(klon) :: d_znebprecip_cld_clr, d_znebprecip_clr_cld
373REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
374REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
375
376! Local variables for autoconversion
377REAL :: zct, zcl, zexpo, ffallv
378REAL :: zchau, zfroi
379REAL :: zrain, zsnow, zprecip
380REAL :: effective_zneb
381REAL, DIMENSION(klon) :: zrho, zvelo
382REAL, DIMENSION(klon) :: zdz, iwc
383
384! Local variables for Bergeron process
385REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
386REAL, DIMENSION(klon) :: zqpreci, zqprecl
387
388! Local variables for calculation of quantity of condensates seen by the radiative scheme
389REAL, DIMENSION(klon) :: zradoliq
390REAL, DIMENSION(klon) :: ziflprev
391
392! Misc
393INTEGER :: i, n
394
395! Initialisation
396smallestreal=1.e-9
397zqprecl(:) = 0.
398zqpreci(:) = 0.
399ziflprev(:) = 0.
400
401
402IF (iflag_evap_prec .GE. 4) THEN
403
404  !Partitionning between precipitation coming from clouds and that coming from CS
405
406  !0) Calculate tot_zneb, total cloud fraction above the cloud
407  !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
408 
409  DO i=1, klon
410    tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i),zneb(i))) &
411            /(1.-min(zneb(i),1.-smallestreal))
412    d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
413    tot_zneb(i) = tot_znebn(i)
414
415
416    !1) Cloudy to clear air
417    d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i),znebprecipcld(i))
418    IF (znebprecipcld(i) .GT. 0.) THEN
419      d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
420      d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
421    ELSE
422      d_zrfl_cld_clr(i) = 0.
423      d_zifl_cld_clr(i) = 0.
424    ENDIF
425
426    !2) Clear to cloudy air
427    d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i) - d_tot_zneb(i) - zneb(i)))
428    IF (znebprecipclr(i) .GT. 0) THEN
429      d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
430      d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
431    ELSE
432      d_zrfl_clr_cld(i) = 0.
433      d_zifl_clr_cld(i) = 0.
434    ENDIF
435
436    !Update variables
437    znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
438    znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
439    zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
440    ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
441    zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
442    ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
443
444    ! c_iso  do the same thing for isotopes precip
445  ENDDO
446ENDIF
447
448
449! Autoconversion
450! -------------------------------------------------------------------------------
451
452! Initialisation
453DO i = 1, klon
454  zrho(i)  = pplay(i) / zt(i) / RD
455  zdz(i)   = (paprsdn(i)-paprsup(i)) / (zrho(i)*RG)
456  iwc(i)   = 0.
457  zneb(i)  = MAX(rneb(i),seuil_neb)
458
459  ! variables for calculation of quantity of condensates seen by the radiative scheme
460  ! NB. only zradocond and zradoice are outputed, but zradoliq is used if ok_radocond_snow
461  ! is TRUE
462  zradocond(i) = zoliq(i)/REAL(niter_lscp+1)
463  zradoliq(i)  = zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
464  zradoice(i)  = zoliq(i)*zfice(i)/REAL(niter_lscp+1)
465ENDDO
466
467   
468DO n = 1, niter_lscp
469
470  DO i=1,klon
471    IF (rneb(i).GT.0.0) THEN
472      iwc(i) = zrho(i) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
473    ENDIF
474  ENDDO
475
476  CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, zvelo)
477
478  DO i = 1, klon
479
480    IF (rneb(i).GT.0.0) THEN
481
482      ! Initialization of zrain, zsnow and zprecip:
483      zrain=0.
484      zsnow=0.
485      zprecip=0.
486      ! c_iso same init for isotopes. Externalisation?
487
488      IF (zneb(i).EQ.seuil_neb) THEN
489        zprecip = 0.0
490        zsnow = 0.0
491        zrain= 0.0
492      ELSE
493
494        IF (ptconv(i)) THEN ! if convective point
495          zcl=cld_lc_con
496          zct=1./cld_tau_con
497          zexpo=cld_expo_con
498          ffallv=ffallv_con
499        ELSE
500          zcl=cld_lc_lsc
501          zct=1./cld_tau_lsc
502          zexpo=cld_expo_lsc
503          ffallv=ffallv_lsc
504        ENDIF
505
506
507        ! if vertical heterogeneity is taken into account, we use
508        ! the "true" volume fraction instead of a modified
509        ! surface fraction (which is larger and artificially
510        ! reduces the in-cloud water).
511        IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
512          effective_zneb=ctot_vol(i)
513        ELSE
514          effective_zneb=zneb(i)
515        ENDIF
516
517
518        ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
519        ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
520        !.........................................................
521        IF (iflag_autoconversion .EQ. 2) THEN
522        ! two-steps resolution with niter_lscp=1 sufficient
523        ! we first treat the second term (with exponential) in an explicit way
524        ! and then treat the first term (-q/tau) in an exact way
525          zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct &
526            *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo))))
527        ELSE
528        ! old explicit resolution with subtimesteps
529          zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) &
530            *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo))
531        ENDIF
532
533
534        ! Ice water quantity to remove (Zender & Kiehl, 1997)
535        ! dqice/dt=1/rho*d(rho*wice*qice)/dz
536        !....................................
537        IF (iflag_autoconversion .EQ. 2) THEN
538        ! exact resolution, niter_lscp=1 is sufficient but works only
539        ! with iflag_vice=0
540          IF (zoliqi(i) .GT. 0.) THEN
541            zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
542              +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
543          ELSE
544            zfroi=0.
545          ENDIF
546        ELSE
547        ! old explicit resolution with subtimesteps
548          zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*zvelo(i)
549        ENDIF
550
551        zrain   = MIN(MAX(zchau,0.0),zoliql(i))
552        zsnow   = MIN(MAX(zfroi,0.0),zoliqi(i))
553        zprecip = MAX(zrain + zsnow,0.0)
554
555      ENDIF
556
557
558      IF (iflag_autoconversion .GE. 1) THEN
559        ! debugged version with phase conservation through the autoconversion process
560        zoliql(i) = MAX(zoliql(i)-1.*zrain  , 0.0)
561        zoliqi(i) = MAX(zoliqi(i)-1.*zsnow  , 0.0)
562        zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
563      ELSE
564        ! bugged version with phase resetting
565        zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
566        zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
567        zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
568      ENDIF
569
570      ! c_iso: call isotope_conversion (for readibility) or duplicate
571
572      ! variables for calculation of quantity of condensates seen by the radiative scheme
573      zradocond(i) = zradocond(i) + zoliq(i)/REAL(niter_lscp+1)
574      zradoliq(i)  = zradoliq(i) + zoliql(i)/REAL(niter_lscp+1)
575      zradoice(i)  = zradoice(i) + zoliqi(i)/REAL(niter_lscp+1)
576
577      !--Diagnostics
578      dqrauto(i) = dqrauto(i) + zrain / dtime
579      dqsauto(i) = dqsauto(i) + zsnow / dtime
580
581    ENDIF ! rneb >0
582
583  ENDDO  ! i = 1,klon
584
585ENDDO      ! n = 1,niter
586
587! Precipitation flux calculation
588
589DO i = 1, klon
590       
591  IF (iflag_evap_prec.GE.4) THEN
592    ziflprev(i)=ziflcld(i)
593  ELSE
594    ziflprev(i)=zifl(i)*zneb(i)
595  ENDIF
596
597  IF (rneb(i) .GT. 0.0) THEN
598
599    ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
600    ! If T<0C, liquid precip are converted into ice, which leads to an increase in
601    ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
602    ! taken into account through a linearization of the equations and by approximating
603    ! the liquid precipitation process with a "threshold" process. We assume that
604    ! condensates are not modified during this operation. Liquid precipitation is
605    ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
606    ! Water vapor increases as well
607    ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
608
609    IF (ok_bug_phase_lscp) THEN
610      zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
611      zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
612    ELSE
613      zqpreci(i)=zcond(i)*zfice(i)-zoliqi(i)
614      zqprecl(i)=zcond(i)*(1.-zfice(i))-zoliql(i)
615    ENDIF
616    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
617    coef1 = rneb(i)*RLSTT/zcp*zdqsdT_raw(i)
618    ! Computation of DT if all the liquid precip freezes
619    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
620    ! T should not exceed the freezing point
621    ! that is Delta > RTT-zt(i)
622    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
623    zt(i)  = zt(i) + DeltaT
624    ! water vaporization due to temp. increase
625    Deltaq = rneb(i)*zdqsdT_raw(i)*DeltaT
626    ! we add this vaporized water to the total vapor and we remove it from the precip
627    zq(i) = zq(i) +  Deltaq
628    ! The three "max" lines herebelow protect from rounding errors
629    zcond(i) = max( zcond(i)- Deltaq, 0. )
630    ! liquid precipitation converted to ice precip
631    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
632    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
633    ! iced water budget
634    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
635
636    ! c_iso : mv here condensation of isotopes + redispatchage en precip
637
638    IF (iflag_evap_prec.GE.4) THEN
639      zrflcld(i) = zrflcld(i)+zqprecl(i) &
640        *(paprsdn(i)-paprsup(i))/(RG*dtime)
641      ziflcld(i) = ziflcld(i)+ zqpreci(i) &
642        *(paprsdn(i)-paprsup(i))/(RG*dtime)
643      znebprecipcld(i) = rneb(i)
644      zrfl(i) = zrflcld(i) + zrflclr(i)
645      zifl(i) = ziflcld(i) + ziflclr(i)
646    ELSE
647      zrfl(i) = zrfl(i)+ zqprecl(i)        &
648        *(paprsdn(i)-paprsup(i))/(RG*dtime)
649      zifl(i) = zifl(i)+ zqpreci(i)        &
650        *(paprsdn(i)-paprsup(i))/(RG*dtime)
651    ENDIF
652    ! c_iso : same for isotopes
653
654  ENDIF ! rneb>0
655
656ENDDO
657
658! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
659! if iflag_evap_prec>=4
660IF (iflag_evap_prec.GE.4) THEN
661
662  DO i=1,klon
663
664    IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
665      znebprecipclr(i) = min(znebprecipclr(i),max( &
666        zrflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), &
667        ziflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
668    ELSE
669      znebprecipclr(i)=0.0
670    ENDIF
671
672    IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
673      znebprecipcld(i) = min(znebprecipcld(i), max( &
674        zrflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), &
675        ziflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
676    ELSE
677      znebprecipcld(i)=0.0
678    ENDIF
679  ENDDO
680
681ENDIF
682
683! we recalculate zradoice to account for contributions from the precipitation flux
684! if ok_radocond_snow is true
685IF ( ok_radocond_snow ) THEN
686  IF ( .NOT. iftop ) THEN
687    ! for the solid phase (crystals + snowflakes)
688    ! we recalculate zradoice by summing
689    ! the ice content calculated in the mesh
690    ! + the contribution of the non-evaporated snowfall
691    ! from the overlying layer
692    DO i=1,klon
693      IF (ziflprev(i) .NE. 0.0) THEN
694        zradoice(i)=zoliqi(i)+zqpreci(i)+ziflprev(i)/zrho_up(i)/zvelo_up(i)
695      ELSE
696        zradoice(i)=zoliqi(i)+zqpreci(i)
697      ENDIF
698      zradocond(i)=zradoliq(i)+zradoice(i)
699    ENDDO
700  ENDIF
701  ! keep in memory air density and ice fallspeed velocity
702  zrho_up(:) = zrho(:)
703  zvelo_up(:) = zvelo(:)
704ENDIF
705
706END SUBROUTINE histprecip_postcld
707
708
709!----------------------------------------------------------------
710! Computes the processes-oriented precipitation formulations for
711! evaporation and sublimation
712!
713SUBROUTINE poprecip_precld( &
714           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
715           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
716           rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub &
717           )
718
719USE lmdz_lscp_ini, ONLY : prt_level, lunout
720USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
721USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
722USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub
723USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
724
725IMPLICIT NONE
726
727
728INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
729REAL,    INTENT(IN)                     :: dtime          !--time step [s]
730LOGICAL, INTENT(IN)                     :: iftop          !--if top of the column
731
732
733REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
734REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
735REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
736
737REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
738REAL,    INTENT(IN),    DIMENSION(klon) :: tempupnew      !--updated temperature of the overlying layer [K]
739
740REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
741REAL,    INTENT(INOUT), DIMENSION(klon) :: qprecip        !--specific humidity in the precipitation falling from the upper layer [kg/kg]
742
743REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
744REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
745
746REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
747REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
748
749REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
750REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
751REAL,    INTENT(IN),    DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
752REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
753REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
754REAL,    INTENT(IN),    DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
755
756REAL,    INTENT(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
757REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
758
759
760!--Integer for interating over klon
761INTEGER :: i
762!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
763REAL, DIMENSION(klon) :: dhum_to_dflux
764!--
765REAL, DIMENSION(klon) :: rho, dz
766
767!--Saturation values
768REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
769!--Vapor in the clear sky
770REAL :: qvapclr
771!--Fluxes tendencies because of evaporation and sublimation
772REAL :: dprecip_evasub_max, draineva, dsnowsub, dprecip_evasub_tot
773!--Specific humidity tendencies because of evaporation and sublimation
774REAL :: dqrevap, dqssubl
775!--Specific heat constant
776REAL :: cpair, cpw
777
778!--Initialisation
779qzero(:)  = 0.
780dqreva(:) = 0.
781dqssub(:) = 0.
782dqrevap   = 0.
783dqssubl   = 0.
784
785!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
786dhum_to_dflux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime
787rho(:) = pplay(:) / temp(:) / RD
788dz(:) = ( paprsdn(:) - paprsup(:) ) / RG / rho(:)
789
790!--Calculation of saturation specific humidity
791!--depending on temperature:
792CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:))
793!--wrt liquid water
794CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
795!--wrt ice
796CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
797
798
799
800!--First step consists in "thermalizing" the layer:
801!--as the flux of precip from layer above "advects" some heat (as the precip is at the temperature
802!--of the overlying layer) we recalculate a mean temperature that both the air and the precip in the
803!--layer have.
804
805IF (iftop) THEN
806
807  DO i = 1, klon
808    qprecip(i) = 0.
809  ENDDO
810
811ELSE
812
813  DO i = 1, klon
814    !--No condensed water so cp=cp(vapor+dry air)
815    !-- RVTMP2=rcpv/rcpd-1
816    cpair = RCPD * ( 1. + RVTMP2 * qvap(i) )
817    cpw = RCPD * RVTMP2
818    !--qprecip has to be thermalized with
819    !--layer's air so that precipitation at the ground has the
820    !--same temperature as the lowermost layer
821    !--we convert the flux into a specific quantity qprecip
822    qprecip(i) = ( rain(i) + snow(i) ) / dhum_to_dflux(i)
823    !-- t(i,k+1) + d_t(i,k+1): new temperature of the overlying layer
824    temp(i) = ( tempupnew(i) * qprecip(i) * cpw + cpair * temp(i) ) &
825            / ( cpair + qprecip(i) * cpw )
826  ENDDO
827
828ENDIF
829
830! TODO Probleme : on utilise qvap total dans la maille pour l'evap / sub
831! alors qu'on n'evap / sub que dans le ciel clair
832! deux options pour cette routine :
833! - soit on diagnostique le nuage AVANT l'evap / sub et on estime donc
834! la fraction precipitante ciel clair dans la maille, ce qui permet de travailler
835! avec des fractions, des fluxs et surtout un qvap dans le ciel clair
836! - soit on pousse la param de Ludo au bout, et on prend un qvap de k+1
837! dans le ciel clair, avec un truc comme :
838!   qvapclr(k) = qvapclr(k+1)/qtot(k+1) * qtot(k)
839! UPDATE : on code la seconde version. A voir si on veut mettre la premiere version.
840
841
842DO i = 1, klon
843
844  !--If there is precipitation from the layer above
845  ! NOTE TODO here we could replace the condition on precipfracclr(i) by a condition
846  ! such as eps or thresh_precip_frac, to remove the senseless barrier in the formulas
847  ! of evap / sublim
848  IF ( ( ( rain(i) + snow(i) ) .GT. 0. ) .AND. ( precipfracclr(i) .GT. 0. ) ) THEN
849
850    IF ( ok_corr_vap_evasub ) THEN
851      !--Corrected version - we use the same water ratio between
852      !--the clear and the cloudy sky as in the layer above. This
853      !--extends the assumption that the cloud fraction is the same
854      !--as the layer above. This is assumed only for the evap / subl
855      !--process
856      !--Note that qvap(i) is the total water in the gridbox, and
857      !--precipfraccld(i) is the cloud fraction in the layer above
858      qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
859    ELSE
860      !--Legacy version from Ludo - we use the total specific humidity
861      !--for the evap / subl process
862      qvapclr = qvap(i)
863    ENDIF
864
865    !--Evaporation of liquid precipitation coming from above
866    !--in the clear sky only
867    !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva)
868    !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
869    !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly)
870    !--which does not need a barrier on rainclr, because included in the formula
871    draineva = precipfracclr(i) * ( MAX(0., &
872             - coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) &
873             + ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_eva ) &
874               ) )**( 1. / ( 1. - expo_eva ) ) - rainclr(i)
875             
876    !--Evaporation is limited by 0
877    draineva = MIN(0., draineva)
878
879
880    !--Sublimation of the solid precipitation coming from above
881    !--(same formula as for liquid precip)
882    !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly)
883    !--which does not need a barrier on snowclr, because included in the formula
884    dsnowsub = precipfracclr(i) * ( MAX(0., &
885             - coef_sub * ( 1. - expo_sub ) * (1. - qvapclr / qsati(i)) * dz(i) &
886             + ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_sub ) &
887             ) )**( 1. / ( 1. - expo_sub ) ) - snowclr(i)
888
889    !--Sublimation is limited by 0
890    ! TODO: change max when we will allow for vapor deposition in supersaturated regions
891    dsnowsub = MIN(0., dsnowsub)
892
893    !--Evaporation limit: we ensure that the layer's fraction below
894    !--the clear sky does not reach saturation. In this case, we
895    !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
896    !--Max evaporation is computed not to saturate the clear sky precip fraction
897    !--(i.e., the fraction where evaporation occurs)
898    !--It is expressed as a max flux dprecip_evasub_max
899   
900    dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr(i)) &
901                     * dhum_to_dflux(i)
902    dprecip_evasub_tot = draineva + dsnowsub
903
904    !--Barriers
905    !--If activates if the total is LOWER than the max because
906    !--everything is negative
907    IF ( dprecip_evasub_tot .LT. dprecip_evasub_max ) THEN
908      draineva = dprecip_evasub_max * draineva / dprecip_evasub_tot
909      dsnowsub = dprecip_evasub_max * dsnowsub / dprecip_evasub_tot
910    ENDIF
911
912
913    !--New solid and liquid precipitation fluxes after evap and sublimation
914    dqrevap = draineva / dhum_to_dflux(i)
915    dqssubl = dsnowsub / dhum_to_dflux(i)
916
917
918    !--Vapor is updated after evaporation/sublimation (it is increased)
919    qvap(i) = qvap(i) - dqrevap - dqssubl
920    !--qprecip is the total condensed water in the precip flux (it is decreased)
921    qprecip(i) = qprecip(i) + dqrevap + dqssubl
922    !--Air and precip temperature (i.e., gridbox temperature)
923    !--is updated due to latent heat cooling
924    temp(i) = temp(i) &
925            + dqrevap * RLVTT / RCPD &
926            / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) &
927            + dqssubl * RLSTT / RCPD &
928            / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
929
930    !--Add tendencies
931    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
932    rainclr(i) = MAX(0., rainclr(i) + draineva)
933    snowclr(i) = MAX(0., snowclr(i) + dsnowsub)
934
935    !--If there is no more precip fluxes, the precipitation fraction in clear
936    !--sky is set to 0
937    IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0.
938
939    !--Calculation of the total fluxes
940    rain(i) = rainclr(i) + raincld(i)
941    snow(i) = snowclr(i) + snowcld(i)
942
943  ELSEIF ( ( rain(i) + snow(i) ) .LE. 0. ) THEN
944    !--If no precip, we reinitialize the cloud fraction used for the precip to 0
945    precipfraccld(i) = 0.
946    precipfracclr(i) = 0.
947
948  ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. )
949
950  !--Diagnostic tendencies
951  dqssub(i) = dqssubl / dtime
952  dqreva(i) = dqrevap / dtime
953
954ENDDO ! loop on klon
955
956END SUBROUTINE poprecip_precld
957
958
959!----------------------------------------------------------------
960! Computes the processes-oriented precipitation formulations for
961! - autoconversion (auto) via a deposition process
962! - aggregation (agg)
963! - riming (rim)
964! - collection (col)
965! - melting (melt)
966! - freezing (freez)
967!
968SUBROUTINE poprecip_postcld( &
969           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
970           temp, qvap, qliq, qice, icefrac, cldfra, &
971           precipfracclr, precipfraccld, &
972           rain, rainclr, raincld, snow, snowclr, snowcld, &
973           qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
974           dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
975
976USE lmdz_lscp_ini, ONLY : prt_level, lunout
977USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
978USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
979
980USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb,    &
981                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, &
982                          thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
983                          rho_rain, r_rain, r_snow, rho_ice,                   &
984                          tau_auto_snow_min, tau_auto_snow_max,                &
985                          thresh_precip_frac, eps,                             &
986                          gamma_melt, alpha_freez, beta_freez, temp_nowater,   &
987                          iflag_cloudth_vert, iflag_rain_incloud_vol,          &
988                          cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez,       &
989                          rain_fallspeed_clr, rain_fallspeed_cld,              &
990                          snow_fallspeed_clr, snow_fallspeed_cld
991
992
993IMPLICIT NONE
994
995INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
996REAL,    INTENT(IN)                     :: dtime          !--time step [s]
997
998REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
999REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
1000REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
1001
1002REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
1003LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
1004
1005REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
1006REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity [kg/kg]
1007REAL,    INTENT(INOUT), DIMENSION(klon) :: qliq           !--current liquid water specific humidity [kg/kg]
1008REAL,    INTENT(INOUT), DIMENSION(klon) :: qice           !--current ice water specific humidity [kg/kg]
1009REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--ice fraction [-]
1010REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction [-]
1011
1012REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
1013REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
1014                                                          !--NB. at the end of the routine, becomes the fraction of precip
1015                                                          !--in the current layer
1016
1017REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
1018REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
1019REAL,    INTENT(INOUT), DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
1020REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
1021REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
1022REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
1023
1024REAL,    INTENT(OUT),   DIMENSION(klon) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
1025REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
1026REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
1027REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
1028REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
1029REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
1030REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsrim         !--snow tendency due to riming [kg/kg/s]
1031REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
1032REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
1033REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
1034REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
1035
1036
1037
1038!--Local variables
1039
1040INTEGER :: i
1041!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
1042REAL, DIMENSION(klon) :: dhum_to_dflux
1043REAL, DIMENSION(klon) :: qtot                             !--includes vap, liq, ice and precip
1044
1045!--Partition of the fluxes
1046REAL :: dcldfra
1047REAL :: precipfractot
1048REAL :: dprecipfracclr, dprecipfraccld
1049REAL :: drainclr, dsnowclr
1050REAL :: draincld, dsnowcld
1051
1052!--Collection, aggregation and riming
1053REAL :: eff_cldfra
1054REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
1055REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
1056REAL :: rho_snow
1057REAL :: dqlcol           !--loss of liquid cloud content due to collection by rain [kg/kg/s]
1058REAL :: dqiagg           !--loss of ice cloud content due to collection by aggregation [kg/kg/s]
1059REAL :: dqlrim           !--loss of liquid cloud content due to riming on snow [kg/kg/s]
1060
1061!--Autoconversion
1062REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
1063REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
1064REAL :: dqlauto          !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
1065REAL :: dqiauto          !--loss of ice cloud content due to autoconversion to snow [kg/kg/s]
1066
1067!--Melting
1068REAL :: dqsmelt_max, air_thermal_conduct
1069REAL :: nb_snowflake_clr, nb_snowflake_cld
1070REAL :: capa_snowflake, temp_wetbulb
1071REAL :: rho, r_ice
1072REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
1073REAL, DIMENSION(klon) :: qzero, qsat, dqsat
1074
1075!--Freezing
1076REAL :: dqrfreez_max
1077REAL :: tau_freez
1078REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2
1079REAL :: coef_freez
1080REAL :: dqifreez        !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
1081REAL :: Eff_rain_ice
1082
1083
1084!--Initialisation of variables
1085
1086
1087qzero(:)    = 0.
1088
1089dqrcol(:)   = 0.
1090dqsagg(:)   = 0.
1091dqsauto(:)  = 0.
1092dqrauto(:)  = 0.
1093dqsrim(:)   = 0.
1094dqrmelt(:)  = 0.
1095dqsmelt(:)  = 0.
1096dqrfreez(:) = 0.
1097dqsfreez(:) = 0.
1098
1099
1100DO i = 1, klon
1101
1102  !--Variables initialisation
1103  dqlcol  = 0.
1104  dqiagg  = 0.
1105  dqiauto = 0.
1106  dqlauto = 0.
1107  dqlrim  = 0.
1108
1109  !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
1110  dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
1111  qtot(i) = qvap(i) + qliq(i) + qice(i) &
1112          + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i)
1113
1114  !------------------------------------------------------------
1115  !--     PRECIPITATION FRACTIONS UPDATE
1116  !------------------------------------------------------------
1117  !--The goal of this routine is to reattribute precipitation fractions
1118  !--and fluxes to clear or cloudy air, depending on the variation of
1119  !--the cloud fraction on the vertical dimension. We assume a
1120  !--maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
1121  !--and LTP thesis, 2021)
1122  !--NB. in fact, we assume a maximum-random overlap of the total precip. frac
1123
1124  !--Initialisation
1125  precipfractot = precipfracclr(i) + precipfraccld(i)
1126
1127  !--Instead of using the cloud cover which was use in LTP thesis, we use the
1128  !--total precip. fraction to compute the maximum-random overlap. This is
1129  !--because all the information of the cloud cover is embedded into
1130  !--precipfractot, and this allows for taking into account the potential
1131  !--reduction of the precipitation fraction because either the flux is too
1132  !--small (see barrier at the end of poprecip_postcld) or the flux is completely
1133  !--evaporated (see barrier at the end of poprecip_precld)
1134  !--NB. precipfraccld(i) is here the cloud fraction of the layer above
1135  !precipfractot = 1. - ( 1. - precipfractot ) * &
1136  !               ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
1137  !             / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
1138
1139
1140  IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN
1141    precipfractot = 1.
1142  ELSE
1143    precipfractot = 1. - ( 1. - precipfractot ) * &
1144                 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
1145               / ( 1. - precipfraccld(i) )
1146  ENDIF
1147
1148  !--precipfraccld(i) is here the cloud fraction of the layer above
1149  dcldfra = cldfra(i) - precipfraccld(i)
1150  !--Tendency of the clear-sky precipitation fraction. We add a MAX on the
1151  !--calculation of the current CS precip. frac.
1152  !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
1153  !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated
1154  !--if precipfractot < cldfra)
1155  dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i)
1156  !--Tendency of the cloudy precipitation fraction. We add a MAX on the
1157  !--calculation of the current CS precip. frac.
1158  !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) )
1159  !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated
1160  !--if cldfra < 0)
1161  dprecipfraccld = dcldfra
1162
1163
1164  !--If the cloud extends
1165  IF ( dprecipfraccld .GT. 0. ) THEN
1166    !--If there is no CS precip, nothing happens.
1167    !--If there is, we reattribute some of the CS precip flux
1168    !--to the cloud precip flux, proportionnally to the
1169    !--decrease of the CS precip fraction
1170    IF ( precipfracclr(i) .LE. 0. ) THEN
1171      drainclr = 0.
1172      dsnowclr = 0.
1173    ELSE
1174      drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i)
1175      dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i)
1176    ENDIF
1177  !--If the cloud narrows
1178  ELSEIF ( dprecipfraccld .LT. 0. ) THEN
1179    !--We reattribute some of the cloudy precip flux
1180    !--to the CS precip flux, proportionnally to the
1181    !--decrease of the cloud precip fraction
1182    draincld = dprecipfraccld / precipfraccld(i) * raincld(i)
1183    dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i)
1184    drainclr = - draincld
1185    dsnowclr = - dsnowcld
1186  !--If the cloud stays the same or if there is no cloud above and
1187  !--in the current layer, nothing happens
1188  ELSE
1189    drainclr = 0.
1190    dsnowclr = 0.
1191  ENDIF
1192
1193  !--We add the tendencies
1194  !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1195  precipfraccld(i) = precipfraccld(i) + dprecipfraccld
1196  precipfracclr(i) = precipfracclr(i) + dprecipfracclr
1197  rainclr(i) = MAX(0., rainclr(i) + drainclr)
1198  snowclr(i) = MAX(0., snowclr(i) + dsnowclr)
1199  raincld(i) = MAX(0., raincld(i) - drainclr)
1200  snowcld(i) = MAX(0., snowcld(i) - dsnowclr)
1201 
1202  !--If vertical heterogeneity is taken into account, we use
1203  !--the "true" volume fraction instead of a modified
1204  !--surface fraction (which is larger and artificially
1205  !--reduces the in-cloud water).
1206  IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
1207    eff_cldfra = ctot_vol(i)
1208  ELSE
1209    eff_cldfra = cldfra(i)
1210  ENDIF
1211
1212
1213  !--Start precipitation growth processes
1214
1215  !--If the cloud is big enough, the precipitation processes activate
1216  ! TODO met on seuil_neb ici ?
1217  IF ( cldfra(i) .GE. seuil_neb ) THEN
1218
1219    !---------------------------------------------------------
1220    !--           COLLECTION AND AGGREGATION
1221    !---------------------------------------------------------
1222    !--Collection: processus through which rain collects small liquid droplets
1223    !--in suspension, and add it to the rain flux
1224    !--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
1225    !--Those processes are treated before autoconversion because we do not
1226    !--want to collect/aggregate the newly formed fluxes, which already
1227    !--"saw" the cloud as they come from it
1228    !--The formulas come from Muench and Lohmann 2020
1229
1230    !--gamma_col: tuning coefficient [-]
1231    !--rho_rain: volumic mass of rain [kg/m3]
1232    !--r_rain: size of the rain droplets [m]
1233    !--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
1234    !--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
1235    !--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
1236    !--get in-cloud mean quantities. The two divisions by eff_cldfra are
1237    !--then simplified.
1238
1239    !--The collection efficiency is perfect.
1240    Eff_rain_liq = 1.
1241    coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
1242    IF ( raincld(i) .GT. 0. ) THEN
1243      !--Exact explicit version, which does not need a barrier because of
1244      !--the exponential decrease
1245      dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
1246
1247      !--Add tendencies
1248      qliq(i) = qliq(i) + dqlcol
1249      raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i)
1250
1251      !--Diagnostic tendencies
1252      dqrcol(i)  = - dqlcol  / dtime
1253    ENDIF
1254
1255    !--Same as for aggregation
1256    !--Eff_snow_liq formula:
1257    !--it s a product of a collection efficiency and a sticking efficiency
1258    ! Milbrandt and Yau formula that gives very low values:
1259    ! Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
1260    ! Lin 1983's formula
1261    Eff_snow_ice = EXP( 0.025 * MIN( ( temp(i) - RTT ), 0.) )
1262    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
1263    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
1264    coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
1265    IF ( snowcld(i) .GT. 0. ) THEN
1266      !--Exact explicit version, which does not need a barrier because of
1267      !--the exponential decrease
1268      dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
1269
1270      !--Add tendencies
1271      qice(i) = qice(i) + dqiagg
1272      snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i)
1273
1274      !--Diagnostic tendencies
1275      dqsagg(i)  = - dqiagg  / dtime
1276    ENDIF
1277
1278
1279    !---------------------------------------------------------
1280    !--                  AUTOCONVERSION
1281    !---------------------------------------------------------
1282    !--Autoconversion converts liquid droplets/ice crystals into
1283    !--rain drops/snowflakes. It relies on the formulations by
1284    !--Sundqvist 1978.
1285
1286    !--If we are in a convective point, we have different parameters
1287    !--for the autoconversion
1288    IF ( ptconv(i) ) THEN
1289        qthresh_auto_rain = cld_lc_con
1290        qthresh_auto_snow = cld_lc_con_snow
1291
1292        tau_auto_rain = cld_tau_con
1293        !--tau for snow depends on the ice fraction in mixed-phase clouds     
1294        tau_auto_snow = tau_auto_snow_max &
1295                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
1296
1297        expo_auto_rain = cld_expo_con
1298        expo_auto_snow = cld_expo_con
1299    ELSE
1300        qthresh_auto_rain = cld_lc_lsc
1301        qthresh_auto_snow = cld_lc_lsc_snow
1302
1303        tau_auto_rain = cld_tau_lsc
1304        !--tau for snow depends on the ice fraction in mixed-phase clouds     
1305        tau_auto_snow = tau_auto_snow_max &
1306                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
1307
1308        expo_auto_rain = cld_expo_lsc
1309        expo_auto_snow = cld_expo_lsc
1310    ENDIF
1311
1312
1313    ! Liquid water quantity to remove according to (Sundqvist, 1978)
1314    ! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) )
1315    !
1316    !--And same formula for ice
1317    !
1318    !--We first treat the second term (with exponential) in an explicit way
1319    !--and then treat the first term (-q/tau) in an exact way
1320
1321    dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
1322               - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
1323
1324    dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
1325               - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
1326
1327
1328    !--Barriers so that we don't create more rain/snow
1329    !--than there is liquid/ice
1330    dqlauto = MAX( - qliq(i), dqlauto )
1331    dqiauto = MAX( - qice(i), dqiauto )
1332
1333    !--Add tendencies
1334    qliq(i) = qliq(i) + dqlauto
1335    qice(i) = qice(i) + dqiauto
1336    raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
1337    snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
1338
1339    !--Diagnostic tendencies
1340    dqsauto(i) = - dqiauto / dtime
1341    dqrauto(i) = - dqlauto / dtime
1342
1343
1344    !---------------------------------------------------------
1345    !--                  RIMING
1346    !---------------------------------------------------------
1347    !--Process which converts liquid droplets in suspension into
1348    !--snow because of the collision between
1349    !--those and falling snowflakes.
1350    !--The formula comes from Muench and Lohmann 2020
1351    !--NB.: this process needs a temperature adjustment
1352
1353    !--Eff_snow_liq formula: following Ferrier 1994,
1354    !--assuming 1
1355    Eff_snow_liq = 1.0
1356    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
1357    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
1358    coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
1359    IF ( snowcld(i) .GT. 0. ) THEN
1360      !--Exact version, which does not need a barrier because of
1361      !--the exponential decrease
1362      dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. )
1363
1364      !--Add tendencies
1365      !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1366      qliq(i) = qliq(i) + dqlrim
1367      snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i)
1368
1369      !--Temperature adjustment with the release of latent
1370      !--heat because of solid condensation
1371      temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
1372                        / ( 1. + RVTMP2 * qtot(i) )
1373
1374      !--Diagnostic tendencies
1375      dqsrim(i)  = - dqlrim  / dtime
1376    ENDIF
1377
1378  ENDIF ! cldfra .GE. seuil_neb
1379
1380ENDDO ! loop on klon
1381
1382
1383!--Re-calculation of saturation specific humidity
1384!--because riming changed temperature
1385CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
1386
1387DO i = 1, klon
1388
1389  !---------------------------------------------------------
1390  !--                  MELTING
1391  !---------------------------------------------------------
1392  !--Process through which snow melts into rain.
1393  !--The formula is homemade.
1394  !--NB.: this process needs a temperature adjustment
1395
1396  !--dqsmelt_max         : maximum snow melting so that temperature
1397  !--                      stays higher than 273 K [kg/kg]
1398  !--capa_snowflake      : capacitance of a snowflake, equal to
1399  !--                      the radius if the snowflake is a sphere [m]
1400  !--temp_wetbulb        : wet-bulb temperature [K]
1401  !--snow_fallspeed      : snow fall velocity (in clear/cloudy sky) [m/s]
1402  !--air_thermal_conduct : thermal conductivity of the air [J/m/K/s]
1403  !--gamma_melt          : tuning parameter for melting [-]
1404  !--nb_snowflake        : number of snowflakes (in clear/cloudy air) [-]
1405
1406  IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
1407    !--Computed according to
1408    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1409    dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
1410                        * ( 1. + RVTMP2 * qtot(i) ))
1411   
1412    !--Initialisation
1413    dqsclrmelt = 0.
1414    dqscldmelt = 0.
1415
1416    !--We assume that the snowflakes are spherical
1417    capa_snowflake = r_snow
1418    !--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
1419    air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184   
1420    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
1421    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
1422
1423    !--In clear air
1424    IF ( ( snowclr(i) .GT. 0. ) .AND.  ( precipfracclr(i) .GT. 0. ) ) THEN
1425      !--Formula for the wet-bulb temperature from ECMWF (IFS)
1426      !--The vapor used is the vapor in the clear sky
1427      temp_wetbulb = temp(i) &
1428                   - ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) &
1429                   * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
1430                   - 40.637 * ( temp(i) - 275. ) )
1431      !--Calculated according to
1432      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
1433      nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr &
1434                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
1435      dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
1436                   * capa_snowflake / RLMLT * gamma_melt &
1437                   * MAX(0., temp_wetbulb - RTT) * dtime
1438     
1439      !--Barrier to limit the melting flux to the clr snow flux in the mesh
1440      dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i))
1441    ENDIF
1442
1443
1444    !--In cloudy air
1445    IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
1446      !--As the air is saturated, the wet-bulb temperature is equal to the
1447      !--temperature
1448      temp_wetbulb = temp(i)
1449      !--Calculated according to
1450      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
1451      nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld &
1452                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
1453      dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
1454                   * capa_snowflake / RLMLT * gamma_melt &
1455                   * MAX(0., temp_wetbulb - RTT) * dtime
1456
1457      !--Barrier to limit the melting flux to the cld snow flux in the mesh
1458      dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i))
1459   ENDIF
1460   
1461
1462    !--Barrier on temperature. If the total melting flux leads to a
1463    !--positive temperature, it is limited to keep temperature above 0 degC.
1464    !--It is activated if the total is LOWER than the max
1465    !--because everything is negative
1466    dqstotmelt = dqsclrmelt + dqscldmelt
1467    IF ( dqstotmelt .LT. dqsmelt_max ) THEN
1468      !--We redistribute the max melted snow keeping
1469      !--the clear/cloud partition of the melted snow
1470      dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
1471      dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
1472      dqstotmelt = dqsmelt_max
1473
1474    ENDIF
1475
1476    !--Add tendencies
1477    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1478    rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i))
1479    raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i))
1480    snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i))
1481    snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i))
1482
1483    !--Temperature adjustment with the release of latent
1484    !--heat because of melting
1485    temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
1486                      / ( 1. + RVTMP2 * qtot(i) )
1487
1488    !--Diagnostic tendencies
1489    dqrmelt(i) = - dqstotmelt / dtime
1490    dqsmelt(i) = dqstotmelt / dtime
1491
1492  ENDIF
1493
1494
1495  !---------------------------------------------------------
1496  !--                  FREEZING
1497  !---------------------------------------------------------
1498  !--Process through which rain freezes into snow.
1499  !-- We parameterize it as a 2 step process:
1500  !--first: freezing following collision with ice crystals
1501  !--second: immersion freezing following (inspired by Bigg 1953)
1502  !--the latter is parameterized as an exponential decrease of the rain
1503  !--water content with a homemade formulya
1504  !--This is based on a caracteritic time of freezing, which
1505  !--exponentially depends on temperature so that it is
1506  !--equal to 1 for temp_nowater (see below) and is close to
1507  !--0 for RTT (=273.15 K).
1508  !--NB.: this process needs a temperature adjustment
1509  !--dqrfreez_max : maximum rain freezing so that temperature
1510  !--              stays lower than 273 K [kg/kg]
1511  !--tau_freez    : caracteristic time of freezing [s]
1512  !--gamma_freez  : tuning parameter [s-1]
1513  !--alpha_freez  : tuning parameter for the shape of the exponential curve [-]
1514  !--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC)
1515
1516  IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
1517
1518 
1519    !--1st step: freezing following collision with ice crystals
1520    !--Sub-process of freezing which quantifies the collision between
1521    !--ice crystals in suspension and falling rain droplets.
1522    !--The rain droplets freeze, becoming graupel, and carrying
1523    !--the ice crystal (which acted as an ice nucleating particle).
1524    !--The formula is adapted from the riming formula.
1525    !--it works only in the cloudy part
1526
1527    dqifreez = 0.
1528    dqrtotfreez_step1 = 0.
1529
1530    IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. &
1531         ( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
1532      dqrclrfreez = 0.
1533      dqrcldfreez = 0.
1534
1535      !--Computed according to
1536      !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1537      dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
1538                         * ( 1. + RVTMP2 * qtot(i) ))
1539 
1540
1541      !--The collision efficiency is assumed unity
1542      Eff_rain_ice = 1.
1543      coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
1544      !--Exact version, which does not need a barrier because of
1545      !--the exponential decrease.
1546      dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
1547
1548      !--We add the part of rain water that freezes, limited by a temperature barrier
1549      !--This quantity is calculated assuming that the number of drop that freeze correspond to the number
1550      !--of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
1551      !--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3
1552      !--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3
1553      !--The assumption above corresponds to dNi = dNr, i.e.,
1554      !-- dqr = dqi * (4/3 PI rho_rain * r_rain**3) / (4/3 PI rho_ice * r_ice**3)
1555      !--Dry density [kg/m3]
1556      rho = pplay(i) / temp(i) / RD
1557      !--r_ice formula from Sun and Rikus (1999)
1558      r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 &
1559            + 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
1560      dqrcldfreez = dqifreez * rho_rain * r_rain**3. / ( rho_ice * r_ice**3. )
1561      dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i))
1562      dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max)
1563      dqrtotfreez_step1 = dqrcldfreez
1564
1565      !--Add tendencies
1566      !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1567      qice(i) = qice(i) + dqifreez
1568      raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
1569      snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i))
1570      temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD &
1571                        / ( 1. + RVTMP2 * qtot(i) )
1572
1573    ENDIF
1574   
1575    !-- Second step immersion freezing of rain drops
1576    !-- with a homemade timeconstant depending on temperature
1577
1578    dqrclrfreez = 0.
1579    dqrcldfreez = 0.
1580    dqrtotfreez_step2 = 0.
1581    !--Computed according to
1582    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1583
1584    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
1585                         * ( 1. + RVTMP2 * qtot(i) ))
1586
1587   
1588    tau_freez = 1. / ( beta_freez &
1589              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
1590
1591
1592    !--In clear air
1593    IF ( rainclr(i) .GT. 0. ) THEN
1594      !--Exact solution of dqrain/dt = -qrain/tau_freez
1595      dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
1596    ENDIF
1597
1598    !--In cloudy air
1599    IF ( raincld(i) .GT. 0. ) THEN
1600      !--Exact solution of dqrain/dt = -qrain/tau_freez
1601      dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
1602    ENDIF
1603
1604    !--temperature barrier step 2
1605    !--It is activated if the total is LOWER than the max
1606    !--because everything is negative
1607    dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
1608 
1609    IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
1610      !--We redistribute the max freezed rain keeping
1611      !--the clear/cloud partition of the freezing rain
1612      dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
1613      dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
1614      dqrtotfreez_step2 = dqrfreez_max
1615    ENDIF
1616
1617
1618    !--Add tendencies
1619    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1620    rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i))
1621    raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
1622    snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i))
1623    snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i))
1624
1625
1626    !--Temperature adjustment with the uptake of latent
1627    !--heat because of freezing
1628    temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
1629                      / ( 1. + RVTMP2 * qtot(i) )
1630
1631    !--Diagnostic tendencies
1632    dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2         
1633    dqrfreez(i) = dqrtotfreez / dtime
1634    dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
1635
1636  ENDIF
1637
1638
1639
1640  !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
1641  !--we reduce the precipiration fraction in the clear/cloudy air so that the new
1642  !--local flux of rain+snow is equal to rain_int_min.
1643  !--Here, rain+snow is the gridbox-mean flux of precip.
1644  !--Therefore, (rain+snow)/precipfrac is the local flux of precip.
1645  !--If the local flux of precip is lower than rain_int_min, i.e.,
1646  !-- (rain+snow)/precipfrac < rain_int_min , i.e.,
1647  !-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
1648  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
1649  !--Note that this is physically different than what is proposed in LTP thesis.
1650  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
1651  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
1652
1653  !--Calculate outputs
1654  rain(i) = rainclr(i) + raincld(i)
1655  snow(i) = snowclr(i) + snowcld(i)
1656
1657  !--Diagnostics
1658  !--BEWARE this is indeed a diagnostic: this is an estimation from
1659  !--the value of the flux at the bottom interface of the mesh and
1660  !--and assuming an upstream numerical calculation
1661  !--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is
1662  !--used for computing the total ice water content in the mesh
1663  !--for radiation only
1664  qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr &
1665                 + raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld )
1666  qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr &
1667                 + snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld )
1668
1669
1670ENDDO ! loop on klon
1671
1672END SUBROUTINE poprecip_postcld
1673
1674END MODULE lmdz_lscp_precip
Note: See TracBrowser for help on using the repository browser.