Context Navigation

lmdz_lscp_precip.f90 @ 5624

Last change on this file since 5624 was 5614, checked in by evignon, 7 months ago

Commission liée à un update majeur de la routine de condensation grande echelle suite au travail
de Lea, Audran et Etienne
Elle inclue une restructuration des routines pour clarifier le role "moniteur" de la routine lscp_main,
une mise à jour de la parametrisation de partitionnement de phase de Lea pour inclure les nuages de couche limite,
ainsi que des corrections des routines de precipitations "poprecip".
Convergence numerique verifiee en prod et debug pour les physiques NPv6.3 et 7.0.1c

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

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

Download in other formats:

Original Format