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

Last change on this file since 5625 was 5625, checked in by aborella, 7 months ago

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