Context Navigation

lmdz_lscp_condensation.f90 @ 5396

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

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

Download in other formats:

Original Format