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