source: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90 @ 5679

Last change on this file since 5679 was 5643, checked in by aborella, 3 months ago

Various minor corrections to ISSR and contrails param
File size: 75.3 KB
Line 
1MODULE lmdz_lscp_condensation
2!----------------------------------------------------------------
3! Module for condensation of clouds routines
4! that are called in LSCP
5
6
7IMPLICIT NONE
8
9CONTAINS
10
11!**********************************************************************************
12SUBROUTINE 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
23USE lmdz_lscp_ini, ONLY: eps
24
25IMPLICIT NONE
26
27!
28! Input
29!
30INTEGER,  INTENT(IN)                     :: klon          ! number of horizontal grid points
31!
32REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
33REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot          ! total specific humidity (without precip) [kg/kg]
34REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
35REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
36REAL,     INTENT(IN)   , DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
37LOGICAL,  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!
44REAL,     INTENT(INOUT), DIMENSION(klon) :: cldfra   ! cloud fraction [-]
45REAL,     INTENT(INOUT), DIMENSION(klon) :: qincld   ! cloud-mean in-cloud total specific water [kg/kg]
46REAL,     INTENT(INOUT), DIMENSION(klon) :: qvc      ! gridbox-mean vapor in the cloud [kg/kg]
47!
48! Local
49!
50INTEGER :: i
51REAL :: pdf_std, pdf_k, pdf_delta
52REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
53!
54!--Loop on klon
55!
56DO 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
90ENDDO     ! end loop on i
91END SUBROUTINE condensation_lognormal
92!**********************************************************************************
93
94!**********************************************************************************
95SUBROUTINE condensation_ice_supersat( &
96      klon, dtime, pplay, paprsdn, paprsup, totfra_in, &
97      cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, &
98      temp, qtot_in, qsat, gamma_cond, ratqs, keepgoing, pt_pron_clds, &
99      cldfra_above, icesed_flux, &
100      cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, dcf_sed, &
101      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqi_sed, &
102      dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, dqvc_sed, &
103      lincontfra_in, circontfra_in, qtl_in, qtc_in, flight_dist, flight_h2o, &
104      lincontfra, circontfra, qlincont, qcircont, &
105      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
106      dcfl_ini, dqil_ini, dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub, &
107      dcfl_cir, dqtl_cir, dcfl_mix, dqil_mix, dqtl_mix, &
108      dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix)
109
110!----------------------------------------------------------------------
111! This subroutine calculates the formation, evolution and dissipation
112! of clouds, using a process-oriented treatment of the cloud properties
113! (cloud fraction, vapor in the cloud, condensed water in the cloud).
114! It allows for ice supersaturation in cold regions, in clear sky.
115! If ok_unadjusted_clouds, it also allows for sub- and supersaturated
116! cloud water vapors.
117! It also allows for the formation and evolution of condensation trails
118! (contrails) from aviation.
119! Authors: Audran Borella, Etienne Vignon, Olivier Boucher
120! April 2024
121!----------------------------------------------------------------------
122
123USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC
124USE lmdz_lscp_ini,   ONLY: RLSTT, RTT, RD, RG, RV, RPI, EPS_W
125USE lmdz_lscp_ini,   ONLY: eps, temp_nowater, ok_unadjusted_clouds
126
127USE lmdz_lscp_ini,   ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
128USE lmdz_lscp_ini,   ONLY: N_ice_volume, corr_incld_depsub, nu_iwc_pdf_lscp
129USE lmdz_lscp_ini,   ONLY: beta_pdf_lscp, temp_thresh_pdf_lscp
130USE lmdz_lscp_ini,   ONLY: std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp
131USE lmdz_lscp_ini,   ONLY: coef_mixing_lscp, coef_shear_lscp
132USE lmdz_lscp_ini,   ONLY: aspect_ratio_cirrus, cooling_rate_ice_thresh
133
134USE lmdz_lscp_ini,   ONLY: ok_plane_contrail, aspect_ratio_lincontrails
135USE lmdz_lscp_ini,   ONLY: coef_mixing_lincontrails, coef_shear_lincontrails
136USE lmdz_lscp_ini,   ONLY: chi_mixing_lincontrails, linear_contrails_lifetime
137USE lmdz_aviation,   ONLY: contrails_formation
138
139IMPLICIT NONE
140
141!
142! Input
143!
144INTEGER,  INTENT(IN)                     :: klon          ! number of horizontal grid points
145REAL,     INTENT(IN)                     :: dtime         ! time step [s]
146!
147REAL,     INTENT(IN)   , DIMENSION(klon) :: pplay         ! layer pressure [Pa]
148REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsdn       ! pressure at the lower interface [Pa]
149REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsup       ! pressure at the upper interface [Pa]
150REAL,     INTENT(IN)   , DIMENSION(klon) :: totfra_in     ! total available fraction for stratiform clouds [-]
151REAL,     INTENT(IN)   , DIMENSION(klon) :: cldfra_in     ! cloud fraction [-]
152REAL,     INTENT(IN)   , DIMENSION(klon) :: qvc_in        ! gridbox-mean water vapor in cloud [kg/kg]
153REAL,     INTENT(IN)   , DIMENSION(klon) :: qliq_in       ! specific liquid water content [kg/kg]
154REAL,     INTENT(IN)   , DIMENSION(klon) :: qice_in       ! specific ice water content [kg/kg]
155REAL,     INTENT(IN)   , DIMENSION(klon) :: shear         ! vertical shear [s-1]
156REAL,     INTENT(IN)   , DIMENSION(klon) :: pbl_eps       ! TKE dissipation [m2/s3]
157REAL,     INTENT(IN)   , DIMENSION(klon) :: cell_area     ! cell area [m2]
158REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
159REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot_in       ! total specific humidity (without precip) [kg/kg]
160REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
161REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
162REAL,     INTENT(IN)   , DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
163LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
164LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: pt_pron_clds  ! .TRUE. if clouds are prognostic in this mesh
165REAL,     INTENT(IN)   , DIMENSION(klon) :: cldfra_above  ! cloud fraction IN THE LAYER ABOVE [-]
166REAL,     INTENT(IN)   , DIMENSION(klon) :: icesed_flux   ! sedimentated ice flux [kg/s/m2]
167!
168! Input for aviation
169!
170REAL,     INTENT(IN)   , DIMENSION(klon) :: lincontfra_in ! input linear contrails fraction [-]
171REAL,     INTENT(IN)   , DIMENSION(klon) :: circontfra_in ! input contrail cirrus fraction [-]
172REAL,     INTENT(IN)   , DIMENSION(klon) :: qtl_in        ! input linear contrails total specific humidity [kg/kg]
173REAL,     INTENT(IN)   , DIMENSION(klon) :: qtc_in        ! input contrail cirrus total specific humidity [kg/kg]
174REAL,     INTENT(IN)   , DIMENSION(klon) :: flight_dist   ! aviation distance flown concentration [m/s/m3]
175REAL,     INTENT(IN)   , DIMENSION(klon) :: flight_h2o    ! aviation emitted H2O concentration [kgH2O/s/m3]
176!
177!  Output
178!  NB. cldfra and qincld should be outputed as cf_seri and qi_seri,
179!      or as tendencies (maybe in the future)
180!  NB. those are in INOUT because of the convergence loop on temperature
181!      (in some cases, the values are not re-computed) but the values
182!      are never used explicitely
183!
184REAL,     INTENT(INOUT), DIMENSION(klon) :: cldfra   ! cloud fraction [-]
185REAL,     INTENT(INOUT), DIMENSION(klon) :: qincld   ! cloud-mean in-cloud total specific water [kg/kg]
186REAL,     INTENT(INOUT), DIMENSION(klon) :: qvc      ! gridbox-mean vapor in the cloud [kg/kg]
187REAL,     INTENT(INOUT), DIMENSION(klon) :: issrfra  ! ISSR fraction [-]
188REAL,     INTENT(INOUT), DIMENSION(klon) :: qissr    ! gridbox-mean ISSR specific water [kg/kg]
189!
190!  Diagnostics for condensation and ice supersaturation
191!  NB. idem for the INOUT
192!
193REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_sub  ! cloud fraction tendency because of sublimation [s-1]
194REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_con  ! cloud fraction tendency because of condensation [s-1]
195REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_mix  ! cloud fraction tendency because of cloud mixing [s-1]
196REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_sed  ! cloud fraction tendency because of sedimentation [s-1]
197REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_adj  ! specific ice content tendency because of temperature adjustment [kg/kg/s]
198REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_sub  ! specific ice content tendency because of sublimation [kg/kg/s]
199REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_con  ! specific ice content tendency because of condensation [kg/kg/s]
200REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_mix  ! specific ice content tendency because of cloud mixing [kg/kg/s]
201REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_sed  ! specific ice content tendency because of sedimentation [kg/kg/s]
202REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
203REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s]
204REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s]
205REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
206REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_sed ! specific cloud water vapor tendency because of sedimentation [kg/kg/s]
207!
208!  Diagnostics for aviation
209!  NB. idem for the INOUT
210!
211REAL,     INTENT(INOUT), DIMENSION(klon) :: lincontfra   ! linear contrail fraction [-]
212REAL,     INTENT(INOUT), DIMENSION(klon) :: circontfra   ! contrail cirrus fraction [-]
213REAL,     INTENT(INOUT), DIMENSION(klon) :: qlincont     ! linear contrail specific humidity [kg/kg]
214REAL,     INTENT(INOUT), DIMENSION(klon) :: qcircont     ! contrail cirrus specific humidity [kg/kg]
215REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
216REAL,     INTENT(INOUT), DIMENSION(klon) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
217REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraP  ! potential persistent contrail fraction [-]
218REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-]
219REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfl_ini     ! linear contrails cloud fraction tendency because of initial formation [s-1]
220REAL,     INTENT(INOUT), DIMENSION(klon) :: dqil_ini     ! linear contrails ice specific humidity tendency because of initial formation [kg/kg/s]
221REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtl_ini     ! linear contrails total specific humidity tendency because of initial formation [kg/kg/s]
222REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfl_sub     ! linear contrails cloud fraction tendency because of sublimation [s-1]
223REAL,     INTENT(INOUT), DIMENSION(klon) :: dqil_sub     ! linear contrails ice specific humidity tendency because of sublimation [kg/kg/s]
224REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtl_sub     ! linear contrails total specific humidity tendency because of sublimation [kg/kg/s]
225REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfl_cir     ! linear contrails cloud fraction tendency because of conversion in cirrus [s-1]
226REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtl_cir     ! linear contrails total specific humidity tendency because of conversion in cirrus [kg/kg/s]
227REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfl_mix     ! linear contrails cloud fraction tendency because of mixing [s-1]
228REAL,     INTENT(INOUT), DIMENSION(klon) :: dqil_mix     ! linear contrails ice specific humidity tendency because of mixing [kg/kg/s]
229REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtl_mix     ! linear contrails total specific humidity tendency because of mixing [kg/kg/s]
230REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_sub     ! contrail cirrus cloud fraction tendency because of sublimation [s-1]
231REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_sub     ! contrail cirrus ice specific humidity tendency because of sublimation [kg/kg/s]
232REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_sub     ! contrail cirrus total specific humidity tendency because of sublimation [kg/kg/s]
233REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_mix     ! contrail cirrus cloud fraction tendency because of mixing [s-1]
234REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_mix     ! contrail cirrus ice specific humidity tendency because of mixing [kg/kg/s]
235REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_mix     ! contrail cirrus total specific humidity tendency because of mixing [kg/kg/s]
236!
237! Local
238!
239INTEGER :: i
240LOGICAL :: ok_warm_cloud
241REAL, DIMENSION(klon) :: qcld, qzero
242REAL, DIMENSION(klon) :: clrfra, qclr
243!
244! for lognormal
245REAL :: pdf_std, pdf_k, pdf_delta
246REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
247!
248! for unadjusted clouds
249REAL :: qiceincld, qvapincld, qvapincld_new
250REAL :: qice_ratio
251!
252! for deposition / sublimation
253REAL :: pres_sat, kappa_depsub, tauinv_depsub
254REAL :: air_thermal_conduct, water_vapor_diff
255!
256! for dissipation
257REAL, DIMENSION(klon) :: temp_diss, qsati_diss
258REAL :: pdf_shape, qiceincld_min
259!
260! for condensation
261REAL, DIMENSION(klon) :: qsatl, dqsat_tmp
262REAL, DIMENSION(klon) :: pdf_alpha, pdf_scale, pdf_gamma
263REAL :: rhl_clr, pdf_loc
264REAL :: pdf_e3, pdf_x, pdf_y
265REAL :: dqt_con
266!
267! for sedimentation
268REAL, DIMENSION(klon) :: qice_sedim
269REAL :: clrfra_sed
270!
271! for mixing
272REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix
273REAL :: cldfra_mix, clrfra_mix, sigma_mix
274REAL :: L_shear, shear_fra
275REAL :: qvapinmix, qiceinmix, qvapinmix_lim, qvapinclr_lim
276REAL :: pdf_fra_above_nuc, pdf_q_above_nuc
277REAL :: pdf_fra_above_lim, pdf_q_above_lim
278REAL :: pdf_fra_below_lim
279REAL :: mixed_fraction
280!
281! for cell properties
282REAL :: rho, rhodz, dz
283!
284! for contrails
285REAL :: contrails_conversion_factor
286
287qzero(:) = 0.
288
289!--Calculation of qsat w.r.t. liquid
290CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsat_tmp)
291!--Calculation of qsat max for dissipation
292temp_diss(:) = temp(:) + cooling_rate_ice_thresh * dtime
293CALL calc_qsat_ecmwf(klon, temp_diss, qzero, pplay, RTT, 2, .FALSE., qsati_diss, dqsat_tmp)
294
295!
296!--Loop on klon
297!
298DO i = 1, klon
299
300  !--If a new calculation of the condensation is needed,
301  !--i.e., temperature has not yet converged (or the cloud is
302  !--formed elsewhere)
303  IF (keepgoing(i)) THEN
304
305    !--If the temperature is higher than the threshold below which
306    !--there is no liquid in the gridbox, we activate the usual scheme
307    !--(generalised lognormal from Bony and Emanuel 2001)
308    !--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for
309    !--all clouds, and the lognormal scheme is not activated
310    IF ( .NOT. pt_pron_clds(i) ) THEN
311
312      pdf_std   = ratqs(i) * qtot_in(i)
313      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot_in(i))**2 ) )
314      pdf_delta = LOG( qtot_in(i) / ( gamma_cond(i) * qsat(i) ) )
315      pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
316      pdf_b     = pdf_k / (2. * SQRT(2.))
317      pdf_e1    = pdf_a - pdf_b
318      pdf_e1    = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
319      pdf_e1    = 1. - ERF(pdf_e1)
320      pdf_e2    = pdf_a + pdf_b
321      pdf_e2    = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
322      pdf_e2    = 1. - ERF(pdf_e2)
323
324
325      IF ( pdf_e1 .LT. eps ) THEN
326        cldfra(i) = 0.
327        qincld(i) = qsat(i)
328        qvc(i) = 0.
329      ELSE
330        cldfra(i) = 0.5 * pdf_e1
331        qincld(i) = qtot_in(i) * pdf_e2 / pdf_e1
332        qvc(i) = qsat(i) * cldfra(i)
333      ENDIF
334
335    !--If the temperature is lower than temp_nowater, we use the new
336    !--condensation scheme that allows for ice supersaturation
337    ELSE
338
339      !--Initialisation
340      !--If the air mass is warm (liquid water can exist),
341      !--all the memory is lost and the scheme becomes statistical,
342      !--i.e., the sublimation and mixing processes are deactivated,
343      !--and the condensation process is slightly adapted
344      !--This can happen only if ok_weibull_warm_clouds = .TRUE.
345      ok_warm_cloud = ( temp(i) .GT. temp_nowater )
346
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(totfra_in(i), cldfra_in(i)))
351      qcld(i)   = MAX(0., MIN(qtot_in(i), qliq_in(i) + qice_in(i) + qvc_in(i)))
352      qvc(i)    = MAX(0., MIN(qcld(i), qvc_in(i)))
353
354      !--Initialise clear fraction properties
355      clrfra(i) = totfra_in(i) - cldfra(i)
356      qclr(i) = qtot_in(i) - qcld(i)
357
358      dcf_sub(i)  = 0.
359      dqi_sub(i)  = 0.
360      dqvc_sub(i) = 0.
361      dqi_adj(i)  = 0.
362      dqvc_adj(i) = 0.
363      dcf_con(i)  = 0.
364      dqi_con(i)  = 0.
365      dqvc_con(i) = 0.
366      dcf_mix(i)  = 0.
367      dqi_mix(i)  = 0.
368      dqvc_mix(i) = 0.
369      dcf_sed(i)  = 0.
370      dqi_sed(i)  = 0.
371      dqvc_sed(i) = 0.
372
373      IF ( icesed_flux(i) .GT. 0. ) THEN
374        !--If ice sedimentation is activated, the quantity of sedimentated ice was added
375        !--to the total water vapor in the precipitation routine. Here we remove it
376        !--(it will be reincluded later)
377        qice_sedim(i) = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
378        qclr(i) = qclr(i) - qice_sedim(i)
379      ENDIF
380
381      !--Initialisation of the cell properties
382      !--Dry density [kg/m3]
383      rho = pplay(i) / temp(i) / RD
384      !--Dry air mass [kg/m2]
385      rhodz = ( paprsdn(i) - paprsup(i) ) / RG
386      !--Cell thickness [m]
387      dz = rhodz / rho
388
389      !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
390      !--hypothesis is lost, and the vapor in the cloud is purely prognostic.
391      !
392      !--The deposition equation is
393      !-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
394      !--from Lohmann et al. (2016), where
395      !--alpha is the deposition coefficient [-]
396      !--mi is the mass of one ice crystal [kg]
397      !--C is the capacitance of an ice crystal [m]
398      !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
399      !--R_v is the specific gas constant for humid air [J/kg/K]
400      !--T is the temperature [K]
401      !--esi is the saturation pressure w.r.t. ice [Pa]
402      !--Dv is the diffusivity of water vapor [m2/s]
403      !--Ls is the specific latent heat of sublimation [J/kg/K]
404      !--ka is the thermal conductivity of dry air [J/m/s/K]
405      !
406      !--alpha is a coefficient to take into account the fact that during deposition, a water
407      !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
408      !--coherent (with the same structure). It has no impact for sublimation.
409      !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
410      !--during deposition, and alpha = 1. during sublimation.
411      !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
412      !-- C = capa_cond_cirrus * rm_ice
413      !
414      !--We have qice = Nice * mi, where Nice is the ice crystal
415      !--number concentration per kg of moist air
416      !--HYPOTHESIS 1: the ice crystals are spherical, therefore
417      !-- mi = 4/3 * pi * rm_ice**3 * rho_ice
418      !--HYPOTHESIS 2: the ice crystals concentration is constant in the cloud
419      !
420      !--The equation in terms of q_ice is valide locally, and the local ice water content
421      !--follows a Gamma distribution with a factor nu_iwc_pdf_lscp. Therefore, by
422      !--integrating the local equation over the PDF (entire cloud), a correcting factor
423      !--must be included, equal to
424      !-- corr_incld_depsub = GAMMA(nu + 1/3) / GAMMA(nu) / nu**(1/3)
425      !--NB. this is equal to about 0.9, hence the correction is not big
426      !--NB. to lighten the calculated, corr_incld_depsub is calculated in lmdz_lscp_ini
427      !
428      !--As the deposition process does not create new ice crystals,
429      !--and because we assume a same rm_ice value for all crystals
430      !--therefore the sublimation process does not destroy ice crystals
431      !--(or, in a limit case, it destroys all ice crystals), then
432      !--Nice is a constant during the sublimation/deposition process
433      !--hence dmi = dqi
434      !
435      !--The deposition equation then reads for qi the in-cloud ice water content:
436      !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat * corr_incld_depsub &
437      !--            / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
438      !-- dqi/dt = alpha*4pi*capa_cond_cirrus*Nice*corr_incld_depsub &
439      !--             / ( 4pi/3 N_ice rho_ice )**(1/3) &
440      !--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
441      !--             qi**(1/3) * (qvc - qsat) / qsat
442      !--and we have
443      !-- dqvc/dt = - alpha * kappa(T) * qi**(1/3) * (qvc - qsat)
444      !-- dqi/dt  =   alpha * kappa(T) * qi**(1/3) * (qvc - qsat)
445      !
446      !--This system of equations can be resolved with an exact
447      !--explicit numerical integration, having one variable resolved
448      !--explicitly, the other exactly. qvc is always the variable solved exactly.
449      !
450      !--kappa is computed as an initialisation constant, as it depends only
451      !--on temperature and other pre-computed values
452      pres_sat = qsat(i) / ( EPS_W + ( 1. - EPS_W ) * qsat(i) ) * pplay(i)
453      !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
454      air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
455      !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
456      water_vapor_diff = 0.211 * ( temp(i) / RTT )**1.94 * ( 101325. / pplay(i) ) * 1.e-4
457      !--NB. the greater kappa_depsub, the more efficient is the
458      !--deposition/sublimation process
459      kappa_depsub = 4. * RPI * capa_cond_cirrus * N_ice_volume / rho * corr_incld_depsub &
460                   / qsat(i) / ( 4. / 3. * RPI * N_ice_volume / rho * rho_ice )**(1./3.) &
461                   / ( RV * temp(i) / water_vapor_diff / pres_sat &
462                     + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
463
464      !--If contrails are activated
465      IF ( ok_plane_contrail ) THEN
466        lincontfra(i) = MAX(0., lincontfra_in(i))
467        circontfra(i) = MAX(0., circontfra_in(i))
468        qlincont(i) = MAX(0., qtl_in(i))
469        qcircont(i) = MAX(0., qtc_in(i))
470        !--The following barriers are needed since the advection scheme does not
471        !--conserve order relations
472        mixed_fraction = lincontfra(i) + circontfra(i)
473        IF ( mixed_fraction .GT. cldfra(i) ) THEN
474          mixed_fraction = cldfra(i) / mixed_fraction
475          lincontfra(i) = lincontfra(i) * mixed_fraction
476          circontfra(i) = circontfra(i) * mixed_fraction
477          qlincont(i)   = qlincont(i)   * mixed_fraction
478          qcircont(i)   = qcircont(i)   * mixed_fraction
479        ENDIF
480        mixed_fraction = qlincont(i) + qcircont(i)
481        IF ( mixed_fraction .GT. qcld(i) ) THEN
482          mixed_fraction = qcld(i) / mixed_fraction
483          lincontfra(i) = lincontfra(i) * mixed_fraction
484          circontfra(i) = circontfra(i) * mixed_fraction
485          qlincont(i)   = qlincont(i)   * mixed_fraction
486          qcircont(i)   = qcircont(i)   * mixed_fraction
487        ENDIF
488
489
490        dcfl_ini(i) = 0.
491        dqil_ini(i) = 0.
492        dqtl_ini(i) = 0.
493        dcfl_sub(i) = 0.
494        dqil_sub(i) = 0.
495        dqtl_sub(i) = 0.
496        dcfl_cir(i) = 0.
497        dqtl_cir(i) = 0.
498        dcfl_mix(i) = 0.
499        dqil_mix(i) = 0.
500        dqtl_mix(i) = 0.
501        dcfc_sub(i) = 0.
502        dqic_sub(i) = 0.
503        dqtc_sub(i) = 0.
504        dcfc_mix(i) = 0.
505        dqic_mix(i) = 0.
506        dqtc_mix(i) = 0.
507      ELSE
508        lincontfra(i) = 0.
509        circontfra(i) = 0.
510        qlincont(i) = 0.
511        qcircont(i) = 0.
512      ENDIF
513
514
515      !----------------------------------------------------------------------
516      !--    SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CONTRAIL    --
517      !----------------------------------------------------------------------
518
519      !--If there is a linear contrail
520      IF ( lincontfra(i) .GT. eps ) THEN
521
522        !--The contrail is always adjusted to saturation
523        qiceincld = ( qlincont(i) / lincontfra(i) - qsat(i) )
524        !--If the ice water content is too low, the cloud is purely sublimated
525        IF ( qiceincld .LT. qiceincld_min ) THEN
526          dcfl_sub(i) = - lincontfra(i)
527          dqil_sub(i) = - qiceincld * lincontfra(i)
528          dqtl_sub(i) = - qlincont(i)
529          lincontfra(i) = 0.
530          qlincont(i) = 0.
531          clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcfl_sub(i))
532          qclr(i) = qclr(i) - dqtl_sub(i)
533        ENDIF   ! qiceincld .LT. eps
534
535        !--We remove contrails from the main class
536        cldfra(i) = MAX(0., cldfra(i) - lincontfra(i))
537        qcld(i) = MAX(0., qcld(i) - qlincont(i))
538        qvc(i) = MAX(0., qvc(i) - qsat(i) * lincontfra(i))
539      ENDIF ! lincontfra(i) .GT. eps
540
541      !--If there is a contrail cirrus
542      IF ( circontfra(i) .GT. eps ) THEN
543
544        !--The contrail is always adjusted to saturation
545        qiceincld = ( qcircont(i) / circontfra(i) - qsat(i) )
546        !--If the ice water content is too low, the cloud is purely sublimated
547        IF ( qiceincld .LT. qiceincld_min ) THEN
548          dcfc_sub(i) = - circontfra(i)
549          dqic_sub(i) = - qiceincld * circontfra(i)
550          dqtc_sub(i) = - qcircont(i)
551          circontfra(i) = 0.
552          qcircont(i) = 0.
553          clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcfc_sub(i))
554          qclr(i) = qclr(i) - dqtc_sub(i)
555        ENDIF   ! qiceincld .LT. eps
556
557        !--We remove contrails from the main class
558        cldfra(i) = MAX(0., cldfra(i) - circontfra(i))
559        qcld(i) = MAX(0., qcld(i) - qcircont(i))
560        qvc(i) = MAX(0., qvc(i) - qsat(i) * circontfra(i))
561      ENDIF ! circontfra(i) .GT. eps
562
563
564      !-------------------------------------------------------------------
565      !--    SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD    --
566      !-------------------------------------------------------------------
567     
568      !--If there is a cloud
569      IF ( cldfra(i) .GT. eps ) THEN
570
571        qvapincld = qvc(i) / cldfra(i)
572        IF ( qvapincld .GT. gamma_cond(i) * qsat(i) ) THEN
573          qvapincld = gamma_cond(i) * qsat(i)
574          qvc(i) = qvapincld * cldfra(i)
575        ENDIF
576        qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
577       
578        !--If the ice water content is too low, the cloud is purely sublimated
579        !--Most probably, we advected a cloud with no ice water content (possible
580        !--if the entire cloud precipited for example)
581        IF ( qiceincld .LT. eps ) THEN
582          dcf_sub(i) = - cldfra(i)
583          dqvc_sub(i) = - qvc(i)
584          dqi_sub(i) = - ( qcld(i) - qvc(i) )
585
586          cldfra(i) = 0.
587          qcld(i) = 0.
588          qvc(i) = 0.
589          clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcf_sub(i))
590          qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i)
591
592        !--Else, the cloud is adjusted and sublimated
593        ELSE
594
595          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
596            IF ( qvapincld .GE. qsat(i) ) THEN
597              !--If the cloud is initially supersaturated
598              !--Exact explicit integration (qvc exact, qice explicit)
599              tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub
600              qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub )
601            ELSE
602              !--If the cloud is initially subsaturated
603              !!--Exact explicit integration (qvc exact, qice explicit)
604              !!--Same but depo_coef_cirrus = 1
605              !tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
606              !qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub )
607              !--Exact explicit integration (qice exact, qvc explicit)
608              !--The barrier is set so that the resulting vapor in cloud
609              !--cannot be greater than qsat
610              !--qice_ratio is the ratio between the new ice content and
611              !--the old one, it is comprised between 0 and 1
612              tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
613              qice_ratio = tauinv_depsub * dtime / 1.5 / qiceincld * ( qsat(i) - qvapincld )
614              !--The new vapor in the cloud is increased with the
615              !--sublimated ice
616              qvapincld_new = qvapincld + qiceincld * ( 1. - MAX(0., 1. - qice_ratio)**1.5 )
617              !--The new vapor in the cloud cannot be greater than qsat
618              qvapincld_new = MIN(qvapincld_new, qsat(i))
619              !--If all the ice is sublimated
620              IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0.
621            ENDIF ! qvapincld .GT. qsat
622          ELSE
623            !--We keep the saturation adjustment hypothesis, and the vapor in the
624            !--cloud is set equal to the saturation vapor
625            IF ( ( qvapincld + qiceincld ) .GT. qsat(i) ) THEN
626              qvapincld_new = qsat(i)
627            ELSE
628              qvapincld_new = 0.
629            ENDIF
630          ENDIF ! ok_unadjusted_clouds
631         
632
633          !------------------------------------
634          !--    DISSIPATION OF THE CLOUD    --
635          !------------------------------------
636          !--Additionally to a minimum in cloud water vapor, we impose a minimum
637          !--on the in-cloud ice water content. It is calculated following
638          !--Marti and Mauersberger (1993), see also Schiller et al. (2008)
639          qiceincld_min = qsati_diss(i) - qsat(i)
640
641          !--If the dissipation process must be activated
642          IF ( ( qvapincld_new + qiceincld_min ) .GT. qvapincld ) THEN
643            !--Gamma distribution starting at qvapincld
644            pdf_shape = nu_iwc_pdf_lscp / qiceincld
645            pdf_y = pdf_shape * ( qvapincld_new + qiceincld_min - qvapincld )
646            pdf_e1 = GAMMAINC ( nu_iwc_pdf_lscp      , pdf_y )
647            pdf_e2 = GAMMAINC ( nu_iwc_pdf_lscp + 1. , pdf_y )
648
649            !--Tendencies and diagnostics
650            dcf_sub(i)  = - cldfra(i) * pdf_e1
651            dqi_sub(i)  = - cldfra(i) * pdf_e2 / pdf_shape
652            dqvc_sub(i) = dcf_sub(i) * qvapincld
653
654            !--Add tendencies
655            cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
656            qcld(i)   = qcld(i) + dqvc_sub(i) + dqi_sub(i)
657            qvc(i)    = qvc(i) + dqvc_sub(i)
658            clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcf_sub(i))
659            qclr(i)   = qclr(i) - dqvc_sub(i) - dqi_sub(i)
660          ELSEIF ( qvapincld_new .EQ. 0. ) THEN
661            !--If all the ice has been sublimated, we sublimate
662            !--completely the cloud and do not activate the dissipation
663            !--process
664            !--Tendencies and diagnostics
665            dcf_sub(i) = - cldfra(i)
666            dqvc_sub(i) = - qvc(i)
667            dqi_sub(i) = - ( qcld(i) - qvc(i) )
668
669            !--Add tendencies
670            cldfra(i) = 0.
671            qcld(i) = 0.
672            qvc(i) = 0.
673            clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcf_sub(i))
674            qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i)
675          ENDIF ! qvapincld_new .GT. qvapincld
676
677
678          !------------------------------------
679          !--       PHASE ADJUSTMENT        --
680          !------------------------------------
681
682          IF ( qvapincld_new .GT. 0. ) THEN
683            !--Adjustment of the IWC to the new vapor in cloud
684            !--(this can be either positive or negative)
685            dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
686            dqi_adj(i) = - dqvc_adj(i)
687
688            !--Add tendencies
689            !--The vapor in the cloud is updated, but not qcld as it is constant
690            !--through this process, as well as cldfra which is unmodified
691            qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
692          ENDIF
693
694        ENDIF   ! qiceincld .LT. eps
695      ENDIF ! cldfra(i) .GT. eps
696
697
698      !--------------------------------------------------------------------------
699      !--    CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS    --
700      !--------------------------------------------------------------------------
701      !--This section relies on a distribution of water in the clear-sky region of
702      !--the mesh.
703
704      !--If there is a clear-sky region
705      IF ( clrfra(i) .GT. eps ) THEN
706
707        !--New PDF
708        rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
709        rhl_clr = MAX(0., MIN(150., rhl_clr))
710
711        !--Calculation of the properties of the PDF
712        !--Parameterization from IAGOS observations
713        !--pdf_alpha, pdf_scale and pdf_gamma will be reused below
714
715        !--Coefficient for standard deviation:
716        !--  tuning coef * (clear sky area**0.25) * (function of temperature)
717        pdf_e1 = beta_pdf_lscp * ( clrfra(i) * cell_area(i) )**0.25 &
718                * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
719        IF ( rhl_clr .GT. 50. ) THEN
720          pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp
721        ELSE
722          pdf_std = pdf_e1 * rhl_clr / 50.
723        ENDIF
724        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
725        pdf_alpha(i) = EXP( rhl_clr / 100. ) * pdf_e3
726        pdf_alpha(i) = MIN(10., pdf_alpha(i)) !--Avoid overflows
727       
728        !IF ( ok_warm_cloud ) THEN
729        !  !--If the statistical scheme is activated, the standard deviation is adapted
730        !  !--to depend on the pressure level. It is multiplied by ratqs, so that near the
731        !  !--surface std is almost 0, and upper than about 450 hPa the std is left untouched
732        !  pdf_std = pdf_std * ratqs(i)
733        !ENDIF
734
735        pdf_gamma(i) = GAMMA(1. + 1. / pdf_alpha(i))
736        !--Barrier to avoid overflows
737        pdf_scale(i) = MAX(eps, MIN(rhl_clr / pdf_gamma(i), pdf_std / SQRT( &
738                                    GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 )))
739        pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
740
741        !--Calculation of the newly condensed water and fraction (pronostic)
742        !--Integration of the clear sky PDF between gamma_cond*qsat and +inf
743        !--NB. the calculated values are clear-sky averaged
744
745        pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
746        pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
747        IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
748          pdf_fra_above_nuc = 0.
749          pdf_q_above_nuc = 0.
750        ELSEIF ( pdf_y .LT. -10. ) THEN
751          pdf_fra_above_nuc = 1.
752          pdf_q_above_nuc = qclr(i) / clrfra(i)
753        ELSE
754          pdf_y = EXP( pdf_y )
755          pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
756          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
757          pdf_fra_above_nuc = EXP( - pdf_y )
758          pdf_q_above_nuc = ( pdf_e3 + pdf_loc * pdf_fra_above_nuc ) * qsatl(i) / 100.
759        ENDIF
760
761        IF ( pdf_fra_above_nuc .GT. eps ) THEN
762
763          dcf_con(i) = clrfra(i) * pdf_fra_above_nuc
764          dqt_con = clrfra(i) * pdf_q_above_nuc
765
766          !--Barriers (should be useless
767          dcf_con(i) = MIN(dcf_con(i), clrfra(i))
768          dqt_con = MIN(dqt_con, qclr(i))
769
770          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
771            !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
772            !--the new vapor qvapincld. The timestep is divided by two because we do not
773            !--know when the condensation occurs
774            qvapincld = gamma_cond(i) * qsat(i)
775            qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
776            tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub
777            qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) &
778                          * EXP( - dtime / 2. * tauinv_depsub )
779          ELSE
780            !--We keep the saturation adjustment hypothesis, and the vapor in the
781            !--newly formed cloud is set equal to the saturation vapor.
782            qvapincld_new = qsat(i)
783          ENDIF
784
785          !--Tendency on cloud vapor and diagnostic
786          dqvc_con(i) = qvapincld_new * dcf_con(i)
787          dqi_con(i) = dqt_con - dqvc_con(i)
788
789          !--Add tendencies
790          cldfra(i)  = cldfra(i) + dcf_con(i)
791          qcld(i)    = qcld(i) + dqt_con
792          qvc(i)     = qvc(i) + dqvc_con(i)
793          clrfra(i)  = clrfra(i) - dcf_con(i)
794          qclr(i)    = qclr(i) - dqt_con
795
796        ENDIF ! pdf_fra_above_nuc .GT. eps
797      ELSE
798        !--Default value for the clear sky distribution: homogeneous distribution
799        pdf_alpha(i) = 1.
800        pdf_gamma(i) = 1.
801        pdf_scale(i) = eps
802      ENDIF ! clrfra(i) .GT. eps
803
804
805      !--------------------------------------
806      !--           CLOUD MIXING           --
807      !--------------------------------------
808      !--This process mixes the cloud with its surroundings: the subsaturated clear sky,
809      !--and the supersaturated clear sky. It is activated if the cloud is big enough,
810      !--but does not cover the entire mesh.
811      !
812      IF ( ( cldfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN
813
814        !-- PART 1 - TURBULENT DIFFUSION
815
816        !--Clouds within the mesh are assumed to be ellipses. The length of the
817        !--semi-major axis is a and the length of the semi-minor axis is b.
818        !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer.
819        !--In particular, it is 0 if cldfra = 1.
820        !--clouds_perim is the total perimeter of the clouds within the mesh,
821        !--not considering interfaces with other meshes (only the interfaces with clear
822        !--sky are taken into account).
823        !--
824        !--The area of each cloud is A = a * b * RPI,
825        !--and the perimeter of each cloud is
826        !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
827        !--
828        !--With cell_area the area of the cell, we have:
829        !-- cldfra = A * N_cld_mix / cell_area
830        !-- clouds_perim = P * N_cld_mix
831        !--
832        !--We assume that the ratio between b and a is a function of
833        !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
834        !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
835        !--are spherical.
836        !-- b / a = bovera = MAX(0.1, cldfra)
837        !bovera = MAX(0.1, cldfra(i))
838        bovera = aspect_ratio_cirrus
839        !--P / a is a function of b / a only, that we can calculate
840        !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
841        Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
842        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
843        !--based on observations.
844        !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
845        !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
846        !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
847        !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
848        a_mix = Povera / coef_mixing_lscp / RPI / ( 1. - cldfra(i) ) / bovera
849        !--and finally,
850        !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
851        N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) ) * cell_area(i) &
852                  / Povera / a_mix
853
854        !--The time required for turbulent diffusion to homogenize a region of size
855        !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
856        !--We compute L_mix and assume that the cloud is mixed over this length
857        L_mix = SQRT( dtime**3 * pbl_eps(i) )
858        !--The mixing length cannot be greater than the semi-minor axis. In this case,
859        !--the entire cloud is mixed.
860        L_mix = MIN(L_mix, a_mix * bovera)
861
862        !--The fraction of clear sky mixed is
863        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
864        clrfra_mix = N_cld_mix * RPI / cell_area(i) &
865                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
866        !--The fraction of clear sky mixed is
867        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
868        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
869                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
870
871
872        !-- PART 2 - SHEARING
873
874        !--The clouds are then sheared. We keep the shape and number
875        !--assumptions from before. The clouds are sheared along their
876        !--semi-major axis (a_mix), on the entire cell heigh dz.
877        !--The increase in size is
878        L_shear = coef_shear_lscp * shear(i) * dz * dtime
879        !--therefore, the fraction of clear sky mixed is
880        !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
881        !--and the fraction of cloud mixed is
882        !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area
883        !--(note that they are equal)
884        shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i)
885        !--and the environment and cloud mixed fractions are the same,
886        !--which we add to the previous calculated mixed fractions.
887        !--We therefore assume that the sheared clouds and the turbulent
888        !--mixed clouds are different.
889        clrfra_mix = clrfra_mix + shear_fra
890        cldfra_mix = cldfra_mix + shear_fra
891
892        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
893
894        clrfra_mix = MIN(clrfra(i), clrfra_mix)
895        cldfra_mix = MIN(cldfra(i), cldfra_mix)
896
897        !--We compute the limit vapor in clear sky where the mixed cloud could not
898        !--survive if all the ice crystals were sublimated. Note that here we assume,
899        !--for growth or reduction of the cloud, that the mixed cloud is adjusted
900        !--to saturation, ie the vapor in the mixed cloud is qsat. This is only a
901        !--diagnostic, and if the cloud size is increased, we add the new vapor to the
902        !--cloud's vapor without condensing or sublimating ice crystals
903        IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
904          qiceinmix = ( qcld(i) - qvc(i) ) / cldfra(i) / ( 1. + clrfra_mix / cldfra_mix )
905          tauinv_depsub = qiceinmix**(1./3.) * kappa_depsub
906          !qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime * tauinv_depsub ) )
907          qvapinmix_lim = qsat(i) - qiceinmix * MAX(1., 1.5 / ( dtime * tauinv_depsub ))
908          qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) &
909                        - qvc(i) / cldfra(i) * cldfra_mix / clrfra_mix
910        ELSE
911          !--NB. if tau_depsub = 0 (ie tauinv_depsub = inf), we get the same result as above
912          qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) &
913                        - qcld(i) / cldfra(i) * cldfra_mix / clrfra_mix
914        ENDIF
915
916        IF ( qvapinclr_lim .LT. 0. ) THEN
917          !--Whatever we do, the cloud will increase in size
918          dcf_mix(i)  = clrfra_mix
919          dqvc_mix(i) = clrfra_mix * qclr(i) / clrfra(i)
920        ELSE
921          !--We then calculate the clear sky part where the humidity is lower than
922          !--qvapinclr_lim, and the part where it is higher than qvapinclr_lim
923          !--This is the clear-sky PDF calculated in the condensation section. Note
924          !--that if we are here, we necessarily went through the condensation part
925          !--because the clear sky fraction can only be reduced by condensation.
926          !--Thus the `pdf_xxx` variables are well defined.
927
928          rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
929          pdf_x = qvapinclr_lim / qsatl(i) * 100.
930          pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
931          pdf_x = qsat(i) / qsatl(i) * 100.
932          pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
933          IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
934            pdf_fra_above_lim = 0.
935            pdf_q_above_lim = 0.
936          ELSEIF ( pdf_y .LT. -10. ) THEN
937            pdf_fra_above_lim = clrfra(i)
938            pdf_q_above_lim = qclr(i)
939          ELSE
940            pdf_y = EXP( pdf_y )
941            pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
942            pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
943            pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
944            pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
945                              + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
946          ENDIF
947
948          pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim
949
950          !--sigma_mix is the ratio of the surroundings of the clouds where mixing
951          !--increases the size of the cloud, to the total surroundings of the clouds.
952          !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing
953          !--decreases the size of the clouds
954          sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim )
955
956          IF ( pdf_fra_above_lim .GT. eps ) THEN
957            dcf_mix(i)  = clrfra_mix * sigma_mix
958            dqvc_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim
959          ENDIF
960
961          IF ( pdf_fra_below_lim .GT. eps ) THEN
962            dcf_mix(i)  = dcf_mix(i)  - cldfra_mix * ( 1. - sigma_mix )
963            dqvc_mix(i) = dqvc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
964                                      * qvc(i) / cldfra(i)
965            dqi_mix(i)  = dqi_mix(i)  - cldfra_mix * ( 1. - sigma_mix ) &
966                                      * ( qcld(i) - qvc(i) ) / cldfra(i)
967          ENDIF
968
969        ENDIF
970      ENDIF ! ( cldfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps )
971
972      !--------------------------------------
973      !--        CONTRAIL MIXING           --
974      !--------------------------------------
975
976      IF ( ( lincontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps ) ) THEN
977
978        !-- PART 1 - TURBULENT DIFFUSION
979       
980        !--Clouds within the mesh are assumed to be ellipses. The length of the
981        !--semi-major axis is a and the length of the semi-minor axis is b.
982        !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer.
983        !--In particular, it is 0 if cldfra = 1.
984        !--clouds_perim is the total perimeter of the clouds within the mesh,
985        !--not considering interfaces with other meshes (only the interfaces with clear
986        !--sky are taken into account).
987        !--
988        !--The area of each cloud is A = a * b * RPI,
989        !--and the perimeter of each cloud is
990        !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
991        !--
992        !--With cell_area the area of the cell, we have:
993        !-- cldfra = A * N_cld_mix / cell_area
994        !-- clouds_perim = P * N_cld_mix
995        !--
996        !--We assume that the ratio between b and a is a function of
997        !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
998        !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
999        !--are spherical.
1000        !-- b / a = bovera = MAX(0.1, cldfra)
1001        bovera = aspect_ratio_lincontrails
1002        !--P / a is a function of b / a only, that we can calculate
1003        !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
1004        Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
1005     
1006        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
1007        !--based on observations.
1008        !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
1009        !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
1010        !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
1011        !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
1012        a_mix = Povera / coef_mixing_lincontrails / RPI / ( 1. - lincontfra(i) ) / bovera
1013        !--and finally,
1014        !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
1015        N_cld_mix = coef_mixing_lincontrails * lincontfra(i) * ( 1. - lincontfra(i) ) &
1016                  * cell_area(i) / Povera / a_mix
1017       
1018        !--The time required for turbulent diffusion to homogenize a region of size
1019        !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
1020        !--We compute L_mix and assume that the cloud is mixed over this length
1021        L_mix = SQRT( dtime**3 * pbl_eps(i) )
1022        !--The mixing length cannot be greater than the semi-minor axis. In this case,
1023        !--the entire cloud is mixed.
1024        L_mix = MIN(L_mix, a_mix * bovera)
1025       
1026        !--The fraction of clear sky mixed is
1027        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
1028        clrfra_mix = N_cld_mix * RPI / cell_area(i) &
1029                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
1030        !--The fraction of clear sky mixed is
1031        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
1032        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
1033                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
1034       
1035       
1036        !-- PART 2 - SHEARING
1037       
1038        !--The clouds are then sheared. We keep the shape and number
1039        !--assumptions from before. The clouds are sheared with a random orientation
1040        !--of the wind, on average we assume that the wind and the semi-major axis
1041        !--make a 45 degrees angle. Moreover, the contrails only mix
1042        !--along their semi-minor axis (b), because it is easier to compute.
1043        !--With this, the clouds increase in size along b only, by a factor
1044        !--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind)
1045        L_shear = coef_shear_lincontrails * shear(i) * dz * dtime
1046        !--therefore, the fraction of clear sky mixed is
1047        !-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area
1048        !--and the fraction of cloud mixed is
1049        !-- N_cld_mix * ( a * b - a * (b - L_shear * SQRT(2.) / 2.) ) * RPI / 2. / cell_area
1050        !--(note that they are equal)
1051        shear_fra = RPI * L_shear * a_mix * SQRT(2.) / 2. / 2. * N_cld_mix / cell_area(i)
1052        !--and the environment and cloud mixed fractions are the same,
1053        !--which we add to the previous calculated mixed fractions.
1054        !--We therefore assume that the sheared clouds and the turbulent
1055        !--mixed clouds are different.
1056        clrfra_mix = clrfra_mix + shear_fra
1057        cldfra_mix = cldfra_mix + shear_fra
1058       
1059       
1060        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
1061       
1062        clrfra_mix = MIN(clrfra(i), clrfra_mix)
1063        cldfra_mix = MIN(lincontfra(i), cldfra_mix)
1064       
1065        !--We compute the limit vapor in clear sky where the mixed cloud could not
1066        !--survive if all the ice crystals were sublimated. Note that here we assume,
1067        !--for growth or reduction of the cloud, that the mixed cloud is adjusted
1068        !--to saturation, ie the vapor in the mixed cloud is qsat. This is only a
1069        !--diagnostic, and if the cloud size is increased, we add the new vapor to the
1070        !--cloud's vapor without condensing or sublimating ice crystals
1071        qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) &
1072                      - qlincont(i) / lincontfra(i) * cldfra_mix / clrfra_mix
1073
1074        IF ( qvapinclr_lim .LT. 0. ) THEN
1075          !--Whatever we do, the cloud will increase in size
1076          !--If the linear contrail increases in size, the increment is considered
1077          !--to be a contrail cirrus
1078          dcfc_mix(i) = dcfc_mix(i) + clrfra_mix
1079          dqtc_mix(i) = dqtc_mix(i) + clrfra_mix * qclr(i) / clrfra(i)
1080          dqic_mix(i) = dqic_mix(i) + clrfra_mix * sigma_mix &
1081                      * ( qclr(i) / clrfra(i) - qsat(i) )
1082        ELSE
1083          !--We then calculate the clear sky part where the humidity is lower than
1084          !--qvapinclr_lim, and the part where it is higher than qvapinclr_lim
1085          !--This is the clear-sky PDF calculated in the condensation section. Note
1086          !--that if we are here, we necessarily went through the condensation part
1087          !--because the clear sky fraction can only be reduced by condensation.
1088          !--Thus the `pdf_xxx` variables are well defined.
1089
1090          rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
1091          pdf_x = qvapinclr_lim / qsatl(i) * 100.
1092          pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
1093          pdf_x = qsat(i) / qsatl(i) * 100.
1094          pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
1095          IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
1096            pdf_fra_above_lim = 0.
1097            pdf_q_above_lim = 0.
1098          ELSEIF ( pdf_y .LT. -10. ) THEN
1099            pdf_fra_above_lim = clrfra(i)
1100            pdf_q_above_lim = qclr(i)
1101          ELSE
1102            pdf_y = EXP( pdf_y )
1103            pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
1104            pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
1105            pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
1106            pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
1107                              + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
1108          ENDIF
1109
1110          pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim
1111
1112          !--sigma_mix is the ratio of the surroundings of the clouds where mixing
1113          !--increases the size of the cloud, to the total surroundings of the clouds.
1114          !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing
1115          !--decreases the size of the clouds
1116          !--For aviation, we increase the chance that the air surrounding contrails
1117          !--is moist. This is quantified with chi_mixing_lincontrails
1118          sigma_mix = chi_mixing_lincontrails * pdf_fra_above_lim &
1119                    / ( pdf_fra_below_lim + chi_mixing_lincontrails * pdf_fra_above_lim )
1120
1121          IF ( pdf_fra_above_lim .GT. eps ) THEN
1122            !--If the linear contrail increases in size, the increment is considered
1123            !--to be a contrail cirrus
1124            qvapinmix = ( pdf_q_above_lim / pdf_fra_above_lim * clrfra_mix &
1125                      + qlincont(i) / lincontfra(i) * cldfra_mix ) &
1126                      / ( clrfra_mix + cldfra_mix )
1127            qiceinmix = ( qlincont(i) / lincontfra(i) - qsat(i) ) * cldfra_mix &
1128                      / ( clrfra_mix + cldfra_mix )
1129            dcfc_mix(i) = dcfc_mix(i) + clrfra_mix * sigma_mix
1130            dqtc_mix(i) = dqtc_mix(i) + clrfra_mix * sigma_mix * qvapinmix
1131            dqtl_mix(i) = dqtl_mix(i) - cldfra_mix * sigma_mix &
1132                        * ( qlincont(i) / lincontfra(i) - qvapinmix )
1133            dqic_mix(i) = dqic_mix(i) + clrfra_mix * sigma_mix * qiceinmix
1134            dqil_mix(i) = dqil_mix(i) - cldfra_mix * sigma_mix &
1135                        * ( qlincont(i) / lincontfra(i) - qsat(i) - qiceinmix )
1136          ENDIF
1137
1138          IF ( pdf_fra_below_lim .GT. eps ) THEN
1139            dcfl_mix(i) = dcfl_mix(i) - cldfra_mix * ( 1. - sigma_mix )
1140            dqtl_mix(i) = dqtl_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
1141                        * qlincont(i) / lincontfra(i)
1142            dqil_mix(i) = dqil_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
1143                        * ( qlincont(i) / lincontfra(i) - qsat(i) )
1144          ENDIF
1145
1146        ENDIF
1147      ENDIF ! ( lincontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps )
1148
1149      IF ( ( circontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps ) ) THEN
1150
1151        !-- PART 1 - TURBULENT DIFFUSION
1152       
1153        !--Clouds within the mesh are assumed to be ellipses. The length of the
1154        !--semi-major axis is a and the length of the semi-minor axis is b.
1155        !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer.
1156        !--In particular, it is 0 if cldfra = 1.
1157        !--clouds_perim is the total perimeter of the clouds within the mesh,
1158        !--not considering interfaces with other meshes (only the interfaces with clear
1159        !--sky are taken into account).
1160        !--
1161        !--The area of each cloud is A = a * b * RPI,
1162        !--and the perimeter of each cloud is
1163        !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
1164        !--
1165        !--With cell_area the area of the cell, we have:
1166        !-- cldfra = A * N_cld_mix / cell_area
1167        !-- clouds_perim = P * N_cld_mix
1168        !--
1169        !--We assume that the ratio between b and a is a function of
1170        !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
1171        !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
1172        !--are spherical.
1173        !-- b / a = bovera = MAX(0.1, cldfra)
1174        bovera = aspect_ratio_cirrus
1175        !--P / a is a function of b / a only, that we can calculate
1176        !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
1177        Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
1178     
1179        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
1180        !--based on observations.
1181        !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
1182        !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
1183        !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
1184        !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
1185        a_mix = Povera / coef_mixing_lscp / RPI / ( 1. - circontfra(i) ) / bovera
1186        !--and finally,
1187        !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
1188        N_cld_mix = coef_mixing_lscp * circontfra(i) * ( 1. - circontfra(i) ) &
1189                  * cell_area(i) / Povera / a_mix
1190       
1191        !--The time required for turbulent diffusion to homogenize a region of size
1192        !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
1193        !--We compute L_mix and assume that the cloud is mixed over this length
1194        L_mix = SQRT( dtime**3 * pbl_eps(i) )
1195        !--The mixing length cannot be greater than the semi-minor axis. In this case,
1196        !--the entire cloud is mixed.
1197        L_mix = MIN(L_mix, a_mix * bovera)
1198       
1199        !--The fraction of clear sky mixed is
1200        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
1201        clrfra_mix = N_cld_mix * RPI / cell_area(i) &
1202                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
1203        !--The fraction of clear sky mixed is
1204        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
1205        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
1206                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
1207       
1208       
1209        !-- PART 2 - SHEARING
1210       
1211        !--The clouds are then sheared. We keep the shape and number
1212        !--assumptions from before. The clouds are sheared along their
1213        !--semi-major axis (a_mix), on the entire cell heigh dz.
1214        !--The increase in size is
1215        L_shear = coef_shear_lscp * shear(i) * dz * dtime
1216        !--therefore, the fraction of clear sky mixed is
1217        !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
1218        !--and the fraction of cloud mixed is
1219        !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area
1220        !--(note that they are equal)
1221        shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i)
1222        !--and the environment and cloud mixed fractions are the same,
1223        !--which we add to the previous calculated mixed fractions.
1224        !--We therefore assume that the sheared clouds and the turbulent
1225        !--mixed clouds are different.
1226        clrfra_mix = clrfra_mix + shear_fra
1227        cldfra_mix = cldfra_mix + shear_fra
1228       
1229       
1230        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
1231       
1232        clrfra_mix = MIN(clrfra(i), clrfra_mix)
1233        cldfra_mix = MIN(circontfra(i), cldfra_mix)
1234       
1235        !--We compute the limit vapor in clear sky where the mixed cloud could not
1236        !--survive if all the ice crystals were sublimated. Note that here we assume,
1237        !--for growth or reduction of the cloud, that the mixed cloud is adjusted
1238        !--to saturation, ie the vapor in the mixed cloud is qsat. This is only a
1239        !--diagnostic, and if the cloud size is increased, we add the new vapor to the
1240        !--cloud's vapor without condensing or sublimating ice crystals
1241        qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) &
1242                      - qcircont(i) / circontfra(i) * cldfra_mix / clrfra_mix
1243
1244        IF ( qvapinclr_lim .LT. 0. ) THEN
1245          !--Whatever we do, the cloud will increase in size
1246          dcfc_mix(i) = dcfc_mix(i) + clrfra_mix
1247          dqtc_mix(i) = dqtc_mix(i) + clrfra_mix * qclr(i) / clrfra(i)
1248          dqic_mix(i) = dqic_mix(i) + clrfra_mix * sigma_mix &
1249                      * ( qclr(i) / clrfra(i) - qsat(i) )
1250        ELSE
1251          !--We then calculate the clear sky part where the humidity is lower than
1252          !--qvapinclr_lim, and the part where it is higher than qvapinclr_lim
1253          !--This is the clear-sky PDF calculated in the condensation section. Note
1254          !--that if we are here, we necessarily went through the condensation part
1255          !--because the clear sky fraction can only be reduced by condensation.
1256          !--Thus the `pdf_xxx` variables are well defined.
1257
1258          rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
1259          pdf_x = qvapinclr_lim / qsatl(i) * 100.
1260          pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
1261          pdf_x = qsat(i) / qsatl(i) * 100.
1262          pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
1263          IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
1264            pdf_fra_above_lim = 0.
1265            pdf_q_above_lim = 0.
1266          ELSEIF ( pdf_y .LT. -10. ) THEN
1267            pdf_fra_above_lim = clrfra(i)
1268            pdf_q_above_lim = qclr(i)
1269          ELSE
1270            pdf_y = EXP( pdf_y )
1271            pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
1272            pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
1273            pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
1274            pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
1275                              + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
1276          ENDIF
1277
1278          pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim
1279
1280          !--sigma_mix is the ratio of the surroundings of the clouds where mixing
1281          !--increases the size of the cloud, to the total surroundings of the clouds.
1282          !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing
1283          !--decreases the size of the clouds
1284          sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim )
1285
1286          IF ( pdf_fra_above_lim .GT. eps ) THEN
1287            dcfc_mix(i) = dcfc_mix(i) + clrfra_mix * sigma_mix
1288            dqtc_mix(i) = dqtc_mix(i) + clrfra_mix * sigma_mix &
1289                        * pdf_q_above_lim / pdf_fra_above_lim
1290            dqic_mix(i) = dqic_mix(i) + clrfra_mix * sigma_mix &
1291                        * ( pdf_q_above_lim / pdf_fra_above_lim - qsat(i) )
1292          ENDIF
1293
1294          IF ( pdf_fra_below_lim .GT. eps ) THEN
1295            dcfc_mix(i) = dcfc_mix(i) - cldfra_mix * ( 1. - sigma_mix )
1296            dqtc_mix(i) = dqtc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
1297                        * qcircont(i) / circontfra(i)
1298            dqic_mix(i) = dqic_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
1299                        * ( qcircont(i) / circontfra(i) - qsat(i) )
1300          ENDIF
1301
1302        ENDIF
1303      ENDIF ! ( circontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps )
1304
1305      !--We balance the mixing tendencies between the different cloud classes
1306      mixed_fraction = dcf_mix(i) + dcfl_mix(i) + dcfc_mix(i)
1307      IF ( mixed_fraction .GT. clrfra(i) ) THEN
1308        mixed_fraction = clrfra(i) / mixed_fraction
1309        dcf_mix(i)  = dcf_mix(i)  * mixed_fraction
1310        dqvc_mix(i) = dqvc_mix(i) * mixed_fraction
1311        dqi_mix(i)  = dqi_mix(i)  * mixed_fraction
1312        dcfl_mix(i) = dcfl_mix(i) * mixed_fraction
1313        dqtl_mix(i) = dqtl_mix(i) * mixed_fraction
1314        dqil_mix(i) = dqil_mix(i) * mixed_fraction
1315        dcfc_mix(i) = dcfc_mix(i) * mixed_fraction
1316        dqtc_mix(i) = dqtc_mix(i) * mixed_fraction
1317        dqic_mix(i) = dqic_mix(i) * mixed_fraction
1318      ENDIF
1319
1320      IF ( lincontfra(i) .GT. eps ) THEN
1321        !--Add tendencies
1322        lincontfra(i) = lincontfra(i) + dcfl_mix(i)
1323        qlincont(i)   = qlincont(i) + dqtl_mix(i)
1324        clrfra(i)  = clrfra(i) - dcfl_mix(i)
1325        qclr(i)    = qclr(i) - dqtl_mix(i)
1326      ENDIF
1327
1328      IF ( circontfra(i) .GT. eps ) THEN
1329        !--Add tendencies
1330        circontfra(i) = circontfra(i) + dcfc_mix(i)
1331        qcircont(i)   = qcircont(i) + dqtc_mix(i)
1332        clrfra(i)  = clrfra(i) - dcfc_mix(i)
1333        qclr(i)    = qclr(i) - dqtc_mix(i)
1334      ENDIF
1335
1336      !--Add tendencies
1337      cldfra(i)  = cldfra(i) + dcf_mix(i)
1338      qcld(i)    = qcld(i) + dqvc_mix(i) + dqi_mix(i)
1339      qvc(i)     = qvc(i) + dqvc_mix(i)
1340      clrfra(i)  = clrfra(i) - dcf_mix(i)
1341      qclr(i)    = qclr(i) - dqvc_mix(i) - dqi_mix(i)
1342
1343
1344      !---------------------------------------
1345      !--         ICE SEDIMENTATION         --
1346      !---------------------------------------
1347      !
1348      !--If ice supersaturation is activated, the cloud properties are prognostic.
1349      !--The falling ice is then partly considered a new cloud in the gridbox.
1350      !--BEWARE with this parameterisation, we can create a new cloud with a much
1351      !--different ice water content and water vapor content than the existing cloud
1352      !--(if it exists). This implies than unphysical fluxes of ice and vapor
1353      !--occur between the existing cloud and the newly formed cloud (from sedimentation).
1354      !--Note also that currently, the precipitation scheme assume a maximum
1355      !--random overlap, meaning that the new formed clouds will not be affected
1356      !--by vertical wind shear.
1357      !
1358      IF ( icesed_flux(i) .GT. 0. ) THEN
1359
1360        clrfra_sed = MIN(clrfra(i), cldfra_above(i) - cldfra(i))
1361
1362        IF ( ( clrfra_sed .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN
1363
1364          qiceincld = qice_sedim(i) / cldfra_above(i)
1365          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
1366            tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
1367            qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) )
1368          ELSE
1369            qvapinclr_lim = qsat(i) - qiceincld
1370          ENDIF
1371
1372          rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
1373          pdf_x = qvapinclr_lim / qsatl(i) * 100.
1374          pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
1375          pdf_x = qsat(i) / qsatl(i) * 100.
1376          pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
1377          IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
1378            pdf_fra_above_lim = 0.
1379            pdf_q_above_lim = 0.
1380          ELSEIF ( pdf_y .LT. -10. ) THEN
1381            pdf_fra_above_lim = clrfra(i)
1382            pdf_q_above_lim = qclr(i)
1383          ELSE
1384            pdf_y = EXP( pdf_y )
1385            pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
1386            pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
1387            pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
1388            pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
1389                              + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
1390          ENDIF
1391
1392          IF ( pdf_fra_above_lim .GT. eps ) THEN
1393            dcf_sed(i)  = clrfra_sed * pdf_fra_above_lim / clrfra(i)
1394            dqvc_sed(i) = dcf_sed(i) * pdf_q_above_lim / pdf_fra_above_lim
1395          ENDIF
1396          !--We include the sedimentated ice:
1397          dqi_sed(i)  = qiceincld &    !--We include the sedimentated ice:
1398                      * ( dcf_sed(i) & !--the part that falls in the newly formed cloud and
1399                        + cldfra(i) )  !--the part that falls in the existing cloud
1400
1401        ELSE
1402
1403          dqi_sed(i) = qice_sedim(i)
1404
1405        ENDIF ! ( clrfra_sed .GT. eps .AND. ( clrfra(i) .GT. eps )
1406
1407        !--Add tendencies
1408        cldfra(i)  = MIN(totfra_in(i), cldfra(i) + dcf_sed(i))
1409        qcld(i)    = qcld(i) + dqvc_sed(i) + dqi_sed(i)
1410        qvc(i)     = qvc(i) + dqvc_sed(i)
1411        clrfra(i)  = MAX(0., clrfra(i) - dcf_sed(i))
1412        !--We re-include sublimated sedimentated ice into clear sky water vapor
1413        qclr(i)    = qclr(i) - dqvc_sed(i) + qice_sedim(i) - dqi_sed(i)
1414
1415      ENDIF ! icesed_flux(i) .GT. 0.
1416
1417
1418      !--We put back contrails in the clouds class
1419      IF ( ( lincontfra(i) + circontfra(i) ) .GT. 0. ) THEN
1420        cldfra(i) = cldfra(i) + lincontfra(i) + circontfra(i)
1421        qcld(i) = qcld(i) + qlincont(i) + qcircont(i)
1422        qvc(i) = qvc(i) + qsat(i) * ( lincontfra(i) + circontfra(i) )
1423      ENDIF
1424
1425
1426      !--Diagnose ISSRs
1427      IF ( clrfra(i) .GT. eps ) THEN
1428        rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
1429        pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
1430        pdf_x = qsat(i) / qsatl(i) * 100.
1431        pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
1432        IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
1433          issrfra(i) = 0.
1434          qissr(i) = 0.
1435        ELSEIF ( pdf_y .LT. -10. ) THEN
1436          issrfra(i) = clrfra(i)
1437          qissr(i) = qclr(i)
1438        ELSE
1439          pdf_y = EXP( pdf_y )
1440          pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
1441          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
1442          issrfra(i) = EXP( - pdf_y ) * clrfra(i)
1443          qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
1444        ENDIF
1445        IF ( issrfra(i) .LT. eps ) THEN
1446          issrfra(i) = 0.
1447          qissr(i) = 0.
1448        ENDIF
1449      ELSE
1450        issrfra(i) = 0.
1451        qissr(i) = 0.
1452      ENDIF
1453
1454      !-------------------------------------------
1455      !--       FINAL BARRIERS AND OUTPUTS      --
1456      !-------------------------------------------
1457
1458      cldfra(i) = MIN(cldfra(i), totfra_in(i))
1459      qcld(i) = MIN(qcld(i), qtot_in(i))
1460      qvc(i) = MIN(qvc(i), qcld(i))
1461
1462      IF ( cldfra(i) .LT. eps ) THEN
1463        !--If the cloud is too small, it is sublimated.
1464        cldfra(i) = 0.
1465        qcld(i)   = 0.
1466        qvc(i)    = 0.
1467        qincld(i) = qsat(i)
1468      ELSE
1469        qincld(i) = qcld(i) / cldfra(i)
1470      ENDIF ! cldfra .LT. eps
1471
1472      !--Diagnostics
1473      dcf_sub(i)  = dcf_sub(i)  / dtime
1474      dcf_con(i)  = dcf_con(i)  / dtime
1475      dcf_mix(i)  = dcf_mix(i)  / dtime
1476      dcf_sed(i)  = dcf_sed(i)  / dtime
1477      dqi_adj(i)  = dqi_adj(i)  / dtime
1478      dqi_sub(i)  = dqi_sub(i)  / dtime
1479      dqi_con(i)  = dqi_con(i)  / dtime
1480      dqi_mix(i)  = dqi_mix(i)  / dtime
1481      dqi_sed(i)  = dqi_sed(i)  / dtime
1482      dqvc_adj(i) = dqvc_adj(i) / dtime
1483      dqvc_sub(i) = dqvc_sub(i) / dtime
1484      dqvc_con(i) = dqvc_con(i) / dtime
1485      dqvc_mix(i) = dqvc_mix(i) / dtime
1486      dqvc_sed(i) = dqvc_sed(i) / dtime
1487
1488    ENDIF ! pt_pron_clds(i)
1489
1490  ENDIF   ! end keepgoing
1491
1492ENDDO     ! end loop on i
1493
1494
1495!----------------------------------------
1496!--       CONTRAILS AND AVIATION       --
1497!----------------------------------------
1498IF ( ok_plane_contrail ) THEN
1499
1500  CALL contrails_formation( &
1501      klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, &
1502      flight_dist, clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, &
1503      keepgoing, pt_pron_clds, &
1504      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
1505      dcfl_ini, dqil_ini, dqtl_ini)
1506
1507  DO i = 1, klon
1508    IF ( keepgoing(i) .AND. pt_pron_clds(i) ) THEN
1509
1510      !--Convert existing contrail fraction into "natural" cirrus cloud fraction
1511      IF ( ( cldfra(i) .GE. ( totfra_in(i) - eps ) ) .OR. ( lincontfra(i) .LE. eps ) ) THEN
1512        contrails_conversion_factor = 1.
1513      ELSE
1514        contrails_conversion_factor = ( 1. - &
1515            !--First, the linear contrails are continuously degraded in induced cirrus
1516            EXP( - dtime / linear_contrails_lifetime ) &
1517            !--Second, if the cloud fills the entire gridbox, the linear contrails
1518            !--cannot exist. The exponent is set so that this only happens for
1519            !--very cloudy gridboxes
1520            * ( 1. - cldfra(i) / totfra_in(i) )**0.1 )
1521      ENDIF
1522      dcfl_cir(i) = - contrails_conversion_factor * lincontfra(i)
1523      dqtl_cir(i) = - contrails_conversion_factor * qlincont(i)
1524
1525      dcfl_ini(i) = MIN(MIN(dcfl_ini(i), issrfra(i)), totfra_in(i) - cldfra(i))
1526      dqtl_ini(i) = MIN(MIN(dqtl_ini(i), qissr(i)), qtot_in(i) - qcld(i))
1527
1528      !--Add tendencies
1529      issrfra(i) = issrfra(i) - dcfl_ini(i)
1530      qissr(i)   = qissr(i) - dqtl_ini(i)
1531      clrfra(i)  = clrfra(i) - dcfl_ini(i)
1532      qclr(i)    = qclr(i) - dqtl_ini(i)
1533
1534      cldfra(i)  = cldfra(i) + dcfl_ini(i)
1535      qcld(i)    = qcld(i) + dqtl_ini(i)
1536      qvc(i)     = qvc(i) + dcfl_ini(i) * qsat(i)
1537      lincontfra(i) = lincontfra(i) + dcfl_cir(i) + dcfl_ini(i)
1538      qlincont(i)   = qlincont(i) + dqtl_cir(i) + dqtl_ini(i)
1539      circontfra(i) = circontfra(i) - dcfl_cir(i)
1540      qcircont(i)   = qcircont(i) - dqtl_cir(i)
1541
1542      !--Diagnostics
1543      dcfl_ini(i) = dcfl_ini(i) / dtime
1544      dqil_ini(i) = dqil_ini(i) / dtime
1545      dqtl_ini(i) = dqtl_ini(i) / dtime
1546      dcfl_sub(i) = dcfl_sub(i) / dtime
1547      dqil_sub(i) = dqil_sub(i) / dtime
1548      dqtl_sub(i) = dqtl_sub(i) / dtime
1549      dcfl_cir(i) = dcfl_cir(i) / dtime
1550      dqtl_cir(i) = dqtl_cir(i) / dtime
1551      dcfl_mix(i) = dcfl_mix(i) / dtime
1552      dqil_mix(i) = dqil_mix(i) / dtime
1553      dqtl_mix(i) = dqtl_mix(i) / dtime
1554      dcfc_sub(i) = dcfc_sub(i) / dtime
1555      dqic_sub(i) = dqic_sub(i) / dtime
1556      dqtc_sub(i) = dqtc_sub(i) / dtime
1557      dcfc_mix(i) = dcfc_mix(i) / dtime
1558      dqic_mix(i) = dqic_mix(i) / dtime
1559      dqtc_mix(i) = dqtc_mix(i) / dtime
1560
1561      !-------------------------------------------
1562      !--       FINAL BARRIERS AND OUTPUTS      --
1563      !-------------------------------------------
1564
1565      cldfra(i) = MIN(cldfra(i), totfra_in(i))
1566      qcld(i) = MIN(qcld(i), qtot_in(i))
1567      qvc(i) = MIN(qvc(i), qcld(i))
1568
1569      IF ( cldfra(i) .LT. eps ) THEN
1570        !--If the cloud is too small, it is sublimated.
1571        cldfra(i) = 0.
1572        qcld(i)   = 0.
1573        qvc(i)    = 0.
1574        lincontfra(i) = 0.
1575        qlincont(i)   = 0.
1576        circontfra(i) = 0.
1577        qcircont(i)   = 0.
1578        qincld(i) = qsat(i)
1579      ELSE
1580        qincld(i) = qcld(i) / cldfra(i)
1581      ENDIF ! cldfra .LT. eps
1582
1583      IF ( (lincontfra(i) .LT. eps) .OR. (qlincont(i) .LT. (qsat(i) * lincontfra(i))) ) THEN
1584        lincontfra(i) = 0.
1585        qlincont(i) = 0.
1586      ENDIF
1587
1588      IF ( (circontfra(i) .LT. eps) .OR. (qcircont(i) .LT. (qsat(i) * circontfra(i))) ) THEN
1589        circontfra(i) = 0.
1590        qcircont(i) = 0.
1591      ENDIF
1592
1593    ENDIF ! keepgoing
1594  ENDDO ! loop on klon
1595ENDIF ! ok_plane_contrail
1596
1597
1598END SUBROUTINE condensation_ice_supersat
1599!**********************************************************************************
1600
1601END MODULE lmdz_lscp_condensation
Note: See TracBrowser for help on using the repository browser.