source: LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.f90 @ 5423

Last change on this file since 5423 was 5406, checked in by evignon, 2 weeks ago

changement de la parametrisation de la dissipation des cirrus et ajouts de commentaires.

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