Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

lmdz_lscp_precip.f90

Last change on this file was 5431, checked in by aborella, 4 hours ago
Merge with trunk
File size: 79.1 KB

Line
1	MODULE 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
10	IMPLICIT NONE
11
12	CONTAINS
13
14	!----------------------------------------------------------------
15	! historical (till CMIP6) version of the pre-cloud formation
16	! precipitation scheme containing precip evaporation and melting
17
18	SUBROUTINE 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
24	USE lmdz_lscp_ini, ONLY : iflag_evap_prec
25	USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, ztfondue
26	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
27	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
28
29	IMPLICIT NONE
30
31
32	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
33	REAL, INTENT(IN) :: dtime !--time step [s]
34	LOGICAL, INTENT(IN) :: iftop !--if top of the column
35
36
37	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
38	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
39	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
40
41	REAL, INTENT(INOUT), DIMENSION(klon) :: zt !--current temperature [K]
42	REAL, INTENT(IN), DIMENSION(klon) :: ztupnew !--updated temperature of the overlying layer [K]
43
44	REAL, INTENT(INOUT), DIMENSION(klon) :: zq !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
45	REAL, INTENT(INOUT), DIMENSION(klon) :: zmqc !--specific humidity in the precipitation falling from the upper layer [kg/kg]
46
47	REAL, INTENT(IN), DIMENSION(klon) :: zneb !--cloud fraction IN THE LAYER ABOVE [-]
48	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecip !--fraction of precipitation IN THE LAYER ABOVE [-]
49	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
50
51	REAL, INTENT(INOUT), DIMENSION(klon) :: zrfl !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
52	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
53	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflcld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
54	REAL, INTENT(INOUT), DIMENSION(klon) :: zifl !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
55	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
56	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
57
58	REAL, INTENT(OUT), DIMENSION(klon) :: dqreva !--rain tendency due to evaporation [kg/kg/s]
59	REAL, INTENT(OUT), DIMENSION(klon) :: dqssub !--snow tendency due to sublimation [kg/kg/s]
60
61
62	REAL :: zmair
63	REAL :: zcpair, zcpeau
64	REAL, DIMENSION(klon) :: qzero
65	REAL, DIMENSION(klon) :: zqs, zdqs
66	REAL, DIMENSION(klon) :: qsl, qsi
67	REAL, DIMENSION(klon) :: dqsl, dqsi
68
69	REAL :: zqev0
70	REAL :: zqev, zqevt
71	REAL :: zqevi, zqevti
72
73	REAL, DIMENSION(klon) :: zrfln, zifln
74
75	REAL :: zmelt
76	INTEGER :: i
77
78	qzero(:) = 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
94	IF (iftop) THEN
95
96	DO i = 1, klon
97	zmqc(i) = 0.
98	ENDDO
99
100	ELSE
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+RVTMP2zq(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
120	ENDIF
121
122	! --------------------------------------------------------------------
123	! P2> Precipitation evaporation/sublimation/melting
124	! --------------------------------------------------------------------
125	! A part of the precipitation coming from above is evaporated/sublimated/melted.
126	! --------------------------------------------------------------------
127
128	IF (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
303	ENDIF ! (iflag_evap_prec>=1)
304
305	END SUBROUTINE histprecip_precld
306
307
308	SUBROUTINE 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
316	USE lmdz_lscp_ini, ONLY : RD, RG, RCPD, RVTMP2, RLSTT, RLMLT, RTT
317	USE 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
320	USE 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
323	USE lmdz_lscp_tools, ONLY : fallice_velocity
324
325	IMPLICIT NONE
326
327
328	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
329	REAL, INTENT(IN) :: dtime !--time step [s]
330	LOGICAL, INTENT(IN) :: iftop !--if top of the column
331
332	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
333	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
334	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
335	REAL, INTENT(IN), DIMENSION(klon) :: ctot_vol !--volumic cloud fraction [-]
336	LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv !--true if we are in a convective point
337	REAL, INTENT(IN), DIMENSION(klon) :: zdqsdT_raw !--derivative of qsat w.r.t. temperature times L/Cp [SI]
338
339	REAL, INTENT(INOUT), DIMENSION(klon) :: zt !--current temperature [K]
340	REAL, INTENT(INOUT), DIMENSION(klon) :: zq !--current water vapor specific humidity [kg/kg]
341	REAL, INTENT(INOUT), DIMENSION(klon) :: zoliq !--current liquid+ice water specific humidity [kg/kg]
342	REAL, INTENT(INOUT), DIMENSION(klon) :: zoliql !--current liquid water specific humidity [kg/kg]
343	REAL, INTENT(INOUT), DIMENSION(klon) :: zoliqi !--current ice water specific humidity [kg/kg]
344	REAL, INTENT(INOUT), DIMENSION(klon) :: zcond !--liquid+ice water specific humidity AFTER CONDENSATION (no precip process) [kg/kg]
345	REAL, INTENT(IN), DIMENSION(klon) :: zfice !--ice fraction AFTER CONDENSATION [-]
346	REAL, INTENT(IN), DIMENSION(klon) :: zmqc !--specific humidity in the precipitation falling from the upper layer [kg/kg]
347
348	REAL, INTENT(IN), DIMENSION(klon) :: rneb !--cloud fraction [-]
349	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
350	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipcld !--fraction of precipitation in the cloud IN THE LAYER ABOVE [-]
351	REAL, INTENT(INOUT), DIMENSION(klon) :: zneb !--cloud fraction IN THE LAYER ABOVE [-]
352	REAL, INTENT(INOUT), DIMENSION(klon) :: tot_zneb !--total cloud cover above the current mesh [-]
353	REAL, INTENT(INOUT), DIMENSION(klon) :: zrho_up !--air density IN THE LAYER ABOVE [kg/m3]
354	REAL, INTENT(INOUT), DIMENSION(klon) :: zvelo_up !--ice fallspeed velocity IN THE LAYER ABOVE [m/s]
355
356	REAL, INTENT(INOUT), DIMENSION(klon) :: zrfl !--flux of rain gridbox-mean [kg/s/m2]
357	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflclr !--flux of rain gridbox-mean in clear sky [kg/s/m2]
358	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflcld !--flux of rain gridbox-mean in cloudy air [kg/s/m2]
359	REAL, INTENT(INOUT), DIMENSION(klon) :: zifl !--flux of snow gridbox-mean [kg/s/m2]
360	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflclr !--flux of snow gridbox-mean in clear sky [kg/s/m2]
361	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflcld !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
362
363	REAL, INTENT(OUT), DIMENSION(klon) :: zradocond !--condensed water used in the radiation scheme [kg/kg]
364	REAL, INTENT(OUT), DIMENSION(klon) :: zradoice !--condensed ice water used in the radiation scheme [kg/kg]
365	REAL, INTENT(OUT), DIMENSION(klon) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
366	REAL, 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
370	REAL :: smallestreal
371	REAL, DIMENSION(klon) :: tot_znebn, d_tot_zneb
372	REAL, DIMENSION(klon) :: d_znebprecip_cld_clr, d_znebprecip_clr_cld
373	REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
374	REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
375
376	! Local variables for autoconversion
377	REAL :: zct, zcl, zexpo, ffallv
378	REAL :: zchau, zfroi
379	REAL :: zrain, zsnow, zprecip
380	REAL :: effective_zneb
381	REAL, DIMENSION(klon) :: zrho, zvelo
382	REAL, DIMENSION(klon) :: zdz, iwc
383
384	! Local variables for Bergeron process
385	REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
386	REAL, DIMENSION(klon) :: zqpreci, zqprecl
387
388	! Local variables for calculation of quantity of condensates seen by the radiative scheme
389	REAL, DIMENSION(klon) :: zradoliq
390	REAL, DIMENSION(klon) :: ziflprev
391
392	! Misc
393	INTEGER :: i, n
394
395	! Initialisation
396	smallestreal=1.e-9
397	zqprecl(:) = 0.
398	zqpreci(:) = 0.
399	ziflprev(:) = 0.
400
401
402	IF (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
446	ENDIF
447
448
449	! Autoconversion
450	! -------------------------------------------------------------------------------
451
452	! Initialisation
453	DO 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)
465	ENDDO
466
467
468	DO 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 = zctdtime/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/rhod(rhowice*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_velodtime/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
585	ENDDO ! n = 1,niter
586
587	! Precipitation flux calculation
588
589	DO 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/zcpzdqsdT_raw(i)
618	! Computation of DT if all the liquid precip freezes
619	DeltaT = RLMLTzqprecl(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))/(RGdtime)
641	ziflcld(i) = ziflcld(i)+ zqpreci(i) &
642	(paprsdn(i)-paprsup(i))/(RGdtime)
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))/(RGdtime)
649	zifl(i) = zifl(i)+ zqpreci(i) &
650	(paprsdn(i)-paprsup(i))/(RGdtime)
651	ENDIF
652	! c_iso : same for isotopes
653
654	ENDIF ! rneb>0
655
656	ENDDO
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
660	IF (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
681	ENDIF
682
683	! we recalculate zradoice to account for contributions from the precipitation flux
684	! if ok_radocond_snow is true
685	IF ( 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(:)
704	ENDIF
705
706	END SUBROUTINE histprecip_postcld
707
708
709	!----------------------------------------------------------------
710	! Computes the processes-oriented precipitation formulations for
711	! evaporation and sublimation
712	!
713	SUBROUTINE poprecip_precld( &
714	klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
715	qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
716	cldfra, rvc_seri, qliq, qice, &
717	rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub &
718	)
719
720	USE lmdz_lscp_ini, ONLY : prt_level, lunout
721	USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
722	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
723	USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds
724	USE lmdz_lscp_ini, ONLY : eps, temp_nowater
725	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
726
727	IMPLICIT NONE
728
729
730	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
731	REAL, INTENT(IN) :: dtime !--time step [s]
732	LOGICAL, INTENT(IN) :: iftop !--if top of the column
733
734
735	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
736	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
737	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
738
739	REAL, INTENT(INOUT), DIMENSION(klon) :: temp !--current temperature [K]
740	REAL, INTENT(IN), DIMENSION(klon) :: tempupnew !--updated temperature of the overlying layer [K]
741
742	REAL, INTENT(INOUT), DIMENSION(klon) :: qvap !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
743	REAL, INTENT(INOUT), DIMENSION(klon) :: qprecip !--specific humidity in the precipitation falling from the upper layer [kg/kg]
744
745	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
746	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
747
748	REAL, INTENT(IN), DIMENSION(klon) :: qvapclrup !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
749	REAL, INTENT(IN), DIMENSION(klon) :: qtotupnew !--total specific humidity IN THE LAYER ABOVE [kg/kg]
750
751	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
752	REAL, INTENT(IN), DIMENSION(klon) :: rvc_seri !--cloud water vapor at the beginning of lscp (ratio wrt total water) - used only if the cloud properties are advected [kg/kg]
753	REAL, INTENT(IN), DIMENSION(klon) :: qliq !--liquid water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
754	REAL, INTENT(IN), DIMENSION(klon) :: qice !--ice water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
755
756
757	REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
758	REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
759	REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
760	REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
761	REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
762	REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
763
764	REAL, INTENT(OUT), DIMENSION(klon) :: dqreva !--rain tendency due to evaporation [kg/kg/s]
765	REAL, INTENT(OUT), DIMENSION(klon) :: dqssub !--snow tendency due to sublimation [kg/kg/s]
766
767
768	!--Integer for interating over klon
769	INTEGER :: i
770	!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
771	REAL, DIMENSION(klon) :: dhum_to_dflux
772	!--Air density [kg/m3] and layer thickness [m]
773	REAL, DIMENSION(klon) :: rho, dz
774	!--Temporary precip fractions and fluxes for the evaporation
775	REAL, DIMENSION(klon) :: precipfracclr_tmp, precipfraccld_tmp
776	REAL, DIMENSION(klon) :: rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp
777
778	!--Saturation values
779	REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
780	!--Vapor in the clear sky and cloud
781	REAL :: qvapclr, qvapcld
782	!--Fluxes tendencies because of evaporation and sublimation
783	REAL :: dprecip_evasub_max, dprecip_evasub_tot
784	REAL :: drainclreva, draincldeva, dsnowclrsub, dsnowcldsub
785	!--Specific humidity tendencies because of evaporation and sublimation
786	REAL :: dqrevap, dqssubl
787	!--Specific heat constant
788	REAL :: cpair, cpw
789
790	!--Initialisation
791	qzero(:) = 0.
792	dqreva(:) = 0.
793	dqssub(:) = 0.
794	dqrevap = 0.
795	dqssubl = 0.
796
797	!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
798	dhum_to_dflux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime
799	rho(:) = pplay(:) / temp(:) / RD
800	dz(:) = ( paprsdn(:) - paprsup(:) ) / RG / rho(:)
801
802	!--Calculation of saturation specific humidity
803	!--depending on temperature:
804	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:))
805	!--wrt liquid water
806	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
807	!--wrt ice
808	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
809
810
811
812	!--First step consists in "thermalizing" the layer:
813	!--as the flux of precip from layer above "advects" some heat (as the precip is at the temperature
814	!--of the overlying layer) we recalculate a mean temperature that both the air and the precip in the
815	!--layer have.
816
817	IF (iftop) THEN
818
819	DO i = 1, klon
820	qprecip(i) = 0.
821	ENDDO
822
823	ELSE
824
825	DO i = 1, klon
826	!--No condensed water so cp=cp(vapor+dry air)
827	!-- RVTMP2=rcpv/rcpd-1
828	cpair = RCPD * ( 1. + RVTMP2 * qvap(i) )
829	cpw = RCPD * RVTMP2
830	!--qprecip has to be thermalized with
831	!--layer's air so that precipitation at the ground has the
832	!--same temperature as the lowermost layer
833	!--we convert the flux into a specific quantity qprecip
834	qprecip(i) = ( rain(i) + snow(i) ) / dhum_to_dflux(i)
835	!-- t(i,k+1) + d_t(i,k+1): new temperature of the overlying layer
836	temp(i) = ( tempupnew(i) * qprecip(i) * cpw + cpair * temp(i) ) &
837	/ ( cpair + qprecip(i) * cpw )
838	ENDDO
839
840	ENDIF
841
842
843	!--Initialise the precipitation fractions and fluxes
844	DO i = 1, klon
845	precipfracclr_tmp(i) = precipfracclr(i)
846	precipfraccld_tmp(i) = precipfraccld(i)
847	rainclr_tmp(i) = rainclr(i)
848	snowclr_tmp(i) = snowclr(i)
849	raincld_tmp(i) = raincld(i)
850	snowcld_tmp(i) = snowcld(i)
851	ENDDO
852
853	!--If we relax the LTP assumption that the cloud is the same than in the
854	!--gridbox above, we can change the tmp quantities
855	IF ( ok_ice_supersat ) THEN
856	!--Update the precipitation fraction using advected cloud fraction
857	CALL poprecip_fracupdate( &
858	klon, cldfra, precipfracclr_tmp, precipfraccld_tmp, &
859	rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp)
860
861	! here we could put an ELSE statement and do one iteration of the condensation scheme
862	! (but it would only diagnose cloud from lognormal - link with thermals?)
863
864	DO i = 1, klon
865	!--We cannot have a total fraction greater than the one in the gridbox above
866	!--because new cloud formation has not yet occured.
867	!--If this happens, we reduce the precip frac in the cloud, because it can
868	!--only mean that the cloud fraction at this level is greater than the total
869	!--precipitation fraction of the level above, and the precip frac in clear sky
870	!--is zero.
871	!--NB. this only works because we assume a max-random cloud and precip overlap
872	precipfraccld_tmp(i) = MIN( precipfraccld_tmp(i), precipfracclr(i) + precipfraccld(i) )
873	ENDDO
874
875	ENDIF
876	!--At this stage, we guarantee that
877	!-- precipfracclr + precipfraccld = precipfracclr_tmp + precipfraccld_tmp
878
879
880	DO i = 1, klon
881
882	!--If there is precipitation from the layer above
883	IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
884
885	!--Init
886	drainclreva = 0.
887	draincldeva = 0.
888	dsnowclrsub = 0.
889	dsnowcldsub = 0.
890
891	!--Init the properties of the air where reevaporation occurs
892	IF ( ok_ice_supersat ) THEN
893	!--We can diagnose the water vapor in the cloud using the advected
894	!--water vapor in the cloud and cloud fraction
895	IF ( ok_unadjusted_clouds .AND. ( cldfra(i) .GT. eps ) ) THEN
896	qvapcld = rvc_seri(i) * qvap(i) / cldfra(i)
897	ENDIF
898	!--We can diagnose completely the water vapor in clear sky, because all
899	!--the needed variables (ice, liq, vapor in cld, cloud fraction) are
900	!--advected
901	!--Note that qvap(i) is the total water in the gridbox
902	IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
903	qvapclr = ( qvap(i) - qice(i) - qliq(i) - rvc_seri(i) * qvap(i) ) / ( 1. - cldfra(i) )
904	ENDIF
905	ELSEIF ( ok_corr_vap_evasub ) THEN
906	!--Corrected version from Ludo - we use the same water ratio between
907	!--the clear and the cloudy sky as in the layer above. This
908	!--extends the assumption that the cloud fraction is the same
909	!--as the layer above. This is assumed only for the evap / subl
910	!--process
911	!--Note that qvap(i) is the total water in the gridbox, and
912	!--precipfraccld(i) is the cloud fraction in the layer above
913	IF ( ( 1. - precipfraccld(i) ) .GT. eps ) THEN
914	qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
915	ENDIF
916	ELSE
917	!--Legacy version from Ludo - we use the total specific humidity
918	!--for the evap / subl process
919	qvapclr = qvap(i)
920	ENDIF
921
922	!---------------------------
923	!-- EVAP / SUBL IN CLEAR SKY
924	!---------------------------
925
926	IF ( precipfracclr_tmp(i) .GT. eps ) THEN
927
928	!--Evaporation of liquid precipitation coming from above
929	!--in the clear sky only
930	!--dprecip/dz = -beta(1-qvap/qsat)(precip**expo_eva)
931	!--formula from Sundqvist 1988, Klemp & Wilhemson 1978
932	!--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly)
933	!--which does not need a barrier on rainclr, because included in the formula
934	drainclreva = precipfracclr_tmp(i) * MAX(0., &
935	- coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) &
936	+ ( rainclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_eva ) &
937	)**( 1. / ( 1. - expo_eva ) ) - rainclr_tmp(i)
938
939	!--Evaporation is limited by 0
940	!--NB. with ok_ice_supersat activated, this barrier should be useless
941	drainclreva = MIN(0., drainclreva)
942
943
944	!--Sublimation of the solid precipitation coming from above
945	!--(same formula as for liquid precip)
946	!--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly)
947	!--which does not need a barrier on snowclr, because included in the formula
948	dsnowclrsub = precipfracclr_tmp(i) * MAX(0., &
949	- coef_sub * ( 1. - expo_sub ) * (1. - qvapclr / qsati(i)) * dz(i) &
950	+ ( snowclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_sub ) &
951	)**( 1. / ( 1. - expo_sub ) ) - snowclr_tmp(i)
952
953	!--If ice supersaturation is activated, we allow vapor to condense on falling snow
954	!--i.e., having a positive dsnowclrsub
955	IF ( .NOT. ok_ice_supersat ) THEN
956	!--Sublimation is limited by 0
957	dsnowclrsub = MIN(0., dsnowclrsub)
958	ENDIF
959
960
961	!--Evaporation limit: we ensure that the layer's fraction below
962	!--the clear sky does not reach saturation. In this case, we
963	!--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
964	!--Max evaporation is computed not to saturate the clear sky precip fraction
965	!--(i.e., the fraction where evaporation occurs)
966	!--It is expressed as a max flux dprecip_evasub_max
967
968	IF ( .NOT. ok_ice_supersat ) THEN
969	dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr_tmp(i)) &
970	* dhum_to_dflux(i)
971	ELSE
972	dprecip_evasub_max = ( qvapclr - qsat(i) ) * precipfracclr_tmp(i) &
973	* dhum_to_dflux(i)
974	ENDIF
975	dprecip_evasub_tot = drainclreva + dsnowclrsub
976
977	!--Barriers
978	!--If activates if the total is LOWER than the max because
979	!--everything is negative
980	IF ( (( dprecip_evasub_max .LT. 0. ) .AND. &
981	( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. &
982	(( dprecip_evasub_max .GT. 0. ) .AND. &
983	( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN
984	drainclreva = dprecip_evasub_max * drainclreva / dprecip_evasub_tot
985	dsnowclrsub = dprecip_evasub_max * dsnowclrsub / dprecip_evasub_tot
986	ENDIF
987
988	ENDIF ! precipfracclr_tmp .GT. eps
989
990
991	!---------------------------
992	!-- EVAP / SUBL IN THE CLOUD
993	!---------------------------
994
995	IF ( ok_unadjusted_clouds .AND. ( temp(i) .LE. temp_nowater ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
996	!--Evaporation of liquid precipitation coming from above
997	!--in the cloud only
998	!--dprecip/dz = -beta(1-qvap/qsat)(precip**expo_eva)
999	!--formula from Sundqvist 1988, Klemp & Wilhemson 1978
1000	!--Exact explicit formulation (raincld is resolved exactly, qvap explicitly)
1001	!--which does not need a barrier on raincld, because included in the formula
1002	!draincldeva = precipfraccld_tmp(i) * MAX(0., &
1003	! - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
1004	! + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
1005	! )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
1006
1007	!--Evaporation is limited by 0
1008	!--NB. with ok_ice_supersat activated, this barrier should be useless
1009	!draincldeva = MIN(0., draincldeva)
1010	draincldeva = 0.
1011
1012
1013	!--Sublimation of the solid precipitation coming from above
1014	!--(same formula as for liquid precip)
1015	!--Exact explicit formulation (snowcld is resolved exactly, qvap explicitly)
1016	!--which does not need a barrier on snowcld, because included in the formula
1017	dsnowcldsub = precipfraccld_tmp(i) * MAX(0., &
1018	- coef_sub * ( 1. - expo_sub ) * (1. - qvapcld / qsati(i)) * dz(i) &
1019	+ ( snowcld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_sub ) &
1020	)**( 1. / ( 1. - expo_sub ) ) - snowcld_tmp(i)
1021
1022	!--There is no barrier because we want deposition to occur if it is possible
1023
1024
1025	!--Evaporation limit: we ensure that the layer's fraction below
1026	!--the clear sky does not reach saturation. In this case, we
1027	!--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
1028	!--Max evaporation is computed not to saturate the clear sky precip fraction
1029	!--(i.e., the fraction where evaporation occurs)
1030	!--It is expressed as a max flux dprecip_evasub_max
1031
1032	dprecip_evasub_max = ( qvapcld - qsat(i) ) * precipfraccld_tmp(i) * dhum_to_dflux(i)
1033	dprecip_evasub_tot = draincldeva + dsnowcldsub
1034
1035	!--Barriers
1036	!--If activates if the total is LOWER than the max because
1037	!--everything is negative
1038	IF ( (( dprecip_evasub_max .LT. 0. ) .AND. &
1039	( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. &
1040	(( dprecip_evasub_max .GT. 0. ) .AND. &
1041	( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN
1042	draincldeva = dprecip_evasub_max * draincldeva / dprecip_evasub_tot
1043	dsnowcldsub = dprecip_evasub_max * dsnowcldsub / dprecip_evasub_tot
1044	ENDIF
1045
1046	ENDIF ! ok_unadjusted_clouds
1047
1048
1049	!--New solid and liquid precipitation fluxes after evap and sublimation
1050	dqrevap = ( drainclreva + draincldeva ) / dhum_to_dflux(i)
1051	dqssubl = ( dsnowclrsub + dsnowcldsub ) / dhum_to_dflux(i)
1052
1053	!--Vapor is updated after evaporation/sublimation (it is increased)
1054	qvap(i) = qvap(i) - dqrevap - dqssubl
1055	!--qprecip is the total condensed water in the precip flux (it is decreased)
1056	qprecip(i) = qprecip(i) + dqrevap + dqssubl
1057	!--Air and precip temperature (i.e., gridbox temperature)
1058	!--is updated due to latent heat cooling
1059	temp(i) = temp(i) &
1060	+ dqrevap * RLVTT / RCPD &
1061	/ ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) &
1062	+ dqssubl * RLSTT / RCPD &
1063	/ ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
1064
1065	!--Add tendencies
1066	!--The MAX is needed because in some cases, the flux can be slightly
1067	!--negative (numerical precision)
1068	rainclr_tmp(i) = MAX(0., rainclr_tmp(i) + drainclreva)
1069	snowclr_tmp(i) = MAX(0., snowclr_tmp(i) + dsnowclrsub)
1070	raincld_tmp(i) = MAX(0., raincld_tmp(i) + draincldeva)
1071	snowcld_tmp(i) = MAX(0., snowcld_tmp(i) + dsnowcldsub)
1072
1073
1074	!--We reattribute fluxes according to initial fractions, using proportionnality
1075	!--Note that trough this process, we conserve total precip fluxes
1076	!-- rainclr_tmp + raincld_tmp = rainclr + raincld (same for snow)
1077	!--If the cloud has increased, a part of the precip in clear sky falls in clear sky,
1078	!--the rest falls in the cloud
1079	IF ( ( precipfraccld(i) .GT. eps ) .AND. ( precipfraccld_tmp(i) .GT. precipfraccld(i) ) ) THEN
1080	!--All the precip from cloud falls in the cloud
1081	raincld(i) = MAX(0., raincld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i))
1082	rainclr(i) = MAX(0., rainclr_tmp(i) + raincld_tmp(i) - raincld(i))
1083	snowcld(i) = MAX(0., snowcld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i))
1084	snowclr(i) = MAX(0., snowclr_tmp(i) + snowcld_tmp(i) - snowcld(i))
1085
1086	!--If the cloud has narrowed, a part of the precip in cloud falls in clear sky,
1087	!--the rest falls in the cloud
1088	ELSEIF ( ( precipfracclr(i) .GT. eps ) .AND. ( precipfracclr_tmp(i) .GT. precipfracclr(i) ) ) THEN
1089	!--All the precip from clear sky falls in clear sky
1090	rainclr(i) = MAX(0., rainclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i))
1091	raincld(i) = MAX(0., raincld_tmp(i) + rainclr_tmp(i) - rainclr(i))
1092	snowclr(i) = MAX(0., snowclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i))
1093	snowcld(i) = MAX(0., snowcld_tmp(i) + snowclr_tmp(i) - snowclr(i))
1094
1095	ELSE
1096	!--If the cloud stays the same or if there is no cloud above and
1097	!--in the current layer, there is no reattribution
1098	rainclr(i) = rainclr_tmp(i)
1099	raincld(i) = raincld_tmp(i)
1100	snowclr(i) = snowclr_tmp(i)
1101	snowcld(i) = snowcld_tmp(i)
1102	ENDIF
1103
1104	!--If there is no more precip fluxes, the precipitation fraction is set to 0
1105	IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0.
1106	IF ( ( raincld(i) + snowcld(i) ) .LE. 0. ) precipfraccld(i) = 0.
1107
1108	!--Calculation of the total fluxes
1109	rain(i) = rainclr(i) + raincld(i)
1110	snow(i) = snowclr(i) + snowcld(i)
1111
1112	ELSEIF ( ( rain(i) + snow(i) ) .LE. 0. ) THEN
1113	!--If no precip, we reinitialize the cloud fraction used for the precip to 0
1114	precipfraccld(i) = 0.
1115	precipfracclr(i) = 0.
1116
1117	ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. )
1118
1119	!--Diagnostic tendencies
1120	dqssub(i) = dqssubl / dtime
1121	dqreva(i) = dqrevap / dtime
1122
1123	ENDDO ! loop on klon
1124
1125	END SUBROUTINE poprecip_precld
1126
1127
1128	!----------------------------------------------------------------
1129	! Computes the precipitation fraction update
1130	! The goal of this routine is to reattribute precipitation fractions
1131	! and fluxes to clear or cloudy air, depending on the variation of
1132	! the cloud fraction on the vertical dimension. We assume a
1133	! maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
1134	! and LTP thesis, 2021)
1135	! NB. in fact, we assume a maximum-random overlap of the total precip. frac
1136	!
1137	SUBROUTINE poprecip_fracupdate( &
1138	klon, cldfra, precipfracclr, precipfraccld, &
1139	rainclr, raincld, snowclr, snowcld)
1140
1141	USE lmdz_lscp_ini, ONLY : eps
1142
1143	IMPLICIT NONE
1144
1145	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
1146
1147	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction [-]
1148	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
1149	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
1150	!--NB. at the end of the routine, becomes the fraction of precip
1151	!--in the current layer
1152
1153	REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
1154	REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
1155	REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
1156	REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
1157
1158	!--Local variables
1159	INTEGER :: i
1160	REAL :: dcldfra
1161	REAL :: precipfractot
1162	REAL :: dprecipfracclr, dprecipfraccld
1163	REAL :: drainclr, dsnowclr
1164	REAL :: draincld, dsnowcld
1165
1166
1167	DO i = 1, klon
1168
1169	!--Initialisation
1170	precipfractot = precipfracclr(i) + precipfraccld(i)
1171
1172	!--Instead of using the cloud cover which was use in LTP thesis, we use the
1173	!--total precip. fraction to compute the maximum-random overlap. This is
1174	!--because all the information of the cloud cover is embedded into
1175	!--precipfractot, and this allows for taking into account the potential
1176	!--reduction of the precipitation fraction because either the flux is too
1177	!--small (see barrier at the end of poprecip_postcld) or the flux is completely
1178	!--evaporated (see barrier at the end of poprecip_precld)
1179	!--NB. precipfraccld(i) is here the cloud fraction of the layer above
1180	!precipfractot = 1. - ( 1. - precipfractot ) * &
1181	! ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
1182	! / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
1183
1184
1185	IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN
1186	precipfractot = 1.
1187	ELSE
1188	precipfractot = 1. - ( 1. - precipfractot ) * &
1189	( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
1190	/ ( 1. - precipfraccld(i) )
1191	ENDIF
1192
1193	!--precipfraccld(i) is here the cloud fraction of the layer above
1194	dcldfra = cldfra(i) - precipfraccld(i)
1195	!--Tendency of the clear-sky precipitation fraction. We add a MAX on the
1196	!--calculation of the current CS precip. frac.
1197	!dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
1198	!--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated
1199	!--if precipfractot < cldfra)
1200	dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i)
1201	!--Tendency of the cloudy precipitation fraction. We add a MAX on the
1202	!--calculation of the current CS precip. frac.
1203	!dprecipfraccld = MAX( dcldfra , - precipfraccld(i) )
1204	!--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated
1205	!--if cldfra < 0)
1206	dprecipfraccld = dcldfra
1207
1208
1209	!--If the cloud extends
1210	IF ( dprecipfraccld .GT. 0. ) THEN
1211	!--If there is no CS precip, nothing happens.
1212	!--If there is, we reattribute some of the CS precip flux
1213	!--to the cloud precip flux, proportionnally to the
1214	!--decrease of the CS precip fraction
1215	IF ( precipfracclr(i) .LT. eps ) THEN
1216	drainclr = 0.
1217	dsnowclr = 0.
1218	ELSE
1219	drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i)
1220	dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i)
1221	ENDIF
1222	!--If the cloud narrows
1223	ELSEIF ( dprecipfraccld .LT. 0. ) THEN
1224	!--We reattribute some of the cloudy precip flux
1225	!--to the CS precip flux, proportionnally to the
1226	!--decrease of the cloud precip fraction
1227	IF ( precipfraccld(i) .LT. eps ) THEN
1228	draincld = 0.
1229	dsnowcld = 0.
1230	ELSE
1231	draincld = dprecipfraccld / precipfraccld(i) * raincld(i)
1232	dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i)
1233	ENDIF
1234	drainclr = - draincld
1235	dsnowclr = - dsnowcld
1236	!--If the cloud stays the same or if there is no cloud above and
1237	!--in the current layer, nothing happens
1238	ELSE
1239	drainclr = 0.
1240	dsnowclr = 0.
1241	ENDIF
1242
1243	!--We add the tendencies
1244	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1245	precipfraccld(i) = precipfraccld(i) + dprecipfraccld
1246	precipfracclr(i) = precipfracclr(i) + dprecipfracclr
1247	rainclr(i) = MAX(0., rainclr(i) + drainclr)
1248	snowclr(i) = MAX(0., snowclr(i) + dsnowclr)
1249	raincld(i) = MAX(0., raincld(i) - drainclr)
1250	snowcld(i) = MAX(0., snowcld(i) - dsnowclr)
1251
1252	ENDDO
1253
1254	END SUBROUTINE poprecip_fracupdate
1255
1256
1257	!----------------------------------------------------------------
1258	! Computes the processes-oriented precipitation formulations for
1259	! - autoconversion (auto) via a deposition process
1260	! - aggregation (agg)
1261	! - riming (rim)
1262	! - collection (col)
1263	! - melting (melt)
1264	! - freezing (freez)
1265	!
1266	SUBROUTINE poprecip_postcld( &
1267	klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
1268	temp, qvap, qliq, qice, icefrac, cldfra, &
1269	precipfracclr, precipfraccld, &
1270	rain, rainclr, raincld, snow, snowclr, snowcld, &
1271	qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
1272	dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
1273
1274	USE lmdz_lscp_ini, ONLY : prt_level, lunout
1275	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
1276	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
1277
1278	USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb, &
1279	cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, &
1280	thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
1281	rho_rain, r_rain, r_snow, rho_ice, &
1282	tau_auto_snow_min, tau_auto_snow_max, &
1283	thresh_precip_frac, eps, &
1284	gamma_melt, alpha_freez, beta_freez, temp_nowater, &
1285	iflag_cloudth_vert, iflag_rain_incloud_vol, &
1286	cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez, &
1287	rain_fallspeed_clr, rain_fallspeed_cld, &
1288	snow_fallspeed_clr, snow_fallspeed_cld
1289
1290
1291	IMPLICIT NONE
1292
1293	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
1294	REAL, INTENT(IN) :: dtime !--time step [s]
1295
1296	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
1297	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
1298	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
1299
1300	REAL, INTENT(IN), DIMENSION(klon) :: ctot_vol !--volumic cloud fraction [-]
1301	LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv !--true if we are in a convective point
1302
1303	REAL, INTENT(INOUT), DIMENSION(klon) :: temp !--current temperature [K]
1304	REAL, INTENT(INOUT), DIMENSION(klon) :: qvap !--current water vapor specific humidity [kg/kg]
1305	REAL, INTENT(INOUT), DIMENSION(klon) :: qliq !--current liquid water specific humidity [kg/kg]
1306	REAL, INTENT(INOUT), DIMENSION(klon) :: qice !--current ice water specific humidity [kg/kg]
1307	REAL, INTENT(IN), DIMENSION(klon) :: icefrac !--ice fraction [-]
1308	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction [-]
1309
1310	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
1311	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
1312	!--NB. at the end of the routine, becomes the fraction of precip
1313	!--in the current layer
1314
1315	REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
1316	REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
1317	REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
1318	REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
1319	REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
1320	REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
1321
1322	REAL, INTENT(OUT), DIMENSION(klon) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg]
1323	REAL, INTENT(OUT), DIMENSION(klon) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg]
1324	REAL, INTENT(OUT), DIMENSION(klon) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
1325	REAL, INTENT(OUT), DIMENSION(klon) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
1326	REAL, INTENT(OUT), DIMENSION(klon) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
1327	REAL, INTENT(OUT), DIMENSION(klon) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
1328	REAL, INTENT(OUT), DIMENSION(klon) :: dqsrim !--snow tendency due to riming [kg/kg/s]
1329	REAL, INTENT(OUT), DIMENSION(klon) :: dqsmelt !--snow tendency due to melting [kg/kg/s]
1330	REAL, INTENT(OUT), DIMENSION(klon) :: dqrmelt !--rain tendency due to melting [kg/kg/s]
1331	REAL, INTENT(OUT), DIMENSION(klon) :: dqsfreez !--snow tendency due to freezing [kg/kg/s]
1332	REAL, INTENT(OUT), DIMENSION(klon) :: dqrfreez !--rain tendency due to freezing [kg/kg/s]
1333
1334
1335
1336	!--Local variables
1337
1338	INTEGER :: i
1339	!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
1340	REAL, DIMENSION(klon) :: dhum_to_dflux
1341	REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip
1342
1343	!--Collection, aggregation and riming
1344	REAL :: eff_cldfra
1345	REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
1346	REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
1347	REAL :: rho_snow
1348	REAL :: dqlcol !--loss of liquid cloud content due to collection by rain [kg/kg/s]
1349	REAL :: dqiagg !--loss of ice cloud content due to collection by aggregation [kg/kg/s]
1350	REAL :: dqlrim !--loss of liquid cloud content due to riming on snow [kg/kg/s]
1351
1352	!--Autoconversion
1353	REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
1354	REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
1355	REAL :: dqlauto !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
1356	REAL :: dqiauto !--loss of ice cloud content due to autoconversion to snow [kg/kg/s]
1357
1358	!--Melting
1359	REAL :: dqsmelt_max, air_thermal_conduct
1360	REAL :: nb_snowflake_clr, nb_snowflake_cld
1361	REAL :: capa_snowflake, temp_wetbulb
1362	REAL :: rho, r_ice
1363	REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
1364	REAL, DIMENSION(klon) :: qzero, qsat, dqsat
1365
1366	!--Freezing
1367	REAL :: dqrfreez_max
1368	REAL :: tau_freez
1369	REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2
1370	REAL :: coef_freez
1371	REAL :: dqifreez !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
1372	REAL :: Eff_rain_ice
1373
1374
1375	!--Initialisation of variables
1376
1377
1378	qzero(:) = 0.
1379
1380	dqrcol(:) = 0.
1381	dqsagg(:) = 0.
1382	dqsauto(:) = 0.
1383	dqrauto(:) = 0.
1384	dqsrim(:) = 0.
1385	dqrmelt(:) = 0.
1386	dqsmelt(:) = 0.
1387	dqrfreez(:) = 0.
1388	dqsfreez(:) = 0.
1389
1390
1391	!--Update the precipitation fraction following cloud formation
1392	CALL poprecip_fracupdate( &
1393	klon, cldfra, precipfracclr, precipfraccld, &
1394	rainclr, raincld, snowclr, snowcld)
1395
1396
1397	DO i = 1, klon
1398
1399	!--Variables initialisation
1400	dqlcol = 0.
1401	dqiagg = 0.
1402	dqiauto = 0.
1403	dqlauto = 0.
1404	dqlrim = 0.
1405
1406	!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
1407	dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
1408	qtot(i) = qvap(i) + qliq(i) + qice(i) &
1409	+ ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i)
1410
1411	!--If vertical heterogeneity is taken into account, we use
1412	!--the "true" volume fraction instead of a modified
1413	!--surface fraction (which is larger and artificially
1414	!--reduces the in-cloud water).
1415	IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
1416	eff_cldfra = ctot_vol(i)
1417	ELSE
1418	eff_cldfra = cldfra(i)
1419	ENDIF
1420
1421
1422	!--Start precipitation growth processes
1423
1424	!--If the cloud is big enough, the precipitation processes activate
1425	! TODO met on seuil_neb ici ?
1426	IF ( cldfra(i) .GE. seuil_neb ) THEN
1427
1428	!---------------------------------------------------------
1429	!-- COLLECTION AND AGGREGATION
1430	!---------------------------------------------------------
1431	!--Collection: processus through which rain collects small liquid droplets
1432	!--in suspension, and add it to the rain flux
1433	!--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
1434	!--Those processes are treated before autoconversion because we do not
1435	!--want to collect/aggregate the newly formed fluxes, which already
1436	!--"saw" the cloud as they come from it
1437	!--The formulas come from Muench and Lohmann 2020
1438
1439	!--gamma_col: tuning coefficient [-]
1440	!--rho_rain: volumic mass of rain [kg/m3]
1441	!--r_rain: size of the rain droplets [m]
1442	!--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
1443	!--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
1444	!--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
1445	!--get in-cloud mean quantities. The two divisions by eff_cldfra are
1446	!--then simplified.
1447
1448	!--The collection efficiency is perfect.
1449	Eff_rain_liq = 1.
1450	coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
1451	IF ( raincld(i) .GT. 0. ) THEN
1452	!--Exact explicit version, which does not need a barrier because of
1453	!--the exponential decrease
1454	dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
1455
1456	!--Add tendencies
1457	qliq(i) = qliq(i) + dqlcol
1458	raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i)
1459
1460	!--Diagnostic tendencies
1461	dqrcol(i) = - dqlcol / dtime
1462	ENDIF
1463
1464	!--Same as for aggregation
1465	!--Eff_snow_liq formula:
1466	!--it s a product of a collection efficiency and a sticking efficiency
1467	! Milbrandt and Yau formula that gives very low values:
1468	! Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
1469	! Lin 1983's formula
1470	Eff_snow_ice = EXP( 0.025 * MIN( ( temp(i) - RTT ), 0.) )
1471	!--rho_snow formula follows Brandes et al. 2007 (JAMC)
1472	rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
1473	coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
1474	IF ( snowcld(i) .GT. 0. ) THEN
1475	!--Exact explicit version, which does not need a barrier because of
1476	!--the exponential decrease
1477	dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
1478
1479	!--Add tendencies
1480	qice(i) = qice(i) + dqiagg
1481	snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i)
1482
1483	!--Diagnostic tendencies
1484	dqsagg(i) = - dqiagg / dtime
1485	ENDIF
1486
1487
1488	!---------------------------------------------------------
1489	!-- AUTOCONVERSION
1490	!---------------------------------------------------------
1491	!--Autoconversion converts liquid droplets/ice crystals into
1492	!--rain drops/snowflakes. It relies on the formulations by
1493	!--Sundqvist 1978.
1494
1495	!--If we are in a convective point, we have different parameters
1496	!--for the autoconversion
1497	IF ( ptconv(i) ) THEN
1498	qthresh_auto_rain = cld_lc_con
1499	qthresh_auto_snow = cld_lc_con_snow
1500
1501	tau_auto_rain = cld_tau_con
1502	!--tau for snow depends on the ice fraction in mixed-phase clouds
1503	tau_auto_snow = tau_auto_snow_max &
1504	+ ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
1505
1506	expo_auto_rain = cld_expo_con
1507	expo_auto_snow = cld_expo_con
1508	ELSE
1509	qthresh_auto_rain = cld_lc_lsc
1510	qthresh_auto_snow = cld_lc_lsc_snow
1511
1512	tau_auto_rain = cld_tau_lsc
1513	!--tau for snow depends on the ice fraction in mixed-phase clouds
1514	tau_auto_snow = tau_auto_snow_max &
1515	+ ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
1516
1517	expo_auto_rain = cld_expo_lsc
1518	expo_auto_snow = cld_expo_lsc
1519	ENDIF
1520
1521
1522	! Liquid water quantity to remove according to (Sundqvist, 1978)
1523	! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) )
1524	!
1525	!--And same formula for ice
1526	!
1527	!--We first treat the second term (with exponential) in an explicit way
1528	!--and then treat the first term (-q/tau) in an exact way
1529
1530	dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
1531	- ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
1532
1533	dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
1534	- ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
1535
1536
1537	!--Barriers so that we don't create more rain/snow
1538	!--than there is liquid/ice
1539	dqlauto = MAX( - qliq(i), dqlauto )
1540	dqiauto = MAX( - qice(i), dqiauto )
1541
1542	!--Add tendencies
1543	qliq(i) = qliq(i) + dqlauto
1544	qice(i) = qice(i) + dqiauto
1545	raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
1546	snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
1547
1548	!--Diagnostic tendencies
1549	dqsauto(i) = - dqiauto / dtime
1550	dqrauto(i) = - dqlauto / dtime
1551
1552
1553	!---------------------------------------------------------
1554	!-- RIMING
1555	!---------------------------------------------------------
1556	!--Process which converts liquid droplets in suspension into
1557	!--snow because of the collision between
1558	!--those and falling snowflakes.
1559	!--The formula comes from Muench and Lohmann 2020
1560	!--NB.: this process needs a temperature adjustment
1561
1562	!--Eff_snow_liq formula: following Ferrier 1994,
1563	!--assuming 1
1564	Eff_snow_liq = 1.0
1565	!--rho_snow formula follows Brandes et al. 2007 (JAMC)
1566	rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
1567	coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
1568	IF ( snowcld(i) .GT. 0. ) THEN
1569	!--Exact version, which does not need a barrier because of
1570	!--the exponential decrease
1571	dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. )
1572
1573	!--Add tendencies
1574	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1575	qliq(i) = qliq(i) + dqlrim
1576	snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i)
1577
1578	!--Temperature adjustment with the release of latent
1579	!--heat because of solid condensation
1580	temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
1581	/ ( 1. + RVTMP2 * qtot(i) )
1582
1583	!--Diagnostic tendencies
1584	dqsrim(i) = - dqlrim / dtime
1585	ENDIF
1586
1587	ENDIF ! cldfra .GE. seuil_neb
1588
1589	ENDDO ! loop on klon
1590
1591
1592	!--Re-calculation of saturation specific humidity
1593	!--because riming changed temperature
1594	CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
1595
1596	DO i = 1, klon
1597
1598	!---------------------------------------------------------
1599	!-- MELTING
1600	!---------------------------------------------------------
1601	!--Process through which snow melts into rain.
1602	!--The formula is homemade.
1603	!--NB.: this process needs a temperature adjustment
1604
1605	!--dqsmelt_max : maximum snow melting so that temperature
1606	!-- stays higher than 273 K [kg/kg]
1607	!--capa_snowflake : capacitance of a snowflake, equal to
1608	!-- the radius if the snowflake is a sphere [m]
1609	!--temp_wetbulb : wet-bulb temperature [K]
1610	!--snow_fallspeed : snow fall velocity (in clear/cloudy sky) [m/s]
1611	!--air_thermal_conduct : thermal conductivity of the air [J/m/K/s]
1612	!--gamma_melt : tuning parameter for melting [-]
1613	!--nb_snowflake : number of snowflakes (in clear/cloudy air) [-]
1614
1615	IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
1616	!--Computed according to
1617	!--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1618	dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
1619	* ( 1. + RVTMP2 * qtot(i) ))
1620
1621	!--Initialisation
1622	dqsclrmelt = 0.
1623	dqscldmelt = 0.
1624
1625	!--We assume that the snowflakes are spherical
1626	capa_snowflake = r_snow
1627	!--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
1628	air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
1629	!--rho_snow formula follows Brandes et al. 2007 (JAMC)
1630	rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
1631
1632	!--In clear air
1633	IF ( ( snowclr(i) .GT. 0. ) .AND. ( precipfracclr(i) .GT. 0. ) ) THEN
1634	!--Formula for the wet-bulb temperature from ECMWF (IFS)
1635	!--The vapor used is the vapor in the clear sky
1636	temp_wetbulb = temp(i) &
1637	- ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) &
1638	* ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
1639	- 40.637 * ( temp(i) - 275. ) )
1640	!--Calculated according to
1641	!-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
1642	nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr &
1643	/ ( 4. / 3. * RPI * r_snow*3. rho_snow )
1644	dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
1645	* capa_snowflake / RLMLT * gamma_melt &
1646	* MAX(0., temp_wetbulb - RTT) * dtime
1647
1648	!--Barrier to limit the melting flux to the clr snow flux in the mesh
1649	dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i))
1650	ENDIF
1651
1652
1653	!--In cloudy air
1654	IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
1655	!--As the air is saturated, the wet-bulb temperature is equal to the
1656	!--temperature
1657	temp_wetbulb = temp(i)
1658	!--Calculated according to
1659	!-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
1660	nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld &
1661	/ ( 4. / 3. * RPI * r_snow*3. rho_snow )
1662	dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
1663	* capa_snowflake / RLMLT * gamma_melt &
1664	* MAX(0., temp_wetbulb - RTT) * dtime
1665
1666	!--Barrier to limit the melting flux to the cld snow flux in the mesh
1667	dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i))
1668	ENDIF
1669
1670
1671	!--Barrier on temperature. If the total melting flux leads to a
1672	!--positive temperature, it is limited to keep temperature above 0 degC.
1673	!--It is activated if the total is LOWER than the max
1674	!--because everything is negative
1675	dqstotmelt = dqsclrmelt + dqscldmelt
1676	IF ( dqstotmelt .LT. dqsmelt_max ) THEN
1677	!--We redistribute the max melted snow keeping
1678	!--the clear/cloud partition of the melted snow
1679	dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
1680	dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
1681	dqstotmelt = dqsmelt_max
1682
1683	ENDIF
1684
1685	!--Add tendencies
1686	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1687	rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i))
1688	raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i))
1689	snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i))
1690	snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i))
1691
1692	!--Temperature adjustment with the release of latent
1693	!--heat because of melting
1694	temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
1695	/ ( 1. + RVTMP2 * qtot(i) )
1696
1697	!--Diagnostic tendencies
1698	dqrmelt(i) = - dqstotmelt / dtime
1699	dqsmelt(i) = dqstotmelt / dtime
1700
1701	ENDIF
1702
1703
1704	!---------------------------------------------------------
1705	!-- FREEZING
1706	!---------------------------------------------------------
1707	!--Process through which rain freezes into snow.
1708	!-- We parameterize it as a 2 step process:
1709	!--first: freezing following collision with ice crystals
1710	!--second: immersion freezing following (inspired by Bigg 1953)
1711	!--the latter is parameterized as an exponential decrease of the rain
1712	!--water content with a homemade formulya
1713	!--This is based on a caracteritic time of freezing, which
1714	!--exponentially depends on temperature so that it is
1715	!--equal to 1 for temp_nowater (see below) and is close to
1716	!--0 for RTT (=273.15 K).
1717	!--NB.: this process needs a temperature adjustment
1718	!--dqrfreez_max : maximum rain freezing so that temperature
1719	!-- stays lower than 273 K [kg/kg]
1720	!--tau_freez : caracteristic time of freezing [s]
1721	!--gamma_freez : tuning parameter [s-1]
1722	!--alpha_freez : tuning parameter for the shape of the exponential curve [-]
1723	!--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC)
1724
1725	IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
1726
1727
1728	!--1st step: freezing following collision with ice crystals
1729	!--Sub-process of freezing which quantifies the collision between
1730	!--ice crystals in suspension and falling rain droplets.
1731	!--The rain droplets freeze, becoming graupel, and carrying
1732	!--the ice crystal (which acted as an ice nucleating particle).
1733	!--The formula is adapted from the riming formula.
1734	!--it works only in the cloudy part
1735
1736	dqifreez = 0.
1737	dqrtotfreez_step1 = 0.
1738
1739	IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. &
1740	( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
1741	dqrclrfreez = 0.
1742	dqrcldfreez = 0.
1743
1744	!--Computed according to
1745	!--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1746	dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
1747	* ( 1. + RVTMP2 * qtot(i) ))
1748
1749
1750	!--The collision efficiency is assumed unity
1751	Eff_rain_ice = 1.
1752	coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
1753	!--Exact version, which does not need a barrier because of
1754	!--the exponential decrease.
1755	dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
1756
1757	!--We add the part of rain water that freezes, limited by a temperature barrier
1758	!--This quantity is calculated assuming that the number of drop that freeze correspond to the number
1759	!--of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
1760	!--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3
1761	!--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3
1762	!--The assumption above corresponds to dNi = dNr, i.e.,
1763	!-- dqr = dqi * (4/3 PI rho_rain * r_rain*3) / (4/3 PI rho_ice r_ice**3)
1764	!--Dry density [kg/m3]
1765	rho = pplay(i) / temp(i) / RD
1766	!--r_ice formula from Sun and Rikus (1999)
1767	r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 &
1768	+ 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )*0.2535 ( temp(i) - RTT + 190. ) ) / 2.
1769	dqrcldfreez = dqifreez * rho_rain * r_rain*3. / ( rho_ice r_ice**3. )
1770	dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i))
1771	dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max)
1772	dqrtotfreez_step1 = dqrcldfreez
1773
1774	!--Add tendencies
1775	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1776	qice(i) = qice(i) + dqifreez
1777	raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
1778	snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i))
1779	temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD &
1780	/ ( 1. + RVTMP2 * qtot(i) )
1781
1782	ENDIF
1783
1784	!-- Second step immersion freezing of rain drops
1785	!-- with a homemade timeconstant depending on temperature
1786
1787	dqrclrfreez = 0.
1788	dqrcldfreez = 0.
1789	dqrtotfreez_step2 = 0.
1790	!--Computed according to
1791	!--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1792
1793	dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
1794	* ( 1. + RVTMP2 * qtot(i) ))
1795
1796
1797	tau_freez = 1. / ( beta_freez &
1798	* EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
1799
1800
1801	!--In clear air
1802	IF ( rainclr(i) .GT. 0. ) THEN
1803	!--Exact solution of dqrain/dt = -qrain/tau_freez
1804	dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
1805	ENDIF
1806
1807	!--In cloudy air
1808	IF ( raincld(i) .GT. 0. ) THEN
1809	!--Exact solution of dqrain/dt = -qrain/tau_freez
1810	dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
1811	ENDIF
1812
1813	!--temperature barrier step 2
1814	!--It is activated if the total is LOWER than the max
1815	!--because everything is negative
1816	dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
1817
1818	IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
1819	!--We redistribute the max freezed rain keeping
1820	!--the clear/cloud partition of the freezing rain
1821	dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
1822	dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
1823	dqrtotfreez_step2 = dqrfreez_max
1824	ENDIF
1825
1826
1827	!--Add tendencies
1828	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1829	rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i))
1830	raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
1831	snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i))
1832	snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i))
1833
1834
1835	!--Temperature adjustment with the uptake of latent
1836	!--heat because of freezing
1837	temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
1838	/ ( 1. + RVTMP2 * qtot(i) )
1839
1840	!--Diagnostic tendencies
1841	dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2
1842	dqrfreez(i) = dqrtotfreez / dtime
1843	dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
1844
1845	ENDIF
1846
1847
1848
1849	!--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
1850	!--we reduce the precipiration fraction in the clear/cloudy air so that the new
1851	!--local flux of rain+snow is equal to rain_int_min.
1852	!--Here, rain+snow is the gridbox-mean flux of precip.
1853	!--Therefore, (rain+snow)/precipfrac is the local flux of precip.
1854	!--If the local flux of precip is lower than rain_int_min, i.e.,
1855	!-- (rain+snow)/precipfrac < rain_int_min , i.e.,
1856	!-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
1857	!--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
1858	!--Note that this is physically different than what is proposed in LTP thesis.
1859	precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
1860	precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
1861
1862	!--Calculate outputs
1863	rain(i) = rainclr(i) + raincld(i)
1864	snow(i) = snowclr(i) + snowcld(i)
1865
1866	!--Diagnostics
1867	!--BEWARE this is indeed a diagnostic: this is an estimation from
1868	!--the value of the flux at the bottom interface of the mesh and
1869	!--and assuming an upstream numerical calculation
1870	!--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is
1871	!--used for computing the total ice water content in the mesh
1872	!--for radiation only
1873	qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr &
1874	+ raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld )
1875	qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr &
1876	+ snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld )
1877
1878
1879	ENDDO ! loop on klon
1880
1881	END SUBROUTINE poprecip_postcld
1882
1883	END MODULE lmdz_lscp_precip

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format