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

Last change on this file since 5242 was 5240, checked in by Laurent Fairhead, 39 hours ago

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