Context Navigation

lmdz_lscp_precip.f90 @ 5624

Last change on this file since 5624 was 5613, checked in by aborella, 2 months ago
Added barriers and bypass overflows
File size: 84.9 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, icesed_flux, &
21	zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, dqreva, dqssub &
22	)
23
24	USE lmdz_lscp_ini, ONLY : iflag_evap_prec, ok_ice_sedim
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	REAL, INTENT(IN), DIMENSION(klon) :: icesed_flux !--sedimentated ice flux [kg/s/m2]
47
48	REAL, INTENT(IN), DIMENSION(klon) :: zneb !--cloud fraction IN THE LAYER ABOVE [-]
49	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecip !--fraction of precipitation IN THE LAYER ABOVE [-]
50	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
51
52	REAL, INTENT(INOUT), DIMENSION(klon) :: zrfl !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
53	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
54	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflcld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
55	REAL, INTENT(INOUT), DIMENSION(klon) :: zifl !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
56	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
57	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
58
59	REAL, INTENT(OUT), DIMENSION(klon) :: dqreva !--rain tendency due to evaporation [kg/kg/s]
60	REAL, INTENT(OUT), DIMENSION(klon) :: dqssub !--snow tendency due to sublimation [kg/kg/s]
61
62
63	REAL :: zmair
64	REAL :: zcpair, zcpeau
65	REAL, DIMENSION(klon) :: qzero
66	REAL, DIMENSION(klon) :: zqs, zdqs
67	REAL, DIMENSION(klon) :: qsl, qsi
68	REAL, DIMENSION(klon) :: dqsl, dqsi
69
70	REAL :: zqev0
71	REAL :: zqev, zqevt
72	REAL :: zqevi, zqevti
73
74	REAL, DIMENSION(klon) :: zrfln, zifln
75
76	REAL :: zmelt
77	REAL :: qice_sedim
78	INTEGER :: i
79
80	qzero(:) = 0.
81
82
83	! --------------------------------------------------------------------
84	! P1> Thermalization of precipitation falling from the overlying layer
85	! --------------------------------------------------------------------
86	! Computes air temperature variation due to enthalpy transported by
87	! precipitation. Precipitation is then thermalized with the air in the
88	! layer.
89	! The precipitation should remain thermalized throughout the different
90	! thermodynamical transformations.
91	! The corresponding water mass should
92	! be added when calculating the layer's enthalpy change with
93	! temperature
94	! ---------------------------------------------------------------------
95
96	IF (iftop) THEN
97
98	DO i = 1, klon
99	zmqc(i) = 0.
100	ENDDO
101
102	ELSE
103
104	DO i = 1, klon
105
106	zmair=(paprsdn(i)-paprsup(i))/RG
107	! no condensed water so cp=cp(vapor+dry air)
108	! RVTMP2=rcpv/rcpd-1
109	zcpair=RCPD(1.0+RVTMP2zq(i))
110	zcpeau=RCPD*RVTMP2
111
112	! zmqc: precipitation mass that has to be thermalized with
113	! layer's air so that precipitation at the ground has the
114	! same temperature as the lowermost layer
115	zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
116	! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
117	zt(i) = ( ztupnew(i)zmqc(i)zcpeau + zcpair*zt(i) ) &
118	/ (zcpair + zmqc(i)*zcpeau)
119
120	ENDDO
121
122	ENDIF
123
124	!--------------------
125	!-- ICE SEDIMENTATION
126	!--------------------
127	IF ( ok_ice_sedim ) THEN
128	DO i =1, klon
129	!--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
130	!--added to the total water content of the gridbox
131	IF ( icesed_flux(i) .GT. 0. ) THEN
132	qice_sedim = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
133	!--Vapor is updated after evaporation/sublimation (it is increased)
134	zq(i) = zq(i) + qice_sedim
135	!--Air and precip temperature (i.e., gridbox temperature)
136	!--is updated due to latent heat cooling
137	zt(i) = zt(i) &
138	- qice_sedim * RLSTT / RCPD &
139	/ ( 1. + RVTMP2 * ( zq(i) + zmqc(i) ) )
140	ENDIF ! ok_ice_sedim
141	ENDDO
142	ENDIF
143
144	! --------------------------------------------------------------------
145	! P2> Precipitation evaporation/sublimation/melting
146	! --------------------------------------------------------------------
147	! A part of the precipitation coming from above is evaporated/sublimated/melted.
148	! --------------------------------------------------------------------
149
150	IF (iflag_evap_prec.GE.1) THEN
151
152	! Calculation of saturation specific humidity
153	! depending on temperature:
154	CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,0,.false.,zqs,zdqs)
155	! wrt liquid water
156	CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,1,.false.,qsl,dqsl)
157	! wrt ice
158	CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,2,.false.,qsi,dqsi)
159
160	DO i = 1, klon
161
162	! if precipitation
163	IF (zrfl(i)+zifl(i).GT.0.) THEN
164
165	! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
166	! c_iso: likely important to distinguish cs from neb isotope precipitation
167
168	IF (iflag_evap_prec.GE.4) THEN
169	zrfl(i) = zrflclr(i)
170	zifl(i) = ziflclr(i)
171	ENDIF
172
173	IF (iflag_evap_prec.EQ.1) THEN
174	znebprecip(i)=zneb(i)
175	ELSE
176	znebprecip(i)=MAX(zneb(i),znebprecip(i))
177	ENDIF
178
179	IF (iflag_evap_prec.GT.4) THEN
180	! Max evaporation not to saturate the clear sky precip fraction
181	! i.e. the fraction where evaporation occurs
182	zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
183	ELSEIF (iflag_evap_prec .EQ. 4) THEN
184	! Max evaporation not to saturate the whole mesh
185	! Pay attention -> lead to unrealistic and excessive evaporation
186	zqev0 = MAX(0.0, zqs(i)-zq(i))
187	ELSE
188	! Max evap not to saturate the fraction below the cloud
189	zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
190	ENDIF
191
192	! Evaporation of liquid precipitation coming from above
193	! dP/dz=beta(1-q/qsat)sqrt(P)
194	! formula from Sundquist 1988, Klemp & Wilhemson 1978
195	! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
196
197	IF (iflag_evap_prec.EQ.3) THEN
198	zqevt = znebprecip(i)coef_eva(1.0-zq(i)/qsl(i)) &
199	*SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
200	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
201	ELSE IF (iflag_evap_prec.GE.4) THEN
202	zqevt = znebprecipclr(i)coef_eva(1.0-zq(i)/qsl(i)) &
203	*SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
204	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
205	ELSE
206	zqevt = 1.coef_eva(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
207	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
208	ENDIF
209
210	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * RG*dtime/(paprsdn(i)-paprsup(i))
211
212	! sublimation of the solid precipitation coming from above
213	IF (iflag_evap_prec.EQ.3) THEN
214	zqevti = znebprecip(i)coef_sub(1.0-zq(i)/qsi(i)) &
215	*SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
216	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
217	ELSE IF (iflag_evap_prec.GE.4) THEN
218	zqevti = znebprecipclr(i)coef_sub(1.0-zq(i)/qsi(i)) &
219	*SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
220	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
221	ELSE
222	zqevti = 1.coef_sub(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
223	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
224	ENDIF
225
226	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) * RG*dtime/(paprsdn(i)-paprsup(i))
227
228	! A. JAM
229	! Evaporation limit: we ensure that the layer's fraction below
230	! the cloud or the whole mesh (depending on iflag_evap_prec)
231	! does not reach saturation. In this case, we
232	! redistribute zqev0 conserving the ratio liquid/ice
233
234	IF (zqevt+zqevti.GT.zqev0) THEN
235	zqev=zqev0*zqevt/(zqevt+zqevti)
236	zqevi=zqev0*zqevti/(zqevt+zqevti)
237	ELSE
238	zqev=zqevt
239	zqevi=zqevti
240	ENDIF
241
242	!--Diagnostics
243	dqreva(i) = - zqev / dtime
244	dqssub(i) = - zqevti / dtime
245
246	! New solid and liquid precipitation fluxes after evap and sublimation
247	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprsdn(i)-paprsup(i))/RG/dtime)
248	zifln(i) = Max(0.,zifl(i) - zqevi*(paprsdn(i)-paprsup(i))/RG/dtime)
249
250
251	! vapor, temperature, precip fluxes update
252	! vapor is updated after evaporation/sublimation (it is increased)
253	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
254	* (RG/(paprsdn(i)-paprsup(i)))*dtime
255	! zmqc is the total condensed water in the precip flux (it is decreased)
256	zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
257	* (RG/(paprsdn(i)-paprsup(i)))*dtime
258	! air and precip temperature (i.e., gridbox temperature)
259	! is updated due to latent heat cooling
260	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
261	* (RG/(paprsdn(i)-paprsup(i)))*dtime &
262	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
263	+ (zifln(i)-zifl(i)) &
264	* (RG/(paprsdn(i)-paprsup(i)))*dtime &
265	* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
266
267	! New values of liquid and solid precipitation
268	zrfl(i) = zrfln(i)
269	zifl(i) = zifln(i)
270
271	! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
272	! due to evap + sublim
273
274
275	IF (iflag_evap_prec.GE.4) THEN
276	zrflclr(i) = zrfl(i)
277	ziflclr(i) = zifl(i)
278	IF(zrflclr(i) + ziflclr(i).LE.0) THEN
279	znebprecipclr(i) = 0.0
280	ENDIF
281	zrfl(i) = zrflclr(i) + zrflcld(i)
282	zifl(i) = ziflclr(i) + ziflcld(i)
283	ENDIF
284
285	! c_iso duplicate for isotopes or loop on isotopes
286
287	! Melting:
288	zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
289	! precip fraction that is melted
290	zmelt = MIN(MAX(zmelt,0.),1.)
291
292	! update of rainfall and snowfall due to melting
293	IF (iflag_evap_prec.GE.4) THEN
294	zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
295	zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
296	zrfl(i)=zrflclr(i)+zrflcld(i)
297	ELSE
298	zrfl(i)=zrfl(i)+zmelt*zifl(i)
299	ENDIF
300
301
302	! c_iso: melting of isotopic precipi with zmelt (no fractionation)
303
304	! Latent heat of melting because of precipitation melting
305	! NB: the air + precip temperature is simultaneously updated
306	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprsdn(i)-paprsup(i)) &
307	RLMLT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
308
309	IF (iflag_evap_prec.GE.4) THEN
310	ziflclr(i)=ziflclr(i)*(1.-zmelt)
311	ziflcld(i)=ziflcld(i)*(1.-zmelt)
312	zifl(i)=ziflclr(i)+ziflcld(i)
313	ELSE
314	zifl(i)=zifl(i)*(1.-zmelt)
315	ENDIF
316
317	ELSE
318	! if no precip, we reinitialize the cloud fraction used for the precip to 0
319	znebprecip(i)=0.
320
321	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
322
323	ENDDO ! loop on klon
324
325	ENDIF ! (iflag_evap_prec>=1)
326
327	END SUBROUTINE histprecip_precld
328
329
330	SUBROUTINE histprecip_postcld( &
331	klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, pt_pron_clds, &
332	zdqsdT_raw, zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
333	rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, &
334	zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
335	zradocond, zradoice, dqrauto, dqsauto, dqised &
336	)
337
338	USE lmdz_lscp_ini, ONLY : RD, RG, RCPD, RVTMP2, RLSTT, RLMLT, RTT
339	USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, ffallv_con, &
340	cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, ffallv_lsc, &
341	seuil_neb, rain_int_min, cice_velo, dice_velo
342	USE lmdz_lscp_ini, ONLY : iflag_evap_prec, iflag_cloudth_vert, iflag_rain_incloud_vol, &
343	iflag_autoconversion, ok_radocond_snow, ok_bug_phase_lscp, &
344	niter_lscp
345	USE lmdz_lscp_ini, ONLY : ok_ice_sedim, fallice_sedim, eps
346	USE lmdz_lscp_tools, ONLY : fallice_velocity
347
348	IMPLICIT NONE
349
350
351	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
352	REAL, INTENT(IN) :: dtime !--time step [s]
353	LOGICAL, INTENT(IN) :: iftop !--if top of the column
354
355	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
356	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
357	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
358	REAL, INTENT(IN), DIMENSION(klon) :: ctot_vol !--volumic cloud fraction [-]
359	LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv !--true if we are in a convective point
360	LOGICAL, INTENT(IN), DIMENSION(klon) :: pt_pron_clds !--true if we are in a prognostic clouds point
361	REAL, INTENT(IN), DIMENSION(klon) :: zdqsdT_raw !--derivative of qsat w.r.t. temperature times L/Cp [SI]
362
363	REAL, INTENT(INOUT), DIMENSION(klon) :: zt !--current temperature [K]
364	REAL, INTENT(INOUT), DIMENSION(klon) :: zq !--current water vapor specific humidity [kg/kg]
365	REAL, INTENT(INOUT), DIMENSION(klon) :: zoliq !--current liquid+ice water specific humidity [kg/kg]
366	REAL, INTENT(INOUT), DIMENSION(klon) :: zoliql !--current liquid water specific humidity [kg/kg]
367	REAL, INTENT(INOUT), DIMENSION(klon) :: zoliqi !--current ice water specific humidity [kg/kg]
368	REAL, INTENT(INOUT), DIMENSION(klon) :: zcond !--liquid+ice water specific humidity AFTER CONDENSATION (no precip process) [kg/kg]
369	REAL, INTENT(IN), DIMENSION(klon) :: zfice !--ice fraction AFTER CONDENSATION [-]
370	REAL, INTENT(IN), DIMENSION(klon) :: zmqc !--specific humidity in the precipitation falling from the upper layer [kg/kg]
371
372	REAL, INTENT(IN), DIMENSION(klon) :: rneb !--cloud fraction [-]
373	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
374	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipcld !--fraction of precipitation in the cloud IN THE LAYER ABOVE [-]
375	REAL, INTENT(INOUT), DIMENSION(klon) :: zneb !--cloud fraction IN THE LAYER ABOVE [-]
376	REAL, INTENT(INOUT), DIMENSION(klon) :: tot_zneb !--total cloud cover above the current mesh [-]
377	REAL, INTENT(INOUT), DIMENSION(klon) :: zrho_up !--air density IN THE LAYER ABOVE [kg/m3]
378	REAL, INTENT(INOUT), DIMENSION(klon) :: zvelo_up !--ice fallspeed velocity IN THE LAYER ABOVE [m/s]
379
380	REAL, INTENT(INOUT), DIMENSION(klon) :: zrfl !--flux of rain gridbox-mean [kg/s/m2]
381	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflclr !--flux of rain gridbox-mean in clear sky [kg/s/m2]
382	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflcld !--flux of rain gridbox-mean in cloudy air [kg/s/m2]
383	REAL, INTENT(INOUT), DIMENSION(klon) :: zifl !--flux of snow gridbox-mean [kg/s/m2]
384	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflclr !--flux of snow gridbox-mean in clear sky [kg/s/m2]
385	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflcld !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
386	REAL, INTENT(OUT), DIMENSION(klon) :: icesed_flux !--flux of sedimentated ice [kg/s/m2]
387
388	REAL, INTENT(OUT), DIMENSION(klon) :: zradocond !--condensed water used in the radiation scheme [kg/kg]
389	REAL, INTENT(OUT), DIMENSION(klon) :: zradoice !--condensed ice water used in the radiation scheme [kg/kg]
390	REAL, INTENT(OUT), DIMENSION(klon) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
391	REAL, INTENT(OUT), DIMENSION(klon) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
392	REAL, INTENT(INOUT), DIMENSION(klon) :: dqised !--ice water content tendency due to sedimentation of ice crystals [kg/kg/s]
393
394
395	! Local variables for precip fraction update
396	REAL :: smallestreal
397	REAL, DIMENSION(klon) :: tot_znebn, d_tot_zneb
398	REAL, DIMENSION(klon) :: d_znebprecip_cld_clr, d_znebprecip_clr_cld
399	REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
400	REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
401
402	! Local variables for autoconversion
403	REAL :: zct, zcl, zexpo, ffallv
404	REAL :: zchau, zfroi
405	REAL :: zrain, zsnow, zprecip
406	REAL :: effective_zneb
407	REAL, DIMENSION(klon) :: zrho, zvelo
408	REAL, DIMENSION(klon) :: zdz, iwc
409
410	! Local variables for Bergeron process
411	REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
412	REAL, DIMENSION(klon) :: zqpreci, zqprecl
413
414	! Local variables for calculation of quantity of condensates seen by the radiative scheme
415	REAL, DIMENSION(klon) :: zradoliq
416	REAL, DIMENSION(klon) :: ziflprev
417
418	! Local variable for sedimentation
419	REAL :: iwcg, tempc, icevelo
420	REAL :: qice_sedim
421
422	! Misc
423	INTEGER :: i, n
424
425	! Initialisation
426	smallestreal=1.e-9
427	zqprecl(:) = 0.
428	zqpreci(:) = 0.
429	ziflprev(:) = 0.
430
431	icesed_flux(:) = 0.
432
433
434	IF (iflag_evap_prec .GE. 4) THEN
435
436	!Partitionning between precipitation coming from clouds and that coming from CS
437
438	!0) Calculate tot_zneb, total cloud fraction above the cloud
439	!assuming a maximum-random overlap (voir Jakob and Klein, 2000)
440
441	DO i=1, klon
442	tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i),zneb(i))) &
443	/(1.-min(zneb(i),1.-smallestreal))
444	d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
445	tot_zneb(i) = tot_znebn(i)
446
447
448	!1) Cloudy to clear air
449	d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i),znebprecipcld(i))
450	IF (znebprecipcld(i) .GT. 0.) THEN
451	d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
452	d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
453	ELSE
454	d_zrfl_cld_clr(i) = 0.
455	d_zifl_cld_clr(i) = 0.
456	ENDIF
457
458	!2) Clear to cloudy air
459	d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i) - d_tot_zneb(i) - zneb(i)))
460	IF (znebprecipclr(i) .GT. 0) THEN
461	d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
462	d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
463	ELSE
464	d_zrfl_clr_cld(i) = 0.
465	d_zifl_clr_cld(i) = 0.
466	ENDIF
467
468	!Update variables
469	znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
470	znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
471	zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
472	ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
473	zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
474	ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
475
476	! c_iso do the same thing for isotopes precip
477	ENDDO
478	ENDIF
479
480
481	! Autoconversion
482	! -------------------------------------------------------------------------------
483
484	! Initialisation
485	DO i = 1, klon
486	zrho(i) = pplay(i) / zt(i) / RD
487	zdz(i) = (paprsdn(i)-paprsup(i)) / (zrho(i)*RG)
488	iwc(i) = 0.
489	zneb(i) = MAX(rneb(i),seuil_neb)
490
491	! variables for calculation of quantity of condensates seen by the radiative scheme
492	! NB. only zradocond and zradoice are outputed, but zradoliq is used if ok_radocond_snow
493	! is TRUE
494	zradocond(i) = zoliq(i)/REAL(niter_lscp+1)
495	zradoliq(i) = zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
496	zradoice(i) = zoliq(i)*zfice(i)/REAL(niter_lscp+1)
497	ENDDO
498
499
500	DO n = 1, niter_lscp
501
502	DO i=1,klon
503	IF (rneb(i).GT.0.0) THEN
504	iwc(i) = zrho(i) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
505	ENDIF
506	ENDDO
507
508	CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, pt_pron_clds, zvelo)
509
510	DO i = 1, klon
511
512	IF (rneb(i).GT.0.0) THEN
513
514	! Initialization of zrain, zsnow and zprecip:
515	zrain=0.
516	zsnow=0.
517	zprecip=0.
518	! c_iso same init for isotopes. Externalisation?
519
520	IF (zneb(i).EQ.seuil_neb) THEN
521	zprecip = 0.0
522	zsnow = 0.0
523	zrain= 0.0
524	ELSE
525
526	IF (ptconv(i)) THEN ! if convective point
527	zcl=cld_lc_con
528	zct=1./cld_tau_con
529	zexpo=cld_expo_con
530	ffallv=ffallv_con
531	ELSE
532	zcl=cld_lc_lsc
533	zct=1./cld_tau_lsc
534	zexpo=cld_expo_lsc
535	ffallv=ffallv_lsc
536	ENDIF
537
538
539	! if vertical heterogeneity is taken into account, we use
540	! the "true" volume fraction instead of a modified
541	! surface fraction (which is larger and artificially
542	! reduces the in-cloud water).
543	IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
544	effective_zneb=ctot_vol(i)
545	ELSE
546	effective_zneb=zneb(i)
547	ENDIF
548
549
550	! Liquid water quantity to remove: zchau (Sundqvist, 1978)
551	! dqliq/dt=-qliq/tau(1-exp(-qcin/clw)*2)
552	!.........................................................
553	IF (iflag_autoconversion .EQ. 2) THEN
554	! two-steps resolution with niter_lscp=1 sufficient
555	! we first treat the second term (with exponential) in an explicit way
556	! and then treat the first term (-q/tau) in an exact way
557	zchau=zoliql(i)(1.-exp(-dtime/REAL(niter_lscp)zct &
558	(1.-exp(-(zoliql(i)/effective_zneb/zcl)*zexpo))))
559	ELSE
560	! old explicit resolution with subtimesteps
561	zchau = zctdtime/REAL(niter_lscp)zoliql(i) &
562	(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)*zexpo))
563	ENDIF
564
565
566	! Ice water quantity to remove (Zender & Kiehl, 1997)
567	! dqice/dt=1/rhod(rhowice*qice)/dz
568	!....................................
569	IF (iflag_autoconversion .EQ. 2) THEN
570	! exact resolution, niter_lscp=1 is sufficient but works only
571	! with iflag_vice=0
572	IF (zoliqi(i) .GT. 0.) THEN
573	zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
574	+dice_velodtime/REAL(niter_lscp)cice_velo/zdz(i)/zneb(i)zrho(i)ffallv)**(-1./dice_velo))
575	ELSE
576	zfroi=0.
577	ENDIF
578	ELSE
579	! old explicit resolution with subtimesteps
580	zfroi = dtime/REAL(niter_lscp)/zdz(i)zoliqi(i)zvelo(i)
581	ENDIF
582
583	zrain = MIN(MAX(zchau,0.0),zoliql(i))
584	zsnow = MIN(MAX(zfroi,0.0),zoliqi(i))
585	zprecip = MAX(zrain + zsnow,0.0)
586
587	ENDIF
588
589
590	IF (iflag_autoconversion .GE. 1) THEN
591	! debugged version with phase conservation through the autoconversion process
592	zoliql(i) = MAX(zoliql(i)-1.*zrain , 0.0)
593	zoliqi(i) = MAX(zoliqi(i)-1.*zsnow , 0.0)
594	zoliq(i) = MAX(zoliq(i)-zprecip , 0.0)
595	ELSE
596	! bugged version with phase resetting
597	zoliql(i) = MAX(zoliq(i)(1.-zfice(i))-1.zrain , 0.0)
598	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zsnow , 0.0)
599	zoliq(i) = MAX(zoliq(i)-zprecip , 0.0)
600	ENDIF
601
602	! c_iso: call isotope_conversion (for readibility) or duplicate
603
604	! variables for calculation of quantity of condensates seen by the radiative scheme
605	zradocond(i) = zradocond(i) + zoliq(i)/REAL(niter_lscp+1)
606	zradoliq(i) = zradoliq(i) + zoliql(i)/REAL(niter_lscp+1)
607	zradoice(i) = zradoice(i) + zoliqi(i)/REAL(niter_lscp+1)
608
609	!--Diagnostics
610	dqrauto(i) = dqrauto(i) + zrain / dtime
611	dqsauto(i) = dqsauto(i) + zsnow / dtime
612
613	ENDIF ! rneb >0
614
615	ENDDO ! i = 1,klon
616
617	ENDDO ! n = 1,niter
618
619	! Precipitation flux calculation
620
621	DO i = 1, klon
622
623	IF (iflag_evap_prec.GE.4) THEN
624	ziflprev(i)=ziflcld(i)
625	ELSE
626	ziflprev(i)=zifl(i)*zneb(i)
627	ENDIF
628
629	IF (rneb(i) .GT. 0.0) THEN
630
631	! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
632	! If T<0C, liquid precip are converted into ice, which leads to an increase in
633	! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
634	! taken into account through a linearization of the equations and by approximating
635	! the liquid precipitation process with a "threshold" process. We assume that
636	! condensates are not modified during this operation. Liquid precipitation is
637	! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
638	! Water vapor increases as well
639	! Note that compared to fisrtilp, we always assume iflag_bergeron=2
640
641	IF (ok_bug_phase_lscp) THEN
642	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
643	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
644	ELSE
645	zqpreci(i)=zcond(i)*zfice(i)-zoliqi(i)
646	zqprecl(i)=zcond(i)*(1.-zfice(i))-zoliql(i)
647	ENDIF
648	zcp=RCPD(1.0+RVTMP2(zq(i)+zmqc(i)+zcond(i)))
649	coef1 = rneb(i)RLSTT/zcpzdqsdT_raw(i)
650	! Computation of DT if all the liquid precip freezes
651	DeltaT = RLMLTzqprecl(i) / (zcp(1.+coef1))
652	! T should not exceed the freezing point
653	! that is Delta > RTT-zt(i)
654	DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
655	zt(i) = zt(i) + DeltaT
656	! water vaporization due to temp. increase
657	Deltaq = rneb(i)zdqsdT_raw(i)DeltaT
658	! we add this vaporized water to the total vapor and we remove it from the precip
659	zq(i) = zq(i) + Deltaq
660	! The three "max" lines herebelow protect from rounding errors
661	zcond(i) = max( zcond(i)- Deltaq, 0. )
662	! liquid precipitation converted to ice precip
663	Deltaqprecl = -zcp/RLMLT(1.+coef1)DeltaT
664	zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
665	! iced water budget
666	zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
667
668	! c_iso : mv here condensation of isotopes + redispatchage en precip
669
670	IF (iflag_evap_prec.GE.4) THEN
671	zrflcld(i) = zrflcld(i)+zqprecl(i) &
672	(paprsdn(i)-paprsup(i))/(RGdtime)
673	ziflcld(i) = ziflcld(i)+ zqpreci(i) &
674	(paprsdn(i)-paprsup(i))/(RGdtime)
675	znebprecipcld(i) = rneb(i)
676	zrfl(i) = zrflcld(i) + zrflclr(i)
677	zifl(i) = ziflcld(i) + ziflclr(i)
678	ELSE
679	zrfl(i) = zrfl(i)+ zqprecl(i) &
680	(paprsdn(i)-paprsup(i))/(RGdtime)
681	zifl(i) = zifl(i)+ zqpreci(i) &
682	(paprsdn(i)-paprsup(i))/(RGdtime)
683	ENDIF
684	! c_iso : same for isotopes
685
686	ENDIF ! rneb>0
687
688	ENDDO
689
690	!---------------------------------------------------------
691	!-- ICE SEDIMENTATION
692	!---------------------------------------------------------
693	IF ( ok_ice_sedim ) THEN
694
695	DO i = 1, klon
696	IF ( ( zoliqi(i) .GT. 0. ) .AND. ( rneb(i) .GT. eps ) .AND. pt_pron_clds(i) ) THEN
697
698	iwcg = MAX( zrho(i) * zoliqi(i) / zneb(i) * 1000., eps ) ! iwc in g/m3
699	!tempc = zt(i) - RTT ! celcius temp
700	!! so-called 'empirical parameterization' in Stubenrauch et al. 2019
701	!IF ( tempc .GE. -60. ) THEN
702	! icevelo = EXP( 1.9714 + 0.00216078 * tempc &
703	! - ( 0.0000414122 * tempc*2 + 0.00538922 tempc + 0.0516344 ) * LOG(iwcg) )
704	!ELSE
705	! icevelo = 65. * iwcg*0.2 ( 15000. / pplay(i) )**0.15
706	!ENDIF
707	!icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s
708
709	icevelo = fallice_sedim * 1.645 * ( iwcg / 1000. )**0.16
710
711	qice_sedim = zoliqi(i) * MIN(1., icevelo * dtime / zdz(i))
712	!--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two
713	!--times in the same timestep)
714	qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime)
715	!--Add tendencies
716	zoliqi(i) = zoliqi(i) - qice_sedim
717	zoliq(i) = zoliq(i) - qice_sedim
718	!--Convert to flux
719	icesed_flux(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime
720
721	!--Diagnostic tendencies
722	dqised(i) = dqised(i) - qice_sedim / dtime
723	ENDIF
724	ENDDO
725	ENDIF
726
727	! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
728	! if iflag_evap_prec>=4
729	IF (iflag_evap_prec.GE.4) THEN
730
731	DO i=1,klon
732
733	IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
734	znebprecipclr(i) = min(znebprecipclr(i),max( &
735	zrflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), &
736	ziflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
737	ELSE
738	znebprecipclr(i)=0.0
739	ENDIF
740
741	IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
742	znebprecipcld(i) = min(znebprecipcld(i), max( &
743	zrflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), &
744	ziflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
745	ELSE
746	znebprecipcld(i)=0.0
747	ENDIF
748	ENDDO
749
750	ENDIF
751
752	! we recalculate zradoice to account for contributions from the precipitation flux
753	! if ok_radocond_snow is true
754	IF ( ok_radocond_snow ) THEN
755	IF ( .NOT. iftop ) THEN
756	! for the solid phase (crystals + snowflakes)
757	! we recalculate zradoice by summing
758	! the ice content calculated in the mesh
759	! + the contribution of the non-evaporated snowfall
760	! from the overlying layer
761	DO i=1,klon
762	IF (ziflprev(i) .NE. 0.0) THEN
763	zradoice(i)=zoliqi(i)+zqpreci(i)+ziflprev(i)/zrho_up(i)/zvelo_up(i)
764	ELSE
765	zradoice(i)=zoliqi(i)+zqpreci(i)
766	ENDIF
767	zradocond(i)=zradoliq(i)+zradoice(i)
768	ENDDO
769	ENDIF
770	! keep in memory air density and ice fallspeed velocity
771	zrho_up(:) = zrho(:)
772	zvelo_up(:) = zvelo(:)
773	ENDIF
774
775	END SUBROUTINE histprecip_postcld
776
777
778	!----------------------------------------------------------------
779	! Computes the processes-oriented precipitation formulations for
780	! evaporation and sublimation
781	!
782	SUBROUTINE poprecip_precld( &
783	klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
784	qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, icesed_flux, &
785	cldfra, qvc, qliq, qice, &
786	rain, rainclr, raincld, snow, snowclr, snowcld, &
787	dqreva, dqssub)
788
789	USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
790	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
791	USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds
792	USE lmdz_lscp_ini, ONLY : ok_growth_precip_deposition, ok_ice_sedim
793	USE lmdz_lscp_ini, ONLY : eps, temp_nowater
794	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
795
796	IMPLICIT NONE
797
798
799	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
800	REAL, INTENT(IN) :: dtime !--time step [s]
801	LOGICAL, INTENT(IN) :: iftop !--if top of the column
802
803
804	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
805	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
806	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
807
808	REAL, INTENT(INOUT), DIMENSION(klon) :: temp !--current temperature [K]
809	REAL, INTENT(IN), DIMENSION(klon) :: tempupnew !--updated temperature of the overlying layer [K]
810
811	REAL, INTENT(INOUT), DIMENSION(klon) :: qvap !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
812	REAL, INTENT(INOUT), DIMENSION(klon) :: qprecip !--specific humidity in the precipitation falling from the upper layer [kg/kg]
813
814	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
815	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
816
817	REAL, INTENT(IN), DIMENSION(klon) :: qvapclrup !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
818	REAL, INTENT(IN), DIMENSION(klon) :: qtotupnew !--total specific humidity IN THE LAYER ABOVE [kg/kg]
819	REAL, INTENT(IN), DIMENSION(klon) :: icesed_flux !--sedimentated ice water flux [kg/s/m2]
820
821	REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
822	REAL, INTENT(INOUT), DIMENSION(klon) :: qvc !--cloud water vapor at the beginning of lscp (ratio wrt total water) - used only if the cloud properties are advected [kg/kg]
823	REAL, INTENT(INOUT), DIMENSION(klon) :: qliq !--liquid water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
824	REAL, INTENT(INOUT), DIMENSION(klon) :: qice !--ice water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
825
826
827	REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
828	REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
829	REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
830	REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
831	REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
832	REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
833
834	REAL, INTENT(OUT), DIMENSION(klon) :: dqreva !--rain tendency due to evaporation [kg/kg/s]
835	REAL, INTENT(OUT), DIMENSION(klon) :: dqssub !--snow tendency due to sublimation [kg/kg/s]
836
837
838	!--Integer for interating over klon
839	INTEGER :: i
840	!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
841	REAL, DIMENSION(klon) :: dhum_to_dflux
842	!--Air density [kg/m3] and layer thickness [m]
843	REAL, DIMENSION(klon) :: rho, dz
844	!--Temporary precip fractions and fluxes for the evaporation
845	REAL, DIMENSION(klon) :: precipfracclr_tmp, precipfraccld_tmp
846	REAL, DIMENSION(klon) :: rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp
847
848	!--Saturation values
849	REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
850	!--Vapor in the clear sky and cloud
851	REAL :: qvapclr, qvapcld
852	!--Fluxes tendencies because of evaporation and sublimation
853	REAL :: dprecip_evasub_max, dprecip_evasub_tot
854	REAL :: drainclreva, draincldeva, dsnowclrsub, dsnowcldsub
855	!--Specific humidity tendencies because of evaporation and sublimation
856	REAL :: dqrevap, dqssubl
857	!--Specific heat constant
858	REAL :: cpair, cpw
859	!--Specific ice water content sedimentated
860	REAL :: qice_sedim
861
862	!--Initialisation
863	qzero(:) = 0.
864	dqreva(:) = 0.
865	dqssub(:) = 0.
866
867	!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
868	dhum_to_dflux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime
869	rho(:) = pplay(:) / temp(:) / RD
870	dz(:) = ( paprsdn(:) - paprsup(:) ) / RG / rho(:)
871
872	!--Calculation of saturation specific humidity
873	!--depending on temperature:
874	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:))
875	!--wrt liquid water
876	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
877	!--wrt ice
878	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
879
880
881
882	!--First step consists in "thermalizing" the layer:
883	!--as the flux of precip from layer above "advects" some heat (as the precip is at the temperature
884	!--of the overlying layer) we recalculate a mean temperature that both the air and the precip in the
885	!--layer have.
886
887	IF (iftop) THEN
888
889	DO i = 1, klon
890	qprecip(i) = 0.
891	ENDDO
892
893	ELSE
894
895	DO i = 1, klon
896	!--No condensed water so cp=cp(vapor+dry air)
897	!-- RVTMP2=rcpv/rcpd-1
898	cpair = RCPD * ( 1. + RVTMP2 * qvap(i) )
899	cpw = RCPD * RVTMP2
900	!--qprecip has to be thermalized with
901	!--layer's air so that precipitation at the ground has the
902	!--same temperature as the lowermost layer
903	!--we convert the flux into a specific quantity qprecip
904	qprecip(i) = ( rain(i) + snow(i) ) / dhum_to_dflux(i)
905	!-- t(i,k+1) + d_t(i,k+1): new temperature of the overlying layer
906	temp(i) = ( tempupnew(i) * qprecip(i) * cpw + cpair * temp(i) ) &
907	/ ( cpair + qprecip(i) * cpw )
908	ENDDO
909
910	ENDIF
911
912
913	!--Initialise the precipitation fractions and fluxes
914	DO i = 1, klon
915	precipfracclr_tmp(i) = precipfracclr(i)
916	precipfraccld_tmp(i) = precipfraccld(i)
917	rainclr_tmp(i) = rainclr(i)
918	snowclr_tmp(i) = snowclr(i)
919	raincld_tmp(i) = raincld(i)
920	snowcld_tmp(i) = snowcld(i)
921	ENDDO
922
923	!--If we relax the LTP assumption that the cloud is the same than in the
924	!--gridbox above, we can change the tmp quantities
925	IF ( ok_ice_supersat ) THEN
926	!--Update the precipitation fraction using advected cloud fraction
927	CALL poprecip_fracupdate( &
928	klon, cldfra, precipfracclr_tmp, precipfraccld_tmp, &
929	rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp)
930
931	! here we could put an ELSE statement and do one iteration of the condensation scheme
932	! (but it would only diagnose cloud from lognormal - link with thermals?)
933
934	DO i = 1, klon
935	!--We cannot have a total fraction greater than the one in the gridbox above
936	!--because new cloud formation has not yet occured.
937	!--If this happens, we reduce the precip frac in the cloud, because it can
938	!--only mean that the cloud fraction at this level is greater than the total
939	!--precipitation fraction of the level above, and the precip frac in clear sky
940	!--is zero.
941	!--NB. this only works because we assume a max-random cloud and precip overlap
942	precipfraccld_tmp(i) = MIN( precipfraccld_tmp(i), precipfracclr(i) + precipfraccld(i) )
943	ENDDO
944
945	ENDIF
946	!--At this stage, we guarantee that
947	!-- precipfracclr + precipfraccld = precipfracclr_tmp + precipfraccld_tmp
948
949
950	DO i = 1, klon
951
952	dqrevap = 0.
953	dqssubl = 0.
954	!--If there is precipitation from the layer above
955	IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
956
957	!--Init
958	drainclreva = 0.
959	draincldeva = 0.
960	dsnowclrsub = 0.
961	dsnowcldsub = 0.
962
963	!--Init the properties of the air where reevaporation occurs
964	IF ( ok_ice_supersat ) THEN
965	!--We can diagnose the water vapor in the cloud using the advected
966	!--water vapor in the cloud and cloud fraction
967	IF ( ok_unadjusted_clouds .AND. ( cldfra(i) .GT. eps ) ) THEN
968	qvapcld = qvc(i) / cldfra(i)
969	ENDIF
970	!--We can diagnose completely the water vapor in clear sky, because all
971	!--the needed variables (ice, liq, vapor in cld, cloud fraction) are
972	!--advected
973	!--Note that qvap(i) is the total water in the gridbox
974	IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
975	qvapclr = ( qvap(i) - qice(i) - qliq(i) - qvc(i) ) / ( 1. - cldfra(i) )
976	ENDIF
977	ELSEIF ( ok_corr_vap_evasub ) THEN
978	!--Corrected version from Ludo - we use the same water ratio between
979	!--the clear and the cloudy sky as in the layer above. This
980	!--extends the assumption that the cloud fraction is the same
981	!--as the layer above. This is assumed only for the evap / subl
982	!--process
983	!--Note that qvap(i) is the total water in the gridbox, and
984	!--precipfraccld(i) is the cloud fraction in the layer above
985	IF ( ( 1. - precipfraccld(i) ) .GT. eps ) THEN
986	qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
987	ELSE
988	qvapclr = 0.
989	ENDIF
990	IF ( precipfraccld(i) .GT. eps ) THEN
991	qvapcld = MAX(qtotupnew(i)-qvapclrup(i) , 0.) / qtotupnew(i) * qvap(i) / precipfraccld(i)
992	ELSE
993	qvapcld = 0.
994	ENDIF
995	ELSE
996	!--Legacy version from Ludo - we use the total specific humidity
997	!--for the evap / subl process
998	qvapclr = qvap(i)
999	qvapcld = qvap(i)
1000	ENDIF
1001
1002
1003	!--------------------
1004	!-- ICE SEDIMENTATION
1005	!--------------------
1006
1007	!--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
1008	!--added to the total water content of the gridbox
1009	IF ( icesed_flux(i) .GT. 0. ) THEN
1010	qice_sedim = icesed_flux(i) / dhum_to_dflux(i)
1011	!--Vapor is updated after evaporation/sublimation (it is increased)
1012	qvap(i) = qvap(i) + qice_sedim
1013	!--Air and precip temperature (i.e., gridbox temperature)
1014	!--is updated due to latent heat cooling
1015	temp(i) = temp(i) &
1016	- qice_sedim * RLSTT / RCPD &
1017	/ ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
1018	ENDIF ! ok_ice_sedim
1019
1020
1021	!---------------------------
1022	!-- EVAP / SUBL IN CLEAR SKY
1023	!---------------------------
1024
1025	IF ( precipfracclr_tmp(i) .GT. eps ) THEN
1026
1027	!--Evaporation of liquid precipitation coming from above
1028	!--in the clear sky only
1029	!--dprecip/dz = -beta(1-qvap/qsat)(precip**expo_eva)
1030	!--formula from Sundqvist 1988, Klemp & Wilhemson 1978
1031	!--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly)
1032	!--which does not need a barrier on rainclr, because included in the formula
1033	drainclreva = precipfracclr_tmp(i) * MAX(0., &
1034	- coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) &
1035	+ ( rainclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_eva ) &
1036	)**( 1. / ( 1. - expo_eva ) ) - rainclr_tmp(i)
1037
1038	!--Evaporation is limited by 0
1039	!--NB. with ok_ice_supersat activated, this barrier should be useless
1040	drainclreva = MIN(0., drainclreva)
1041
1042
1043	!--Sublimation of the solid precipitation coming from above
1044	!--(same formula as for liquid precip)
1045	!--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly)
1046	!--which does not need a barrier on snowclr, because included in the formula
1047	dsnowclrsub = precipfracclr_tmp(i) * MAX(0., &
1048	- coef_sub * ( 1. - expo_sub ) * (1. - qvapclr / qsati(i)) * dz(i) &
1049	+ ( snowclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_sub ) &
1050	)**( 1. / ( 1. - expo_sub ) ) - snowclr_tmp(i)
1051
1052	!--If ice supersaturation is activated, we allow vapor to condense on falling snow
1053	!--i.e., having a positive dsnowclrsub
1054	IF ( .NOT. ok_ice_supersat ) THEN
1055	!--Sublimation is limited by 0
1056	dsnowclrsub = MIN(0., dsnowclrsub)
1057	ENDIF
1058
1059
1060	!--Evaporation limit: we ensure that the layer's fraction below
1061	!--the clear sky does not reach saturation. In this case, we
1062	!--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
1063	!--Max evaporation is computed not to saturate the clear sky precip fraction
1064	!--(i.e., the fraction where evaporation occurs)
1065	!--It is expressed as a max flux dprecip_evasub_max
1066
1067	IF ( .NOT. ok_ice_supersat ) THEN
1068	dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr_tmp(i)) &
1069	* dhum_to_dflux(i)
1070	ELSE
1071	dprecip_evasub_max = ( qvapclr - qsat(i) ) * precipfracclr_tmp(i) &
1072	* dhum_to_dflux(i)
1073	ENDIF
1074	dprecip_evasub_tot = drainclreva + dsnowclrsub
1075
1076	!--Barriers
1077	!--If activates if the total is LOWER than the max because
1078	!--everything is negative
1079	IF ( (( dprecip_evasub_max .LT. 0. ) .AND. &
1080	( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. &
1081	(( dprecip_evasub_max .GT. 0. ) .AND. &
1082	( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN
1083	drainclreva = dprecip_evasub_max * drainclreva / dprecip_evasub_tot
1084	dsnowclrsub = dprecip_evasub_max * dsnowclrsub / dprecip_evasub_tot
1085	ENDIF
1086
1087	ELSE
1088	!--All the precipitation is sublimated if the fraction is zero
1089	drainclreva = - rainclr_tmp(i)
1090	dsnowclrsub = - snowclr_tmp(i)
1091
1092	ENDIF ! precipfracclr_tmp .GT. eps
1093
1094
1095	!---------------------------
1096	!-- EVAP / SUBL IN THE CLOUD
1097	!---------------------------
1098
1099	IF ( ok_growth_precip_deposition .AND. ( temp(i) .LE. RTT ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
1100	!--Evaporation of liquid precipitation coming from above
1101	!--in the cloud only
1102	!--dprecip/dz = -beta(1-qvap/qsat)(precip**expo_eva)
1103	!--formula from Sundqvist 1988, Klemp & Wilhemson 1978
1104	!--Exact explicit formulation (raincld is resolved exactly, qvap explicitly)
1105	!--which does not need a barrier on raincld, because included in the formula
1106	draincldeva = precipfraccld_tmp(i) * MAX(0., &
1107	- coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
1108	+ ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
1109	)**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
1110
1111	!--Evaporation is limited by 0
1112	!--NB. with ok_ice_supersat activated, this barrier should be useless
1113	draincldeva = MIN(0., draincldeva)
1114
1115
1116	!--Sublimation of the solid precipitation coming from above
1117	!--(same formula as for liquid precip)
1118	!--Exact explicit formulation (snowcld is resolved exactly, qvap explicitly)
1119	!--which does not need a barrier on snowcld, because included in the formula
1120	dsnowcldsub = precipfraccld_tmp(i) * MAX(0., &
1121	- coef_sub * ( 1. - expo_sub ) * (1. - qvapcld / qsati(i)) * dz(i) &
1122	+ ( snowcld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_sub ) &
1123	)**( 1. / ( 1. - expo_sub ) ) - snowcld_tmp(i)
1124
1125	!--There is no barrier because we want deposition to occur if it is possible
1126
1127
1128	!--Evaporation limit: we ensure that the layer's fraction below
1129	!--the clear sky does not reach saturation. In this case, we
1130	!--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
1131	!--Max evaporation is computed not to saturate the clear sky precip fraction
1132	!--(i.e., the fraction where evaporation occurs)
1133	!--It is expressed as a max flux dprecip_evasub_max
1134
1135	dprecip_evasub_max = ( qvapcld - qsat(i) ) * precipfraccld_tmp(i) * dhum_to_dflux(i)
1136	dprecip_evasub_tot = draincldeva + dsnowcldsub
1137
1138	!--Barriers
1139	!--If activates if the total is LOWER than the max because
1140	!--everything is negative
1141	IF ( (( dprecip_evasub_max .LT. 0. ) .AND. &
1142	( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. &
1143	(( dprecip_evasub_max .GT. 0. ) .AND. &
1144	( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN
1145	draincldeva = dprecip_evasub_max * draincldeva / dprecip_evasub_tot
1146	dsnowcldsub = dprecip_evasub_max * dsnowcldsub / dprecip_evasub_tot
1147	ENDIF
1148
1149	ENDIF ! ok_unadjusted_clouds
1150
1151
1152	!--New solid and liquid precipitation fluxes after evap and sublimation
1153	dqrevap = ( drainclreva + draincldeva ) / dhum_to_dflux(i)
1154	dqssubl = ( dsnowclrsub + dsnowcldsub ) / dhum_to_dflux(i)
1155
1156	!--Vapor is updated after evaporation/sublimation (it is increased)
1157	qvap(i) = qvap(i) - dqrevap - dqssubl
1158	!--qprecip is the total condensed water in the precip flux (it is decreased)
1159	qprecip(i) = qprecip(i) + dqrevap + dqssubl
1160	!--Air and precip temperature (i.e., gridbox temperature)
1161	!--is updated due to latent heat cooling
1162	temp(i) = temp(i) &
1163	+ dqrevap * RLVTT / RCPD &
1164	/ ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) &
1165	+ dqssubl * RLSTT / RCPD &
1166	/ ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
1167
1168	!--Add tendencies
1169	!--The MAX is needed because in some cases, the flux can be slightly
1170	!--negative (numerical precision)
1171	rainclr_tmp(i) = MAX(0., rainclr_tmp(i) + drainclreva)
1172	snowclr_tmp(i) = MAX(0., snowclr_tmp(i) + dsnowclrsub)
1173	raincld_tmp(i) = MAX(0., raincld_tmp(i) + draincldeva)
1174	snowcld_tmp(i) = MAX(0., snowcld_tmp(i) + dsnowcldsub)
1175
1176
1177	!--We reattribute fluxes according to initial fractions, using proportionnality
1178	!--Note that trough this process, we conserve total precip fluxes
1179	!-- rainclr_tmp + raincld_tmp = rainclr + raincld (same for snow)
1180	!--If the cloud has increased, a part of the precip in clear sky falls in clear sky,
1181	!--the rest falls in the cloud
1182	IF ( ( precipfraccld(i) .GT. eps ) .AND. ( precipfraccld_tmp(i) .GT. precipfraccld(i) ) ) THEN
1183	!--All the precip from cloud falls in the cloud
1184	raincld(i) = MAX(0., raincld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i))
1185	rainclr(i) = MAX(0., rainclr_tmp(i) + raincld_tmp(i) - raincld(i))
1186	snowcld(i) = MAX(0., snowcld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i))
1187	snowclr(i) = MAX(0., snowclr_tmp(i) + snowcld_tmp(i) - snowcld(i))
1188
1189	!--If the cloud has narrowed, a part of the precip in cloud falls in clear sky,
1190	!--the rest falls in the cloud
1191	ELSEIF ( ( precipfracclr(i) .GT. eps ) .AND. ( precipfracclr_tmp(i) .GT. precipfracclr(i) ) ) THEN
1192	!--All the precip from clear sky falls in clear sky
1193	rainclr(i) = MAX(0., rainclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i))
1194	raincld(i) = MAX(0., raincld_tmp(i) + rainclr_tmp(i) - rainclr(i))
1195	snowclr(i) = MAX(0., snowclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i))
1196	snowcld(i) = MAX(0., snowcld_tmp(i) + snowclr_tmp(i) - snowclr(i))
1197
1198	ELSE
1199	!--If the cloud stays the same or if there is no cloud above and
1200	!--in the current layer, there is no reattribution
1201	rainclr(i) = rainclr_tmp(i)
1202	raincld(i) = raincld_tmp(i)
1203	snowclr(i) = snowclr_tmp(i)
1204	snowcld(i) = snowcld_tmp(i)
1205	ENDIF
1206
1207	!--If there is no more precip fluxes, the precipitation fraction is set to 0
1208	IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0.
1209	IF ( ( raincld(i) + snowcld(i) ) .LE. 0. ) precipfraccld(i) = 0.
1210
1211	!--Calculation of the total fluxes
1212	rain(i) = rainclr(i) + raincld(i)
1213	snow(i) = snowclr(i) + snowcld(i)
1214
1215	ELSEIF ( ( rain(i) + snow(i) ) .LE. 0. ) THEN
1216	!--If no precip, we reinitialize the cloud fraction used for the precip to 0
1217	precipfraccld(i) = 0.
1218	precipfracclr(i) = 0.
1219
1220	ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. )
1221
1222	!--Diagnostic tendencies
1223	dqssub(i) = dqssubl / dtime
1224	dqreva(i) = dqrevap / dtime
1225
1226	ENDDO ! loop on klon
1227
1228	END SUBROUTINE poprecip_precld
1229
1230
1231	!----------------------------------------------------------------
1232	! Computes the precipitation fraction update
1233	! The goal of this routine is to reattribute precipitation fractions
1234	! and fluxes to clear or cloudy air, depending on the variation of
1235	! the cloud fraction on the vertical dimension. We assume a
1236	! maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
1237	! and LTP thesis, 2021)
1238	! NB. in fact, we assume a maximum-random overlap of the total precip. frac
1239	!
1240	SUBROUTINE poprecip_fracupdate( &
1241	klon, cldfra, precipfracclr, precipfraccld, &
1242	rainclr, raincld, snowclr, snowcld)
1243
1244	USE lmdz_lscp_ini, ONLY : eps
1245
1246	IMPLICIT NONE
1247
1248	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
1249
1250	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction [-]
1251	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
1252	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
1253	!--NB. at the end of the routine, becomes the fraction of precip
1254	!--in the current layer
1255
1256	REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
1257	REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
1258	REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
1259	REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
1260
1261	!--Local variables
1262	INTEGER :: i
1263	REAL :: dcldfra
1264	REAL :: precipfractot
1265	REAL :: dprecipfracclr, dprecipfraccld
1266	REAL :: drainclr, dsnowclr
1267	REAL :: draincld, dsnowcld
1268
1269
1270	DO i = 1, klon
1271
1272	!--Initialisation
1273	precipfractot = precipfracclr(i) + precipfraccld(i)
1274
1275	!--Instead of using the cloud cover which was use in LTP thesis, we use the
1276	!--total precip. fraction to compute the maximum-random overlap. This is
1277	!--because all the information of the cloud cover is embedded into
1278	!--precipfractot, and this allows for taking into account the potential
1279	!--reduction of the precipitation fraction because either the flux is too
1280	!--small (see barrier at the end of poprecip_postcld) or the flux is completely
1281	!--evaporated (see barrier at the end of poprecip_precld)
1282	!--NB. precipfraccld(i) is here the cloud fraction of the layer above
1283	!precipfractot = 1. - ( 1. - precipfractot ) * &
1284	! ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
1285	! / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
1286
1287
1288	IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN
1289	precipfractot = 1.
1290	ELSE
1291	precipfractot = 1. - ( 1. - precipfractot ) * &
1292	( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
1293	/ ( 1. - precipfraccld(i) )
1294	ENDIF
1295
1296	!--precipfraccld(i) is here the cloud fraction of the layer above
1297	dcldfra = cldfra(i) - precipfraccld(i)
1298	!--Tendency of the clear-sky precipitation fraction. We add a MAX on the
1299	!--calculation of the current CS precip. frac.
1300	!dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
1301	!--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated
1302	!--if precipfractot < cldfra)
1303	dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i)
1304	!--Tendency of the cloudy precipitation fraction. We add a MAX on the
1305	!--calculation of the current CS precip. frac.
1306	!dprecipfraccld = MAX( dcldfra , - precipfraccld(i) )
1307	!--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated
1308	!--if cldfra < 0)
1309	dprecipfraccld = dcldfra
1310
1311
1312	!--If the cloud extends
1313	IF ( dprecipfraccld .GT. 0. ) THEN
1314	!--If there is no CS precip, nothing happens.
1315	!--If there is, we reattribute some of the CS precip flux
1316	!--to the cloud precip flux, proportionnally to the
1317	!--decrease of the CS precip fraction
1318	IF ( precipfracclr(i) .LT. eps ) THEN
1319	drainclr = 0.
1320	dsnowclr = 0.
1321	ELSE
1322	drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i)
1323	dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i)
1324	ENDIF
1325	!--If the cloud narrows
1326	ELSEIF ( dprecipfraccld .LT. 0. ) THEN
1327	!--We reattribute some of the cloudy precip flux
1328	!--to the CS precip flux, proportionnally to the
1329	!--decrease of the cloud precip fraction
1330	IF ( precipfraccld(i) .LT. eps ) THEN
1331	draincld = 0.
1332	dsnowcld = 0.
1333	ELSE
1334	draincld = dprecipfraccld / precipfraccld(i) * raincld(i)
1335	dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i)
1336	ENDIF
1337	drainclr = - draincld
1338	dsnowclr = - dsnowcld
1339	!--If the cloud stays the same or if there is no cloud above and
1340	!--in the current layer, nothing happens
1341	ELSE
1342	drainclr = 0.
1343	dsnowclr = 0.
1344	ENDIF
1345
1346	!--We add the tendencies
1347	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1348	precipfraccld(i) = precipfraccld(i) + dprecipfraccld
1349	precipfracclr(i) = precipfracclr(i) + dprecipfracclr
1350	rainclr(i) = MAX(0., rainclr(i) + drainclr)
1351	snowclr(i) = MAX(0., snowclr(i) + dsnowclr)
1352	raincld(i) = MAX(0., raincld(i) - drainclr)
1353	snowcld(i) = MAX(0., snowcld(i) - dsnowclr)
1354
1355	ENDDO
1356
1357	END SUBROUTINE poprecip_fracupdate
1358
1359
1360	!----------------------------------------------------------------
1361	! Computes the processes-oriented precipitation formulations for
1362	! - autoconversion (auto) via a deposition process
1363	! - aggregation (agg)
1364	! - riming (rim)
1365	! - collection (col)
1366	! - melting (melt)
1367	! - freezing (freez)
1368	!
1369	SUBROUTINE poprecip_postcld( &
1370	klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
1371	temp, qvap, qliq, qice, icefrac, cldfra, icesed_flux, &
1372	precipfracclr, precipfraccld, &
1373	rain, rainclr, raincld, snow, snowclr, snowcld, &
1374	qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
1375	dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez, dqised)
1376
1377	USE lmdz_lscp_ini, ONLY : prt_level, lunout
1378	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
1379	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
1380
1381	USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb, &
1382	cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, &
1383	thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
1384	rho_rain, r_rain, r_snow, rho_ice, &
1385	tau_auto_snow_min, tau_auto_snow_max, &
1386	expo_tau_auto_snow, thresh_precip_frac, eps, &
1387	gamma_melt, alpha_freez, beta_freez, temp_nowater, &
1388	iflag_cloudth_vert, iflag_rain_incloud_vol, &
1389	cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez, &
1390	rain_fallspeed_clr, rain_fallspeed_cld, &
1391	snow_fallspeed_clr, snow_fallspeed_cld, &
1392	ok_ice_sedim, fallice_sedim
1393
1394
1395	IMPLICIT NONE
1396
1397	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
1398	REAL, INTENT(IN) :: dtime !--time step [s]
1399
1400	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
1401	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
1402	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
1403
1404	REAL, INTENT(IN), DIMENSION(klon) :: ctot_vol !--volumic cloud fraction [-]
1405	LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv !--true if we are in a convective point
1406
1407	REAL, INTENT(INOUT), DIMENSION(klon) :: temp !--current temperature [K]
1408	REAL, INTENT(INOUT), DIMENSION(klon) :: qvap !--current water vapor specific humidity [kg/kg]
1409	REAL, INTENT(INOUT), DIMENSION(klon) :: qliq !--current liquid water specific humidity [kg/kg]
1410	REAL, INTENT(INOUT), DIMENSION(klon) :: qice !--current ice water specific humidity [kg/kg]
1411	REAL, INTENT(IN), DIMENSION(klon) :: icefrac !--ice fraction [-]
1412	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction [-]
1413
1414	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
1415	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
1416	!--NB. at the end of the routine, becomes the fraction of precip
1417	!--in the current layer
1418
1419	REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
1420	REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
1421	REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
1422	REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
1423	REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
1424	REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
1425
1426	REAL, INTENT(OUT), DIMENSION(klon) :: icesed_flux !--flux of sedimentated ice [kg/s/m2]
1427	REAL, INTENT(OUT), DIMENSION(klon) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg]
1428	REAL, INTENT(OUT), DIMENSION(klon) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg]
1429	REAL, INTENT(OUT), DIMENSION(klon) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
1430	REAL, INTENT(OUT), DIMENSION(klon) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
1431	REAL, INTENT(OUT), DIMENSION(klon) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
1432	REAL, INTENT(OUT), DIMENSION(klon) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
1433	REAL, INTENT(OUT), DIMENSION(klon) :: dqsrim !--snow tendency due to riming [kg/kg/s]
1434	REAL, INTENT(OUT), DIMENSION(klon) :: dqsmelt !--snow tendency due to melting [kg/kg/s]
1435	REAL, INTENT(OUT), DIMENSION(klon) :: dqrmelt !--rain tendency due to melting [kg/kg/s]
1436	REAL, INTENT(OUT), DIMENSION(klon) :: dqsfreez !--snow tendency due to freezing [kg/kg/s]
1437	REAL, INTENT(OUT), DIMENSION(klon) :: dqrfreez !--rain tendency due to freezing [kg/kg/s]
1438	REAL, INTENT(INOUT), DIMENSION(klon) :: dqised !--ice water content tendency due to sedimentation of ice crystals [kg/kg/s]
1439
1440
1441
1442	!--Local variables
1443
1444	INTEGER :: i
1445	!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
1446	REAL, DIMENSION(klon) :: dhum_to_dflux
1447	REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip
1448
1449	!--Collection, aggregation and riming
1450	REAL :: eff_cldfra
1451	REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
1452	REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
1453	REAL :: rho_snow
1454	REAL :: dqlcol !--loss of liquid cloud content due to collection by rain [kg/kg/s]
1455	REAL :: dqiagg !--loss of ice cloud content due to collection by aggregation [kg/kg/s]
1456	REAL :: dqlrim !--loss of liquid cloud content due to riming on snow [kg/kg/s]
1457
1458	!--Autoconversion
1459	REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
1460	REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
1461	REAL :: dqlauto !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
1462	REAL :: dqiauto !--loss of ice cloud content due to autoconversion to snow [kg/kg/s]
1463
1464	!--Melting
1465	REAL :: dqsmelt_max, air_thermal_conduct
1466	REAL :: nb_snowflake_clr, nb_snowflake_cld
1467	REAL :: capa_snowflake, temp_wetbulb
1468	REAL :: rho, r_ice
1469	REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
1470	REAL, DIMENSION(klon) :: qzero, qsat, dqsat
1471
1472	!--Freezing
1473	REAL :: dqrfreez_max
1474	REAL :: tau_freez
1475	REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2
1476	REAL :: coef_freez
1477	REAL :: dqifreez !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
1478	REAL :: Eff_rain_ice
1479
1480	!--Ice sedimentation
1481	REAL :: iwcg, tempc, icevelo
1482	REAL :: dz, qice_sedim
1483
1484
1485	!--Initialisation of variables
1486
1487
1488	qzero(:) = 0.
1489	icesed_flux(:) = 0.
1490
1491	dqrcol(:) = 0.
1492	dqsagg(:) = 0.
1493	dqsauto(:) = 0.
1494	dqrauto(:) = 0.
1495	dqsrim(:) = 0.
1496	dqrmelt(:) = 0.
1497	dqsmelt(:) = 0.
1498	dqrfreez(:) = 0.
1499	dqsfreez(:) = 0.
1500
1501
1502	!--Update the precipitation fraction following cloud formation
1503	CALL poprecip_fracupdate( &
1504	klon, cldfra, precipfracclr, precipfraccld, &
1505	rainclr, raincld, snowclr, snowcld)
1506
1507
1508	DO i = 1, klon
1509
1510	!--Variables initialisation
1511	dqlcol = 0.
1512	dqiagg = 0.
1513	dqiauto = 0.
1514	dqlauto = 0.
1515	dqlrim = 0.
1516
1517	!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
1518	dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
1519	qtot(i) = qvap(i) + qliq(i) + qice(i) &
1520	+ ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i)
1521
1522	!--If vertical heterogeneity is taken into account, we use
1523	!--the "true" volume fraction instead of a modified
1524	!--surface fraction (which is larger and artificially
1525	!--reduces the in-cloud water).
1526	IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
1527	eff_cldfra = ctot_vol(i)
1528	ELSE
1529	eff_cldfra = cldfra(i)
1530	ENDIF
1531
1532
1533	!--Start precipitation growth processes
1534
1535	!--If the cloud is big enough, the precipitation processes activate
1536	! TODO met on seuil_neb ici ?
1537	IF ( cldfra(i) .GE. seuil_neb ) THEN
1538
1539	!---------------------------------------------------------
1540	!-- COLLECTION AND AGGREGATION
1541	!---------------------------------------------------------
1542	!--Collection: processus through which rain collects small liquid droplets
1543	!--in suspension, and add it to the rain flux
1544	!--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
1545	!--Those processes are treated before autoconversion because we do not
1546	!--want to collect/aggregate the newly formed fluxes, which already
1547	!--"saw" the cloud as they come from it
1548	!--The formulas come from Muench and Lohmann 2020
1549
1550	!--gamma_col: tuning coefficient [-]
1551	!--rho_rain: volumic mass of rain [kg/m3]
1552	!--r_rain: size of the rain droplets [m]
1553	!--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
1554	!--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
1555	!--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
1556	!--get in-cloud mean quantities. The two divisions by eff_cldfra are
1557	!--then simplified.
1558
1559	!--The collection efficiency is perfect.
1560	Eff_rain_liq = 1.
1561	coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
1562	IF ( raincld(i) .GT. 0. ) THEN
1563	!--Exact explicit version, which does not need a barrier because of
1564	!--the exponential decrease
1565	dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
1566
1567	!--Add tendencies
1568	qliq(i) = qliq(i) + dqlcol
1569	raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i)
1570
1571	!--Diagnostic tendencies
1572	dqrcol(i) = - dqlcol / dtime
1573	ENDIF
1574
1575	!--Same as for aggregation
1576	!--Eff_snow_liq formula:
1577	!--it s a product of a collection efficiency and a sticking efficiency
1578	! Milbrandt and Yau formula that gives very low values:
1579	! Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
1580	! Lin 1983's formula
1581	Eff_snow_ice = EXP( 0.025 * MIN( ( temp(i) - RTT ), 0.) )
1582	!--rho_snow formula follows Brandes et al. 2007 (JAMC)
1583	rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
1584	coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
1585	IF ( snowcld(i) .GT. 0. ) THEN
1586	!--Exact explicit version, which does not need a barrier because of
1587	!--the exponential decrease
1588	dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
1589
1590	!--Add tendencies
1591	qice(i) = qice(i) + dqiagg
1592	snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i)
1593
1594	!--Diagnostic tendencies
1595	dqsagg(i) = - dqiagg / dtime
1596	ENDIF
1597
1598
1599	!---------------------------------------------------------
1600	!-- AUTOCONVERSION
1601	!---------------------------------------------------------
1602	!--Autoconversion converts liquid droplets/ice crystals into
1603	!--rain drops/snowflakes. It relies on the formulations by
1604	!--Sundqvist 1978.
1605
1606	!--If we are in a convective point, we have different parameters
1607	!--for the autoconversion
1608	IF ( ptconv(i) ) THEN
1609	qthresh_auto_rain = cld_lc_con
1610	qthresh_auto_snow = cld_lc_con_snow
1611
1612	tau_auto_rain = cld_tau_con
1613	!--tau for snow depends on the ice fraction in mixed-phase clouds
1614	tau_auto_snow = tau_auto_snow_max &
1615	+ ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow)
1616
1617	expo_auto_rain = cld_expo_con
1618	expo_auto_snow = cld_expo_con
1619	ELSE
1620	qthresh_auto_rain = cld_lc_lsc
1621	qthresh_auto_snow = cld_lc_lsc_snow
1622
1623	tau_auto_rain = cld_tau_lsc
1624	!--tau for snow depends on the ice fraction in mixed-phase clouds
1625	tau_auto_snow = tau_auto_snow_max &
1626	+ ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow)
1627
1628	expo_auto_rain = cld_expo_lsc
1629	expo_auto_snow = cld_expo_lsc
1630	ENDIF
1631
1632
1633	! Liquid water quantity to remove according to (Sundqvist, 1978)
1634	! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) )
1635	!
1636	!--And same formula for ice
1637	!
1638	!--We first treat the second term (with exponential) in an explicit way
1639	!--and then treat the first term (-q/tau) in an exact way
1640
1641	dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
1642	- ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
1643
1644	dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
1645	- ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
1646
1647
1648	!--Barriers so that we don't create more rain/snow
1649	!--than there is liquid/ice
1650	dqlauto = MAX( - qliq(i), dqlauto )
1651	dqiauto = MAX( - qice(i), dqiauto )
1652
1653	!--Add tendencies
1654	qliq(i) = qliq(i) + dqlauto
1655	qice(i) = qice(i) + dqiauto
1656	raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
1657	snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
1658
1659	!--Diagnostic tendencies
1660	dqsauto(i) = - dqiauto / dtime
1661	dqrauto(i) = - dqlauto / dtime
1662
1663
1664	!---------------------------------------------------------
1665	!-- RIMING
1666	!---------------------------------------------------------
1667	!--Process which converts liquid droplets in suspension into
1668	!--snow because of the collision between
1669	!--those and falling snowflakes.
1670	!--The formula comes from Muench and Lohmann 2020
1671	!--NB.: this process needs a temperature adjustment
1672
1673	!--Eff_snow_liq formula: following Ferrier 1994,
1674	!--assuming 1
1675	Eff_snow_liq = 1.0
1676	!--rho_snow formula follows Brandes et al. 2007 (JAMC)
1677	rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
1678	coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
1679	IF ( snowcld(i) .GT. 0. ) THEN
1680	!--Exact version, which does not need a barrier because of
1681	!--the exponential decrease
1682	dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. )
1683
1684	!--Add tendencies
1685	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1686	qliq(i) = qliq(i) + dqlrim
1687	snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i)
1688
1689	!--Temperature adjustment with the release of latent
1690	!--heat because of solid condensation
1691	temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
1692	/ ( 1. + RVTMP2 * qtot(i) )
1693
1694	!--Diagnostic tendencies
1695	dqsrim(i) = - dqlrim / dtime
1696	ENDIF
1697
1698	ENDIF ! cldfra .GE. seuil_neb
1699
1700	ENDDO ! loop on klon
1701
1702
1703	!--Re-calculation of saturation specific humidity
1704	!--because riming changed temperature
1705	CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
1706
1707	DO i = 1, klon
1708
1709	!---------------------------------------------------------
1710	!-- MELTING
1711	!---------------------------------------------------------
1712	!--Process through which snow melts into rain.
1713	!--The formula is homemade.
1714	!--NB.: this process needs a temperature adjustment
1715
1716	!--dqsmelt_max : maximum snow melting so that temperature
1717	!-- stays higher than 273 K [kg/kg]
1718	!--capa_snowflake : capacitance of a snowflake, equal to
1719	!-- the radius if the snowflake is a sphere [m]
1720	!--temp_wetbulb : wet-bulb temperature [K]
1721	!--snow_fallspeed : snow fall velocity (in clear/cloudy sky) [m/s]
1722	!--air_thermal_conduct : thermal conductivity of the air [J/m/K/s]
1723	!--gamma_melt : tuning parameter for melting [-]
1724	!--nb_snowflake : number of snowflakes (in clear/cloudy air) [-]
1725
1726	IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
1727	!--Computed according to
1728	!--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1729	dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
1730	* ( 1. + RVTMP2 * qtot(i) ))
1731
1732	!--Initialisation
1733	dqsclrmelt = 0.
1734	dqscldmelt = 0.
1735
1736	!--We assume that the snowflakes are spherical
1737	capa_snowflake = r_snow
1738	!--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
1739	air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
1740	!--rho_snow formula follows Brandes et al. 2007 (JAMC)
1741	rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
1742
1743	!--In clear air
1744	IF ( ( snowclr(i) .GT. 0. ) .AND. ( precipfracclr(i) .GT. 0. ) ) THEN
1745	!--Formula for the wet-bulb temperature from ECMWF (IFS)
1746	!--The vapor used is the vapor in the clear sky
1747	temp_wetbulb = temp(i) &
1748	- ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) &
1749	* ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
1750	- 40.637 * ( temp(i) - 275. ) )
1751	!--Calculated according to
1752	!-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
1753	nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr &
1754	/ ( 4. / 3. * RPI * r_snow*3. rho_snow )
1755	dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
1756	* capa_snowflake / RLMLT * gamma_melt &
1757	* MAX(0., temp_wetbulb - RTT) * dtime
1758
1759	!--Barrier to limit the melting flux to the clr snow flux in the mesh
1760	dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i))
1761	ENDIF
1762
1763
1764	!--In cloudy air
1765	IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
1766	!--As the air is saturated, the wet-bulb temperature is equal to the
1767	!--temperature
1768	temp_wetbulb = temp(i)
1769	!--Calculated according to
1770	!-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
1771	nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld &
1772	/ ( 4. / 3. * RPI * r_snow*3. rho_snow )
1773	dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
1774	* capa_snowflake / RLMLT * gamma_melt &
1775	* MAX(0., temp_wetbulb - RTT) * dtime
1776
1777	!--Barrier to limit the melting flux to the cld snow flux in the mesh
1778	dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i))
1779	ENDIF
1780
1781
1782	!--Barrier on temperature. If the total melting flux leads to a
1783	!--positive temperature, it is limited to keep temperature above 0 degC.
1784	!--It is activated if the total is LOWER than the max
1785	!--because everything is negative
1786	dqstotmelt = dqsclrmelt + dqscldmelt
1787	IF ( dqstotmelt .LT. dqsmelt_max ) THEN
1788	!--We redistribute the max melted snow keeping
1789	!--the clear/cloud partition of the melted snow
1790	dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
1791	dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
1792	dqstotmelt = dqsmelt_max
1793
1794	ENDIF
1795
1796	!--Add tendencies
1797	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1798	rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i))
1799	raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i))
1800	snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i))
1801	snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i))
1802
1803	!--Temperature adjustment with the release of latent
1804	!--heat because of melting
1805	temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
1806	/ ( 1. + RVTMP2 * qtot(i) )
1807
1808	!--Diagnostic tendencies
1809	dqrmelt(i) = - dqstotmelt / dtime
1810	dqsmelt(i) = dqstotmelt / dtime
1811
1812	ENDIF
1813
1814
1815	!---------------------------------------------------------
1816	!-- FREEZING
1817	!---------------------------------------------------------
1818	!--Process through which rain freezes into snow.
1819	!-- We parameterize it as a 2 step process:
1820	!--first: freezing following collision with ice crystals
1821	!--second: immersion freezing following (inspired by Bigg 1953)
1822	!--the latter is parameterized as an exponential decrease of the rain
1823	!--water content with a homemade formulya
1824	!--This is based on a caracteritic time of freezing, which
1825	!--exponentially depends on temperature so that it is
1826	!--equal to 1 for temp_nowater (see below) and is close to
1827	!--0 for RTT (=273.15 K).
1828	!--NB.: this process needs a temperature adjustment
1829	!--dqrfreez_max : maximum rain freezing so that temperature
1830	!-- stays lower than 273 K [kg/kg]
1831	!--tau_freez : caracteristic time of freezing [s]
1832	!--gamma_freez : tuning parameter [s-1]
1833	!--alpha_freez : tuning parameter for the shape of the exponential curve [-]
1834	!--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC)
1835
1836	IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
1837
1838
1839	!--1st step: freezing following collision with ice crystals
1840	!--Sub-process of freezing which quantifies the collision between
1841	!--ice crystals in suspension and falling rain droplets.
1842	!--The rain droplets freeze, becoming graupel, and carrying
1843	!--the ice crystal (which acted as an ice nucleating particle).
1844	!--The formula is adapted from the riming formula.
1845	!--it works only in the cloudy part
1846
1847	dqifreez = 0.
1848	dqrtotfreez_step1 = 0.
1849
1850	IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. &
1851	( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
1852	dqrclrfreez = 0.
1853	dqrcldfreez = 0.
1854
1855	!--Computed according to
1856	!--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1857	dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
1858	* ( 1. + RVTMP2 * qtot(i) ))
1859
1860
1861	!--The collision efficiency is assumed unity
1862	Eff_rain_ice = 1.
1863	coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
1864	!--Exact version, which does not need a barrier because of
1865	!--the exponential decrease.
1866	dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
1867
1868	!--We add the part of rain water that freezes, limited by a temperature barrier
1869	!--This quantity is calculated assuming that the number of drop that freeze correspond to the number
1870	!--of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
1871	!--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3
1872	!--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3
1873	!--The assumption above corresponds to dNi = dNr, i.e.,
1874	!-- dqr = dqi * (4/3 PI rho_rain * r_rain*3) / (4/3 PI rho_ice r_ice**3)
1875	!--Dry density [kg/m3]
1876	rho = pplay(i) / temp(i) / RD
1877	!--r_ice formula from Sun and Rikus (1999)
1878	r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 &
1879	+ 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )*0.2535 ( temp(i) - RTT + 190. ) ) / 2.
1880	dqrcldfreez = dqifreez * rho_rain * r_rain*3. / ( rho_ice r_ice**3. )
1881	dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i))
1882	dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max)
1883	dqrtotfreez_step1 = dqrcldfreez
1884
1885	!--Add tendencies
1886	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1887	qice(i) = qice(i) + dqifreez
1888	raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
1889	snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i))
1890	temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD &
1891	/ ( 1. + RVTMP2 * qtot(i) )
1892
1893	ENDIF
1894
1895	!-- Second step immersion freezing of rain drops
1896	!-- with a homemade timeconstant depending on temperature
1897
1898	dqrclrfreez = 0.
1899	dqrcldfreez = 0.
1900	dqrtotfreez_step2 = 0.
1901	!--Computed according to
1902	!--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
1903
1904	dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
1905	* ( 1. + RVTMP2 * qtot(i) ))
1906
1907
1908	tau_freez = 1. / ( beta_freez &
1909	* EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
1910
1911
1912	!--In clear air
1913	IF ( rainclr(i) .GT. 0. ) THEN
1914	!--Exact solution of dqrain/dt = -qrain/tau_freez
1915	dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
1916	ENDIF
1917
1918	!--In cloudy air
1919	IF ( raincld(i) .GT. 0. ) THEN
1920	!--Exact solution of dqrain/dt = -qrain/tau_freez
1921	dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
1922	ENDIF
1923
1924	!--temperature barrier step 2
1925	!--It is activated if the total is LOWER than the max
1926	!--because everything is negative
1927	dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
1928
1929	IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
1930	!--We redistribute the max freezed rain keeping
1931	!--the clear/cloud partition of the freezing rain
1932	dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
1933	dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
1934	dqrtotfreez_step2 = dqrfreez_max
1935	ENDIF
1936
1937
1938	!--Add tendencies
1939	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
1940	rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i))
1941	raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
1942	snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i))
1943	snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i))
1944
1945
1946	!--Temperature adjustment with the uptake of latent
1947	!--heat because of freezing
1948	temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
1949	/ ( 1. + RVTMP2 * qtot(i) )
1950
1951	!--Diagnostic tendencies
1952	dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2
1953	dqrfreez(i) = dqrtotfreez / dtime
1954	dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
1955
1956	ENDIF
1957
1958
1959	!---------------------------------------------------------
1960	!-- ICE SEDIMENTATION
1961	!---------------------------------------------------------
1962	IF ( ok_ice_sedim ) THEN
1963
1964	IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. eps ) ) THEN
1965
1966	rho = pplay(i) / temp(i) / RD
1967	iwcg = MAX( rho * qice(i) / cldfra(i) * 1000., eps ) ! iwc in g/m3
1968	!tempc = temp(i) - RTT ! celcius temp
1969	!! so-called 'empirical parameterization' in Stubenrauch et al. 2019
1970	!IF ( tempc .GE. -60. ) THEN
1971	! icevelo = EXP( 1.9714 + 0.00216078 * tempc &
1972	! - ( 0.0000414122 * tempc*2 + 0.00538922 tempc + 0.0516344 ) * LOG(iwcg) )
1973	!ELSE
1974	! icevelo = 65. * iwcg*0.2 ( 15000. / pplay(i) )**0.15
1975	!ENDIF
1976	!icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s
1977
1978	icevelo = fallice_sedim * 1.645 * ( iwcg / 1000. )**0.16
1979
1980	dz = ( paprsdn(i) - paprsup(i) ) / RG / rho
1981	qice_sedim = qice(i) * MIN(1., icevelo * dtime / dz)
1982	!--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two
1983	!--times in the same timestep)
1984	qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime)
1985	!--Add tendencies
1986	qice(i) = qice(i) - qice_sedim
1987	!--Convert to flux
1988	icesed_flux(i) = qice_sedim * dhum_to_dflux(i)
1989
1990	!--Diagnostic tendencies
1991	dqised(i) = dqised(i) - qice_sedim / dtime
1992	ENDIF
1993	ENDIF
1994
1995
1996	!--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
1997	!--we reduce the precipiration fraction in the clear/cloudy air so that the new
1998	!--local flux of rain+snow is equal to rain_int_min.
1999	!--Here, rain+snow is the gridbox-mean flux of precip.
2000	!--Therefore, (rain+snow)/precipfrac is the local flux of precip.
2001	!--If the local flux of precip is lower than rain_int_min, i.e.,
2002	!-- (rain+snow)/precipfrac < rain_int_min , i.e.,
2003	!-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
2004	!--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
2005	!--Note that this is physically different than what is proposed in LTP thesis.
2006	precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
2007
2008	!--Calculate outputs
2009	rain(i) = rainclr(i) + raincld(i)
2010	snow(i) = snowclr(i) + snowcld(i)
2011
2012	!--Diagnostics
2013	!--BEWARE this is indeed a diagnostic: this is an estimation from
2014	!--the value of the flux at the bottom interface of the mesh and
2015	!--and assuming an upstream numerical calculation
2016	!--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is
2017	!--used for computing the total ice water content in the mesh
2018	!--for radiation only
2019	qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr &
2020	+ raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld )
2021	qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr &
2022	+ snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld )
2023
2024
2025	ENDDO ! loop on klon
2026
2027	END SUBROUTINE poprecip_postcld
2028
2029	END MODULE lmdz_lscp_precip

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

Download in other formats:

Original Format