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

Last change on this file since 5691 was 5691, checked in by aborella, 3 weeks ago

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