Context Navigation

lmdz_lscp_condensation.F90 @ 5041

Last change on this file since 5041 was 5041, checked in by aborella, 2 months ago
Added capa_cond_cirrus tuning parameter
File size: 45.2 KB

Line
1	MODULE lmdz_lscp_condensation
2	!----------------------------------------------------------------
3	! Module for condensation of clouds routines
4	! that are called in LSCP
5
6
7	IMPLICIT NONE
8
9	CONTAINS
10
11	!**********************************************************************************
12	SUBROUTINE condensation_lognormal( &
13	klon, temp, qtot, qsat, gamma_cond, ratqs, &
14	keepgoing, cldfra, qincld, qvc)
15
16	!----------------------------------------------------------------------
17	! This subroutine calculates the formation of clouds, using a
18	! statistical cloud scheme. It uses a generalised lognormal distribution
19	! of total water in the gridbox
20	! See Bony and Emanuel, 2001
21	!----------------------------------------------------------------------
22
23	USE lmdz_lscp_ini, ONLY: eps
24
25	IMPLICIT NONE
26
27	!
28	! Input
29	!
30	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
31	!
32	REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K]
33	REAL, INTENT(IN) , DIMENSION(klon) :: qtot ! total specific humidity (without precip) [kg/kg]
34	REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg]
35	REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-]
36	REAL, INTENT(IN) , DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-]
37	LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed
38	!
39	! Output
40	! NB. those are in INOUT because of the convergence loop on temperature
41	! (in some cases, the values are not re-computed) but the values
42	! are never used explicitely
43	!
44	REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra ! cloud fraction [-]
45	REAL, INTENT(INOUT), DIMENSION(klon) :: qincld ! cloud-mean in-cloud total specific water [kg/kg]
46	REAL, INTENT(INOUT), DIMENSION(klon) :: qvc ! gridbox-mean vapor in the cloud [kg/kg]
47	!
48	! Local
49	!
50	INTEGER :: i
51	REAL :: pdf_std, pdf_k, pdf_delta
52	REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
53	!
54	!--Loop on klon
55	!
56	DO i = 1, klon
57
58	!--If a new calculation of the condensation is needed,
59	!--i.e., temperature has not yet converged (or the cloud is
60	!--formed elsewhere)
61	IF (keepgoing(i)) THEN
62
63	pdf_std = ratqs(i) * qtot(i)
64	pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) )
65	pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
66	pdf_a = pdf_delta / ( pdf_k * SQRT(2.) )
67	pdf_b = pdf_k / (2. * SQRT(2.))
68	pdf_e1 = pdf_a - pdf_b
69	pdf_e1 = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
70	pdf_e1 = 1. - ERF(pdf_e1)
71	pdf_e2 = pdf_a + pdf_b
72	pdf_e2 = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
73	pdf_e2 = 1. - ERF(pdf_e2)
74
75
76	IF ( pdf_e1 .LT. eps ) THEN
77	cldfra(i) = 0.
78	qincld(i) = qsat(i)
79	!--AB grid-mean vapor in the cloud - we assume saturation adjustment
80	qvc(i) = 0.
81	ELSE
82	cldfra(i) = 0.5 * pdf_e1
83	qincld(i) = qtot(i) * pdf_e2 / pdf_e1
84	!--AB grid-mean vapor in the cloud - we assume saturation adjustment
85	qvc(i) = qsat(i) * cldfra(i)
86	ENDIF
87
88	ENDIF ! end keepgoing
89
90	ENDDO ! end loop on i
91	END SUBROUTINE condensation_lognormal
92	!**********************************************************************************
93
94	!**********************************************************************************
95	SUBROUTINE condensation_ice_supersat( &
96	klon, dtime, missing_val, pplay, paprsdn, paprsup, &
97	cf_seri, rvc_seri, ratio_qi_qtot, shear, pbl_eps, cell_area, &
98	temp, qtot, qsat, gamma_cond, ratqs, keepgoing, &
99	cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, &
100	dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
101	Tcontr, qcontr, qcontr2, fcontrN, fcontrP, flight_dist, flight_h2o, &
102	dcf_avi, dqi_avi, dqvc_avi)
103
104	!----------------------------------------------------------------------
105	! This subroutine calculates the formation, evolution and dissipation
106	! of clouds, using a process-oriented treatment of the cloud properties
107	! (cloud fraction, vapor in the cloud, condensed water in the cloud).
108	! It allows for ice supersaturation in cold regions, in clear sky.
109	! If ok_unadjusted_clouds, it also allows for sub- and supersaturated
110	! cloud water vapors.
111	! It also allows for the formation and evolution of condensation trails
112	! (contrails) from aviation.
113	! Authors: Audran Borella, Etienne Vignon, Olivier Boucher
114	! April 2024
115	!----------------------------------------------------------------------
116
117	USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC
118	USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
119	USE lmdz_lscp_ini, ONLY: eps, temp_nowater, ok_weibull_warm_clouds
120	USE lmdz_lscp_ini, ONLY: ok_unadjusted_clouds, iflag_cloud_sublim_pdf
121	USE lmdz_lscp_ini, ONLY: lunout
122
123	USE lmdz_lscp_ini, ONLY: capa_cond_cirrus, mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
124	rhlmid_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, rhl0_pdf_lscp, &
125	cond_thresh_pdf_lscp, coef_mixing_lscp, coef_shear_lscp, &
126	chi_mixing_lscp, rho_ice
127
128	IMPLICIT NONE
129
130	!
131	! Input
132	!
133	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
134	REAL, INTENT(IN) :: dtime ! time step [s]
135	REAL, INTENT(IN) :: missing_val ! missing value for output
136	!
137	REAL, INTENT(IN) , DIMENSION(klon) :: pplay ! layer pressure [Pa]
138	REAL, INTENT(IN) , DIMENSION(klon) :: paprsdn ! pressure at the lower interface [Pa]
139	REAL, INTENT(IN) , DIMENSION(klon) :: paprsup ! pressure at the upper interface [Pa]
140	REAL, INTENT(IN) , DIMENSION(klon) :: cf_seri ! cloud fraction [-]
141	REAL, INTENT(IN) , DIMENSION(klon) :: rvc_seri ! gridbox-mean water vapor in cloud [kg/kg]
142	REAL, INTENT(IN) , DIMENSION(klon) :: ratio_qi_qtot ! specific ice water content to total specific water ratio [-]
143	REAL, INTENT(IN) , DIMENSION(klon) :: shear ! vertical shear [s-1]
144	REAL, INTENT(IN) , DIMENSION(klon) :: pbl_eps !
145	REAL, INTENT(IN) , DIMENSION(klon) :: cell_area !
146	REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K]
147	REAL, INTENT(IN) , DIMENSION(klon) :: qtot ! total specific humidity (without precip) [kg/kg]
148	REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg]
149	REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-]
150	REAL, INTENT(INOUT), DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-]
151	LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed
152	!
153	! Input for aviation
154	!
155	REAL, INTENT(IN), DIMENSION(klon) :: flight_dist !
156	REAL, INTENT(IN), DIMENSION(klon) :: flight_h2o !
157	!
158	! Output
159	! NB. cldfra and qincld should be outputed as cf_seri and qi_seri,
160	! or as tendencies (maybe in the future)
161	! NB. those are in INOUT because of the convergence loop on temperature
162	! (in some cases, the values are not re-computed) but the values
163	! are never used explicitely
164	!
165	REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra ! cloud fraction [-]
166	REAL, INTENT(INOUT), DIMENSION(klon) :: qincld ! cloud-mean in-cloud total specific water [kg/kg]
167	REAL, INTENT(INOUT), DIMENSION(klon) :: qvc ! gridbox-mean vapor in the cloud [kg/kg]
168	REAL, INTENT(INOUT), DIMENSION(klon) :: issrfra ! ISSR fraction [-]
169	REAL, INTENT(INOUT), DIMENSION(klon) :: qissr ! gridbox-mean ISSR specific water [kg/kg]
170	!
171	! Diagnostics for condensation and ice supersaturation
172	! NB. idem for the INOUT
173	!
174	REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_sub ! cloud fraction tendency because of sublimation [s-1]
175	REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_con ! cloud fraction tendency because of condensation [s-1]
176	REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_mix ! cloud fraction tendency because of cloud mixing [s-1]
177	REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_adj ! specific ice content tendency because of temperature adjustment [kg/kg/s]
178	REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_sub ! specific ice content tendency because of sublimation [kg/kg/s]
179	REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_con ! specific ice content tendency because of condensation [kg/kg/s]
180	REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_mix ! specific ice content tendency because of cloud mixing [kg/kg/s]
181	REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
182	REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s]
183	REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s]
184	REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
185	!
186	! Diagnostics for aviation
187	! NB. idem for the INOUT
188	!
189	REAL, INTENT(INOUT), DIMENSION(klon) :: Tcontr ! critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) [K]
190	REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr !
191	REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr2 !
192	REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrN !
193	REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrP !
194	REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_avi ! cloud fraction tendency because of aviation [s-1]
195	REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_avi ! specific ice content tendency because of aviation [kg/kg/s]
196	REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_avi ! specific cloud water vapor tendency because of aviation [kg/kg/s]
197	!
198	! Local
199	!
200	INTEGER :: i
201	LOGICAL :: ok_warm_cloud
202	REAL, DIMENSION(klon) :: qcld, qzero
203	!
204	! for lognormal
205	REAL :: pdf_std, pdf_k, pdf_delta
206	REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
207	!
208	! for unadjusted clouds
209	REAL :: qvapincld, qvapincld_new
210	REAL :: qiceincld, qice_ratio
211	REAL :: pres_sat, kappa
212	REAL :: air_thermal_conduct, water_vapor_diff
213	REAL :: r_ice
214	!
215	! for sublimation
216	REAL :: pdf_alpha
217	REAL :: dqt_sub
218	!
219	! for condensation
220	REAL, DIMENSION(klon) :: qsatl, dqsatl
221	REAL :: clrfra, qclr, sl_clr, rhl_clr
222	REAL :: pdf_ratqs, pdf_skew, pdf_scale, pdf_loc
223	REAL :: pdf_x, pdf_y, pdf_T
224	REAL :: pdf_e3, pdf_e4
225	REAL :: cf_cond, qt_cond, dqt_con
226	!
227	! for mixing
228	REAL, DIMENSION(klon) :: subfra, qsub
229	REAL :: dqt_mix_sub, dqt_mix_issr
230	REAL :: dcf_mix_sub, dcf_mix_issr
231	REAL :: dqvc_mix_sub, dqvc_mix_issr
232	REAL :: dqt_mix
233	REAL :: a_mix, bovera, N_cld_mix, L_mix
234	REAL :: envfra_mix, cldfra_mix
235	REAL :: L_shear, shear_fra
236	REAL :: sigma_mix, issrfra_mix, subfra_mix, qvapinmix
237	!
238	! for cell properties
239	REAL :: rho, rhodz, dz
240	!REAL :: V_cell, M_cell
241	!
242	! for aviation and cell properties
243	!REAL :: dqt_avi
244	!REAL :: contrail_fra
245	!
246	!
247	!--more local variables for diagnostics
248	!--imported from YOMCST.h
249	!--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1)
250	!--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air
251	!--values from Schumann, Meteorol Zeitschrift, 1996
252	!--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
253	!--Qheat = 43. / 50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
254	!REAL, PARAMETER :: EiH2O=1.25 !--emission index of water vapour for kerosene (kg kg-1)
255	!REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1)
256	!REAL, PARAMETER :: eta=0.3 !--average propulsion efficiency of the aircraft
257	!--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
258	!--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
259	!REAL :: Gcontr
260	!--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2)
261	!--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1)
262	!REAL :: qsatliqcontr
263
264
265	!-----------------------------------------------
266	! Initialisations
267	!-----------------------------------------------
268
269	! Ajout des émissions de H2O dues à l'aviation
270	! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
271	! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O )
272	! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
273	! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
274	! flight_h2O is in kg H2O / s / cell
275	!
276	!IF (ok_plane_h2o) THEN
277	! q = ( M_cellq + flight_h2o(i,k)dtime(1.-q) ) / (M_cell + flight_h2o(i,k)dtime*(1.-q) )
278	!ENDIF
279
280
281	qzero(:) = 0.
282
283	!--Calculation of qsat w.r.t. liquid
284	CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsatl)
285
286	!
287	!--Loop on klon
288	!
289	DO i = 1, klon
290
291	!--If a new calculation of the condensation is needed,
292	!--i.e., temperature has not yet converged (or the cloud is
293	!--formed elsewhere)
294	IF (keepgoing(i)) THEN
295
296	!--If the temperature is higher than the threshold below which
297	!--there is no liquid in the gridbox, we activate the usual scheme
298	!--(generalised lognormal from Bony and Emanuel 2001)
299	!--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for
300	!--all clouds, and the lognormal scheme is not activated
301	IF ( ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ) THEN
302
303	pdf_std = ratqs(i) * qtot(i)
304	pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) )
305	pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
306	pdf_a = pdf_delta / ( pdf_k * SQRT(2.) )
307	pdf_b = pdf_k / (2. * SQRT(2.))
308	pdf_e1 = pdf_a - pdf_b
309	pdf_e1 = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
310	pdf_e1 = 1. - ERF(pdf_e1)
311	pdf_e2 = pdf_a + pdf_b
312	pdf_e2 = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
313	pdf_e2 = 1. - ERF(pdf_e2)
314
315
316	IF ( pdf_e1 .LT. eps ) THEN
317	cldfra(i) = 0.
318	qincld(i) = qsat(i)
319	qvc(i) = 0.
320	ELSE
321	cldfra(i) = 0.5 * pdf_e1
322	qincld(i) = qtot(i) * pdf_e2 / pdf_e1
323	qvc(i) = qsat(i) * cldfra(i)
324	ENDIF
325
326	!--If the temperature is lower than temp_nowater, we use the new
327	!--condensation scheme that allows for ice supersaturation
328	ELSE
329
330	!--Initialisation
331	IF ( temp(i) .GT. temp_nowater ) THEN
332	!--If the air mass is warm (liquid water can exist),
333	!--all the memory is lost and the scheme becomes statistical,
334	!--i.e., the sublimation and mixing processes are deactivated,
335	!--and the condensation process is slightly adapted
336	!--This can happen only if ok_weibull_warm_clouds = .TRUE.
337	cldfra(i) = 0.
338	qvc(i) = 0.
339	qcld(i) = 0.
340	ok_warm_cloud = .TRUE.
341	ELSE
342	!--The following barriers ensure that the traced cloud properties
343	!--are consistent. In some rare cases, i.e. the cloud water vapor
344	!--can be greater than the total water in the gridbox
345	cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
346	qcld(i) = MAX(0., MIN(qtot(i), ( ratio_qi_qtot(i) + rvc_seri(i) ) * qtot(i)))
347	qvc(i) = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
348	ok_warm_cloud = .FALSE.
349	ENDIF
350
351	dcf_sub(i) = 0.
352	dqi_sub(i) = 0.
353	dqvc_sub(i) = 0.
354	dqi_adj(i) = 0.
355	dqvc_adj(i) = 0.
356	dcf_con(i) = 0.
357	dqi_con(i) = 0.
358	dqvc_con(i) = 0.
359	dcf_mix(i) = 0.
360	dqi_mix(i) = 0.
361	dqvc_mix(i) = 0.
362
363	issrfra(i) = 0.
364	qissr(i) = 0.
365	subfra(i) = 0.
366	qsub(i) = 0.
367
368	!--Initialisation of the cell properties
369	!--Dry density [kg/m3]
370	rho = pplay(i) / temp(i) / RD
371	!--Dry air mass [kg/m2]
372	rhodz = ( paprsdn(i) - paprsup(i) ) / RG
373	!--Cell thickness [m]
374	dz = rhodz / rho
375	!--Cell volume [m3]
376	!V_cell = dz * cell_area(i)
377	!--Cell dry air mass [kg]
378	!M_cell = rhodz * cell_area(i)
379
380
381	IF ( ok_unadjusted_clouds ) THEN
382	!--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
383	!--hypothesis is lost, and the vapor in the cloud is purely prognostic.
384	!
385	!--The deposition equation is
386	!-- dmi/dt = alpha4piCSvi / ( R_vT/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
387	!--from Lohmann et al. (2016), where
388	!--alpha is the deposition coefficient [-]
389	!--mi is the mass of one ice crystal [kg]
390	!--C is the capacitance of an ice crystal [m]
391	!--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
392	!--R_v is the specific gas constant for humid air [J/kg/K]
393	!--T is the temperature [K]
394	!--esi is the saturation pressure w.r.t. ice [Pa]
395	!--Dv is the diffusivity of water vapor [m2/s]
396	!--Ls is the specific latent heat of sublimation [J/kg/K]
397	!--ka is the thermal conductivity of dry air [J/m/s/K]
398	!
399	!--We fix alpha = 0.5 following Lohmann et al. (2016)
400	!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
401	!-- C = capa_cond_cirrus * r_ice
402	!
403	!--We have qice = Nice * mi, where Nice is the ice crystal
404	!--number concentration per kg of moist air
405	!--HYPOTHESIS 1: the ice crystals are spherical, therefore
406	!-- mi = 4/3 * pi * r_ice*3 rho_ice
407	!--HYPOTHESIS 2: the ice crystals are monodisperse with the
408	!--initial radius r_ice_0.
409	!--NB. this is notably different than the assumption
410	!--of a distributed qice in the cloud made in the sublimation process;
411	!--should it be consistent?
412	!
413	!--As the deposition process does not create new ice crystals,
414	!--and because we assume a same r_ice value for all crystals
415	!--therefore the sublimation process does not destroy ice crystals
416	!--(or, in a limit case, it destroys all ice crystals), then
417	!--Nice is a constant during the sublimation/deposition process.
418	!-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI r_ice_0**3 rho_ice )
419	!
420	!--The deposition equation then reads:
421	!-- dqi/dt = alpha4picapa_cond_cirrusr_ice(qvc-qsat)/qsat / ( R_vT/esi/Dv + Ls/ka/T (Ls/R_v/T - 1) ) * Nice
422	!-- dqi/dt = alpha4picapa_cond_cirrus* (qi / qi_0)*(1/3) r_ice_0*(qvc-qsat)/qsat &
423	!-- / ( R_vT/esi/Dv + Ls/ka/T (Ls*R_v/T - 1) ) &
424	!-- * qi_0 / ( 4/3 RPI r_ice_0**3 rho_ice )
425	!-- dqi/dt = qi*(1/3) (qvc - qsat) * qi_0**(2/3) &
426	!-- alpha/qsatcapa_cond_cirrus/ (R_vT/esi/Dv + Ls/ka/T(LsR_v/T - 1)) / ( 1/3 r_ice_0*2 rho_ice )
427	!--and we have
428	!-- dqvc/dt = - qi*(1/3) (qvc - qsat) / kappa * qi_0(2/3) / r_ice_02
429	!-- dqi/dt = qi*(1/3) (qvc - qsat) / kappa * qi_0(2/3) / r_ice_02
430	!--where kappa = (2/3rho_ice)/capa_cond_cirrusqsat(R_vT/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
431	!--having replaced alpha = 0.5 in the formula
432	!
433	!--This system of equations can be resolved with an exact
434	!--explicit numerical integration, having one variable resolved
435	!--explicitly, the other exactly. The exactly resolved variable is
436	!--the one decreasing, so it is qvc if the process is
437	!--condensation, qi if it is sublimation.
438	!
439	!--kappa is computed as an initialisation constant, as it depends only
440	!--on temperature and other pre-computed values
441	pres_sat = qsat(i) / ( EPS_W + ( 1. - EPS_W ) * qsat(i) ) * pplay(i)
442	!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
443	air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
444	!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
445	water_vapor_diff = 0.211 * ( temp(i) / RTT )*1.94 ( 101325. / pplay(i) ) * 1.e-4
446	kappa = 2. / 3. * rho_ice / capa_cond_cirrus * qsat(i) &
447	* ( RV * temp(i) / water_vapor_diff / pres_sat &
448	+ RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
449	!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
450	ENDIF
451
452
453	!-------------------------------------------------------------------
454	!-- SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD --
455	!-------------------------------------------------------------------
456
457	!--If there is a cloud
458	IF ( cldfra(i) .GT. eps ) THEN
459
460	qvapincld = qvc(i) / cldfra(i)
461	qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
462
463	!--If the ice water content is too low, the cloud is purely sublimated
464	!--Most probably, we advected a cloud with no ice water content (possible
465	!--if the entire cloud precipited for example)
466	IF ( qiceincld .LT. eps ) THEN
467	dcf_sub(i) = - cldfra(i)
468	dqvc_sub(i) = - qvc(i)
469	dqi_sub(i) = - ( qcld(i) - qvc(i) )
470
471	cldfra(i) = 0.
472	qcld(i) = 0.
473	qvc(i) = 0.
474
475	!--Else, the cloud is adjusted and sublimated
476	ELSE
477
478	!--The vapor in cloud cannot be higher than the
479	!--condensation threshold
480	qvapincld = MIN(qvapincld, gamma_cond(i) * qsat(i))
481	qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
482
483	IF ( ok_unadjusted_clouds ) THEN
484	!--Here, the initial vapor in the cloud is qvapincld, and we compute
485	!--the new vapor qvapincld_new
486
487	!--r_ice formula from Sun and Rikus (1999)
488	r_ice = 1.e-6 * ( 45.8966 * ( qiceincld * rho * 1e3 )**0.2214 &
489	+ 0.7957 * ( qiceincld * rho * 1e3 )*0.2535 ( temp(i) - RTT + 190. ) ) / 2.
490
491	IF ( qvapincld .GE. qsat(i) ) THEN
492	!--If the cloud is initially supersaturated
493	!--Exact explicit integration (qvc exact, qice explicit)
494	qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) &
495	* EXP( - dtime * qiceincld / kappa / r_ice**2. )
496	ELSE
497	!--If the cloud is initially subsaturated
498	!--Exact explicit integration (qice exact, qvc explicit)
499	!--The barrier is set so that the resulting vapor in cloud
500	!--cannot be greater than qsat
501	!--qice_ratio is the ratio between the new ice content and
502	!--the old one, it is comprised between 0 and 1
503	qice_ratio = ( 1. - 2. / 3. / kappa / r_ice*2. dtime * ( qsat(i) - qvapincld ) )
504
505	IF ( qice_ratio .LT. 0. ) THEN
506	!--If all the ice has been sublimated, we sublimate
507	!--completely the cloud and do not activate the sublimation
508	!--process
509	!--Tendencies and diagnostics
510	dcf_sub(i) = - cldfra(i)
511	dqvc_sub(i) = - qvc(i)
512	dqi_sub(i) = - ( qcld(i) - qvc(i) )
513
514	cldfra(i) = 0.
515	qcld(i) = 0.
516	qvc(i) = 0.
517
518	!--The new vapor in cloud is set to 0 so that the
519	!--sublimation process does not activate
520	qvapincld_new = 0.
521	ELSE
522	!--Else, the sublimation process is activated with the
523	!--diagnosed new cloud water vapor
524	!--The new vapor in the cloud is increased with the
525	!--sublimated ice
526	qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**(3./2.) )
527	!--The new vapor in the cloud cannot be greater than qsat
528	qvapincld_new = MIN(qvapincld_new, qsat(i))
529	ENDIF ! qice_ratio .LT. 0.
530	ENDIF ! qvapincld .GT. qsat(i)
531	ELSE
532	!--We keep the saturation adjustment hypothesis, and the vapor in the
533	!--cloud is set equal to the saturation vapor
534	qvapincld_new = qsat(i)
535	ENDIF ! ok_unadjusted_clouds
536
537	!--If the vapor in cloud is below vapor needed for the cloud to survive
538	IF ( qvapincld .LT. qvapincld_new ) THEN
539	!--Sublimation of the subsaturated cloud
540	!--iflag_cloud_sublim_pdf selects the PDF of the ice water content
541	!--to use.
542	!--iflag = 1 --> uniform distribution
543	!--iflag = 2 --> exponential distribution
544	!--iflag = 3 --> gamma distribution (Karcher et al 2018)
545
546	IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN
547	!--Uniform distribution starting at qvapincld
548	pdf_e1 = 1. / ( 2. * qiceincld )
549
550	dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
551	dqt_sub = - cldfra(i) * ( qvapincld_new2. - qvapincld2. ) &
552	* pdf_e1 / 2.
553
554	ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN
555	!--Exponential distribution starting at qvapincld
556	pdf_alpha = 1. / qiceincld
557	pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) )
558
559	dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 )
560	dqt_sub = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha &
561	+ qvapincld - qvapincld_new * pdf_e1 )
562
563	ELSEIF ( iflag_cloud_sublim_pdf .GE. 3 ) THEN
564	!--Gamma distribution starting at qvapincld
565	pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
566	pdf_y = pdf_alpha * ( qvapincld_new - qvapincld )
567	pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
568	pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
569
570	dcf_sub(i) = - cldfra(i) * pdf_e1
571	dqt_sub = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
572	ENDIF
573
574	!--Tendencies and diagnostics
575	dqvc_sub(i) = dcf_sub(i) * qvapincld
576	dqi_sub(i) = dqt_sub - dqvc_sub(i)
577
578	!--Add tendencies
579	cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
580	qcld(i) = MAX(0., qcld(i) + dqt_sub)
581	qvc(i) = MAX(0., qvc(i) + dqvc_sub(i))
582
583	ENDIF ! qvapincld .LT. qvapincld_new
584
585	!--Adjustment of the IWC of the remaining cloud
586	!--to the new vapor in cloud (this can be either positive or negative)
587	dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
588	dqi_adj(i) = - dqvc_adj(i)
589
590	!--Add tendencies
591	!--The vapor in the cloud is updated, but not qcld as it is constant
592	!--through this process, as well as cldfra which is unmodified
593	qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
594
595	ENDIF ! qiceincld .LT. eps
596	ENDIF ! cldfra(i) .GT. eps
597
598
599	!--------------------------------------------------------------------------
600	!-- CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS --
601	!--------------------------------------------------------------------------
602	!--This section relies on a distribution of water in the clear-sky region of
603	!--the mesh.
604
605	!--If there is a clear-sky region
606	IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
607
608	!--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
609	qclr = qtot(i) - qcld(i)
610
611	!--New PDF
612	rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.
613	rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
614
615	!--Calculation of the properties of the PDF
616	!--Parameterization from IAGOS observations
617	!--pdf_e1 and pdf_e2 will be reused below
618
619	pdf_e1 = beta_pdf_lscp * cell_area(i)**( 1. / 4. ) &
620	* MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
621	IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
622	pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
623	ELSE
624	pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
625	ENDIF
626	pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
627	pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
628
629	pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
630	pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
631	pdf_loc = rhl_clr - pdf_scale * pdf_e2
632
633	!--Diagnostics of ratqs
634	ratqs(i) = pdf_std / ( qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100. )
635
636	!--Calculation of the newly condensed water and fraction (pronostic)
637	!--Integration of the clear sky PDF between gamma_cond*qsat and +inf
638	!--NB. the calculated values are clear-sky averaged
639
640	pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
641	pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
642	pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
643	pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
644	cf_cond = EXP( - pdf_y )
645	qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.
646
647	IF ( ok_warm_cloud ) THEN
648	!--If the statistical scheme is activated, the first
649	!--calculation of the tendencies is kept for the diagnosis of
650	!--the cloud properties, and the condensation process stops
651	cldfra(i) = cf_cond
652	qcld(i) = qt_cond
653	qvc(i) = cldfra(i) * qsat(i)
654
655	ELSEIF ( cf_cond .GT. cond_thresh_pdf_lscp ) THEN
656	!--We check that there is something to condense, if not the
657	!--condensation process stops
658
659	!--Calculation of the limit qclr in cloud value (stored in pdf_e4)
660
661	pdf_e3 = ( LOG( 1. / cond_thresh_pdf_lscp ) ** ( 1. / pdf_alpha ) - pdf_e2 ) &
662	/ SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. )
663	pdf_e4 = ( 2. * pdf_e3 * pdf_e1 - pdf_x ) &
664	/ ( pdf_e3 * pdf_e1 / rhlmid_pdf_lscp - 1. )
665	IF ( ( MIN(pdf_e4, rhl_clr) .LT. rhlmid_pdf_lscp ) .OR. ( pdf_e4 .GT. rhl_clr ) ) THEN
666	pdf_e4 = pdf_x / ( pdf_e3 * pdf_e1 / rhlmid_pdf_lscp + 1. )
667	ENDIF
668	pdf_e4 = pdf_e4 * qsatl(i) / 100.
669
670	!--Calculation of the final tendencies
671	dcf_con(i) = ( 1. - cldfra(i) ) * ( pdf_e4 - qclr / ( 1. - cldfra(i) ) ) &
672	/ ( pdf_e4 - qt_cond / cf_cond )
673	dqt_con = qclr - pdf_e4 * ( 1. - cldfra(i) - dcf_con(i) )
674
675
676	!--Barriers
677	dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i))
678	dqt_con = MIN(dqt_con, qclr)
679
680
681	IF ( ok_unadjusted_clouds ) THEN
682	!--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
683	!--the new vapor qvapincld. The timestep is divided by two because we do not
684	!--know when the condensation occurs
685	qvapincld = gamma_cond(i) * qsat(i)
686	qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
687	!--r_ice formula from Sun and Rikus (1999)
688	r_ice = 1.e-6 * ( 45.8966 * ( qiceincld * rho * 1e3 )**0.2214 &
689	+ 0.7957 * ( qiceincld * rho * 1e3 )*0.2535 ( temp(i) - RTT + 190. ) ) / 2.
690	!--As qvapincld is necessarily greater than qsat, we only
691	!--use the exact explicit formulation
692	!--Exact explicit version
693	qvapincld = qsat(i) + ( qvapincld - qsat(i) ) &
694	* EXP( - dtime / 2. * qiceincld / kappa / r_ice**2. )
695	ELSE
696	!--We keep the saturation adjustment hypothesis, and the vapor in the
697	!--newly formed cloud is set equal to the saturation vapor.
698	qvapincld = qsat(i)
699	ENDIF
700
701	!--Tendency on cloud vapor and diagnostic
702	dqvc_con(i) = qvapincld * dcf_con(i)
703	dqi_con(i) = dqt_con - dqvc_con(i)
704
705	!--Note that the tendencies are NOT added because they are
706	!--added after the mixing process. In the following, the gridbox fraction is
707	!-- 1. - dcf_con(i), and the total water in the gridbox is
708	!-- qtot(i) - dqi_con(i) - dqvc_con(i)
709
710	ENDIF ! ok_warm_cloud, cf_cond .GT. cond_thresh_pdf_lscp
711	ENDIF ! ( 1. - cldfra(i) ) .GT. eps
712
713	!--If there is still clear sky, we diagnose the ISSR
714	!--We recalculte the PDF properties (after the condensation process)
715	IF ( ( ( 1. - dcf_con(i) - cldfra(i) ) .GT. eps ) .AND. .NOT. ok_warm_cloud ) THEN
716	!--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
717	qclr = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i)
718
719	!--New PDF
720	rhl_clr = qclr / ( 1. - dcf_con(i) - cldfra(i) ) / qsatl(i) * 100.
721	rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
722
723	!--Calculation of the properties of the PDF
724	!--Parameterization from IAGOS observations
725	!--pdf_e1 and pdf_e2 will be reused below
726
727	pdf_e1 = beta_pdf_lscp * cell_area(i)**( 1. / 4. ) &
728	* MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
729	IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
730	pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
731	ELSE
732	pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
733	ENDIF
734	pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
735	pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
736
737	pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
738	pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
739	pdf_loc = rhl_clr - pdf_scale * pdf_e2
740
741	!--We then calculate the part that is greater than qsat
742	!--and consider it supersaturated
743
744	pdf_x = qsat(i) / qsatl(i) * 100.
745	pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
746	pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
747	pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
748	issrfra(i) = EXP( - pdf_y ) * ( 1. - dcf_con(i) - cldfra(i) )
749	qissr(i) = ( pdf_e3 * ( 1. - dcf_con(i) - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
750	ENDIF
751
752	!--Calculation of the subsaturated clear sky fraction and water
753	subfra(i) = 1. - dcf_con(i) - cldfra(i) - issrfra(i)
754	qsub(i) = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i) - qissr(i)
755
756
757	!--------------------------------------
758	!-- CLOUD MIXING --
759	!--------------------------------------
760	!--This process mixes the cloud with its surroundings: the subsaturated clear sky,
761	!--and the supersaturated clear sky. It is activated if the cloud is big enough,
762	!--but does not cover the entire mesh.
763	!
764	IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
765	.AND. .NOT. ok_warm_cloud ) THEN
766
767	!--Initialisation
768	dcf_mix_sub = 0.
769	dqt_mix_sub = 0.
770	dqvc_mix_sub = 0.
771	dcf_mix_issr = 0.
772	dqt_mix_issr = 0.
773	dqvc_mix_issr = 0.
774
775
776	!-- PART 1 - TURBULENT DIFFUSION
777
778	!--Clouds within the mesh are assumed to be ellipses. The length of the
779	!--semi-major axis is a and the length of the semi-minor axis is b.
780	!--N_cld_mix is the number of clouds within the mesh, and
781	!--clouds_perim is the total perimeter of the clouds within the mesh,
782	!--not considering interfaces with other meshes (only the interfaces with clear
783	!--sky are taken into account).
784	!--
785	!--The area of each cloud is A = a * b * RPI,
786	!--and the perimeter of each cloud is
787	!-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
788	!--
789	!--With cell_area the area of the cell, we have:
790	!-- cldfra = A * N_cld_mix / cell_area
791	!-- clouds_perim = P * N_cld_mix
792	!--
793	!--We assume that the ratio between b and a is a function of
794	!--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
795	!--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
796	!--are spherical.
797	!-- b / a = bovera = MAX(0.1, cldfra)
798	bovera = MAX(0.1, cldfra(i))
799	!--The clouds perimeter is imposed using the formula from Morcrette 2012,
800	!--based on observations.
801	!-- clouds_perim_normalized = alpha * cldfra * ( 1. - cldfra ) = clouds_perim / ( norm_constant * cell_area )
802	!--
803	!--With all this, we have
804	!-- a = SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) )
805	!-- P = RPI * a * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
806	!--and therefore, using the perimeter
807	!-- alpha * cldfra * ( 1. - cldfra ) * norm_constant * cell_area
808	!-- = N_cld_mix * RPI &
809	!-- * SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) ) &
810	!-- * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
811	!--and finally
812	N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - dcf_con(i) - cldfra(i) )**2. &
813	* cell_area(i) * ( 1. - dcf_con(i) ) * bovera / RPI &
814	/ ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2.
815	!--where coef_mix_lscp = ( alpha * norm_constant )**2.
816	!--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer
817	!--In particular, it is 0 if cldfra = 1
818	a_mix = SQRT( cell_area(i) * ( 1. - dcf_con(i) ) * cldfra(i) / bovera / N_cld_mix / RPI )
819
820	!--The time required for turbulent diffusion to homogenize a region of size
821	!--L_mix is defined as (L_mix2./tke_dissip)(1./3.) (Pope, 2000; Field et al., 2014)
822	!--We compute L_mix and assume that the cloud is mixed over this length
823	L_mix = SQRT( dtime*3. pbl_eps(i) )
824	!--The mixing length cannot be greater than the semi-minor axis. In this case,
825	!--the entire cloud is mixed.
826	L_mix = MIN(L_mix, a_mix * bovera)
827
828	!--The fraction of clear sky mixed is
829	!-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
830	envfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
831	* ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2. )
832	!--The fraction of cloudy sky mixed is
833	!-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
834	cldfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
835	* ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2. )
836
837
838	!-- PART 2 - SHEARING
839
840	!--The clouds are then sheared. We keep the shape and number
841	!--assumptions from before. The clouds are sheared along their
842	!--semi-major axis (a_mix), on the entire cell heigh dz.
843	!--The increase in size is
844	L_shear = coef_shear_lscp * shear(i) * dz * dtime
845	!--therefore, the total increase in fraction is
846	shear_fra = RPI * L_shear * a_mix * bovera * N_cld_mix &
847	/ cell_area(i) / ( 1. - dcf_con(i) )
848	!--and the environment and cloud mixed fractions are the same,
849	!--which we add to the previous calculated mixed fractions.
850	!--We therefore assume that the sheared clouds and the turbulent
851	!--mixed clouds are different.
852	envfra_mix = envfra_mix + shear_fra
853	cldfra_mix = cldfra_mix + shear_fra
854
855
856	!-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
857
858	!--The environment fraction is allocated to subsaturated sky or supersaturated sky,
859	!--according to the factor sigma_mix. This is computed as the ratio of the
860	!--subsaturated sky fraction to the environment fraction, corrected by a factor
861	!--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
862	!--supersaturated sky is favoured. Physically, this means that it is more likely
863	!--to have supersaturated sky around the cloud than subsaturated sky.
864	sigma_mix = subfra(i) / ( subfra(i) + chi_mixing_lscp * issrfra(i) )
865	subfra_mix = MIN( sigma_mix * envfra_mix, subfra(i) )
866	issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra(i) )
867	cldfra_mix = MIN( cldfra_mix, cldfra(i) )
868
869	!--First, we mix the subsaturated sky (subfra_mix) and the cloud close
870	!--to this fraction (sigma_mix * cldfra_mix).
871	IF ( subfra(i) .GT. eps ) THEN
872
873	IF ( ok_unadjusted_clouds ) THEN
874	!--The subsaturated air is simply added to the cloud,
875	!--with the corresponding cloud fraction
876	!--If the cloud is too subsaturated, the sublimation process
877	!--activated in the following timestep will reduce the cloud
878	!--fraction
879	dcf_mix_sub = subfra_mix
880	dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
881	dqvc_mix_sub = dqt_mix_sub
882
883	ELSE
884	!--We compute the total humidity in the mixed air, which
885	!--can be either sub- or supersaturated.
886	qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
887	+ qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
888	/ ( subfra_mix + cldfra_mix * sigma_mix )
889
890	IF ( qvapinmix .GT. qsat(i) ) THEN
891	!--If the mixed air is supersaturated, we condense the subsaturated
892	!--region which was mixed.
893	dcf_mix_sub = subfra_mix
894	dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
895	dqvc_mix_sub = dcf_mix_sub * qsat(i)
896	ELSE
897	!--Else, we sublimate the cloud which was mixed.
898	dcf_mix_sub = - sigma_mix * cldfra_mix
899	dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
900	dqvc_mix_sub = dcf_mix_sub * qsat(i)
901	ENDIF
902	ENDIF ! ok_unadjusted_clouds
903	ENDIF ! subfra .GT. eps
904
905	!--We then mix the supersaturated sky (issrfra_mix) and the cloud,
906	!--for which the mixed air is always supersatured, therefore
907	!--the cloud necessarily expands
908	IF ( issrfra(i) .GT. eps ) THEN
909
910	IF ( ok_unadjusted_clouds ) THEN
911	!--The ice supersaturated air is simply added to the
912	!--cloud, and supersaturated vapor will be deposited on the
913	!--cloud ice crystals by the deposition process in the
914	!--following timestep
915	dcf_mix_issr = issrfra_mix
916	dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
917	dqvc_mix_issr = dqt_mix_issr
918	ELSE
919	!--In this case, the additionnal vapor condenses
920	dcf_mix_issr = issrfra_mix
921	dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
922	dqvc_mix_issr = dcf_mix_issr * qsat(i)
923	ENDIF ! ok_unadjusted_clouds
924
925
926	ENDIF ! issrfra .GT. eps
927
928	!--Sum up the tendencies from subsaturated sky and supersaturated sky
929	dcf_mix(i) = dcf_mix_sub + dcf_mix_issr
930	dqt_mix = dqt_mix_sub + dqt_mix_issr
931	dqvc_mix(i) = dqvc_mix_sub + dqvc_mix_issr
932	dqi_mix(i) = dqt_mix - dqvc_mix(i)
933
934	!--Add tendencies
935	issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
936	qissr(i) = MAX(0., qissr(i) - dqt_mix_issr)
937	cldfra(i) = MAX(0., MIN(1. - dcf_con(i), cldfra(i) + dcf_mix(i)))
938	qcld(i) = MAX(0., MIN(qtot(i) - dqi_con(i) - dqvc_con(i), qcld(i) + dqt_mix))
939	qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))
940
941	ENDIF ! ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
942
943	!--Finally, we add the tendencies of condensation
944	cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
945	qcld(i) = MIN(qtot(i), qcld(i) + dqvc_con(i) + dqi_con(i))
946	qvc(i) = MIN(qcld(i), qvc(i) + dqvc_con(i))
947
948
949	!----------------------------------------
950	!-- CONTRAILS AND AVIATION --
951	!----------------------------------------
952
953	!--Add a source of cirrus from aviation contrails
954	!IF ( ok_plane_contrail ) THEN
955	! dcf_avi(i) = 0.
956	! dqi_avi(i) = 0.
957	! dqvc_avi(i) = 0.
958	! ! TODO implement ok_unadjusted_clouds
959	! IF ( issrfra(i) .GT. eps ) THEN
960	! contrail_fra = MIN(1., flight_m(i,k) * dtime * contrail_cross_section / V_cell)
961	! dcf_avi(i) = issrfra(i) * contrail_fra
962	! dqt_avi = dcf_avi(i) * qissr(i) / issrfra(i)
963	! dqvc_avi(i) = qsat(i) * dcf_avi(i)
964	!
965	! !--Add tendencies
966	! cldfra(i) = cldfra(i) + dcf_avi(i)
967	! issrfra(i) = issrfra(i) - dcf_avi(i)
968	! qcld(i) = qcld(i) + dqt_avi
969	! qvc(i) = qvc(i) + dqvc_avi(i)
970	! qissr(i) = qissr(i) - dqt_avi
971
972	! !--Diagnostics
973	! dqi_avi(i) = dqt_avi - qsat(i) * dcf_avi(i)
974	! ENDIF
975	! dcf_avi(i) = dcf_avi(i) / dtime
976	! dqi_avi(i) = dqi_avi(i) / dtime
977	! dqvc_avi(i) = dqvc_avi(i) / dtime
978	!ENDIF
979
980
981
982	!-------------------------------------------
983	!-- FINAL BARRIERS AND OUTPUTS --
984	!-------------------------------------------
985
986	IF ( cldfra(i) .LT. eps ) THEN
987	!--If the cloud is too small, it is sublimated.
988	cldfra(i) = 0.
989	qcld(i) = 0.
990	qvc(i) = 0.
991	qincld(i) = qsat(i)
992	ELSE
993	qincld(i) = qcld(i) / cldfra(i)
994	ENDIF ! cldfra .LT. eps
995
996	!--Diagnostics
997	dcf_sub(i) = dcf_sub(i) / dtime
998	dcf_con(i) = dcf_con(i) / dtime
999	dcf_mix(i) = dcf_mix(i) / dtime
1000	dqi_adj(i) = dqi_adj(i) / dtime
1001	dqi_sub(i) = dqi_sub(i) / dtime
1002	dqi_con(i) = dqi_con(i) / dtime
1003	dqi_mix(i) = dqi_mix(i) / dtime
1004	dqvc_adj(i) = dqvc_adj(i) / dtime
1005	dqvc_sub(i) = dqvc_sub(i) / dtime
1006	dqvc_con(i) = dqvc_con(i) / dtime
1007	dqvc_mix(i) = dqvc_mix(i) / dtime
1008
1009	ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
1010
1011	ENDIF ! end keepgoing
1012
1013	ENDDO ! end loop on i
1014
1015	END SUBROUTINE condensation_ice_supersat
1016	!**********************************************************************************
1017
1018	END MODULE lmdz_lscp_condensation

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

Download in other formats:

Original Format