source: LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.F90 @ 5209

Last change on this file since 5209 was 5204, checked in by Laurent Fairhead, 7 weeks ago

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