Context Navigation

lmdz_lscp_condensation.F90 @ 5213

Last change on this file since 5213 was 5210, checked in by evignon, 24 hours ago

changes in ice radius treatment in lmdz_lscp_condebsation + minor changes

Borella

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

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

Download in other formats:

Original Format