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

Last change on this file since 5450 was 5450, checked in by aborella, 7 hours ago

Merge with trunk
File size: 46.4 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, Povera, 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    !--Initialisation
296    issrfra(i)  = 0.
297    qissr(i)    = 0.
298
299    !--If the temperature is higher than the threshold below which
300    !--there is no liquid in the gridbox, we activate the usual scheme
301    !--(generalised lognormal from Bony and Emanuel 2001)
302    !--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for
303    !--all clouds, and the lognormal scheme is not activated
304    IF ( ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ) THEN
305
306      pdf_std   = ratqs(i) * qtot(i)
307      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) )
308      pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
309      pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
310      pdf_b     = pdf_k / (2. * SQRT(2.))
311      pdf_e1    = pdf_a - pdf_b
312      pdf_e1    = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
313      pdf_e1    = 1. - ERF(pdf_e1)
314      pdf_e2    = pdf_a + pdf_b
315      pdf_e2    = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
316      pdf_e2    = 1. - ERF(pdf_e2)
317
318
319      IF ( pdf_e1 .LT. eps ) THEN
320        cldfra(i) = 0.
321        qincld(i) = qsat(i)
322        qvc(i) = 0.
323      ELSE
324        cldfra(i) = 0.5 * pdf_e1
325        qincld(i) = qtot(i) * pdf_e2 / pdf_e1
326        qvc(i) = qsat(i) * cldfra(i)
327      ENDIF
328
329    !--If the temperature is lower than temp_nowater, we use the new
330    !--condensation scheme that allows for ice supersaturation
331    ELSE
332
333      !--Initialisation
334      IF ( temp(i) .GT. temp_nowater ) THEN
335        !--If the air mass is warm (liquid water can exist),
336        !--all the memory is lost and the scheme becomes statistical,
337        !--i.e., the sublimation and mixing processes are deactivated,
338        !--and the condensation process is slightly adapted
339        !--This can happen only if ok_weibull_warm_clouds = .TRUE.
340        ! AB WARM CLOUD
341        !cldfra(i) = 0.
342        !qcld(i)   = 0.
343        !qvc(i)    = 0.
344        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
345        qcld(i)   = MAX(0., MIN(qtot(i), ql_seri(i) + qi_seri(i) + rvc_seri(i) * qtot(i)))
346        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
347        ok_warm_cloud = .TRUE.
348      ELSE
349        !--The following barriers ensure that the traced cloud properties
350        !--are consistent. In some rare cases, i.e. the cloud water vapor
351        !--can be greater than the total water in the gridbox
352        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
353        qcld(i)   = MAX(0., MIN(qtot(i), rvc_seri(i) * qtot(i) + qi_seri(i)))
354        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
355        ok_warm_cloud = .FALSE.
356      ENDIF
357
358      dcf_sub(i)  = 0.
359      dqi_sub(i)  = 0.
360      dqvc_sub(i) = 0.
361      dqi_adj(i)  = 0.
362      dqvc_adj(i) = 0.
363      dcf_con(i)  = 0.
364      dqi_con(i)  = 0.
365      dqvc_con(i) = 0.
366      dcf_mix(i)  = 0.
367      dqi_mix(i)  = 0.
368      dqvc_mix(i) = 0.
369
370      !--Initialisation of the cell properties
371      !--Dry density [kg/m3]
372      rho = pplay(i) / temp(i) / RD
373      !--Dry air mass [kg/m2]
374      rhodz = ( paprsdn(i) - paprsup(i) ) / RG
375      !--Cell thickness [m]
376      dz = rhodz / rho
377      !--Cell volume [m3]
378      !V_cell = dz * cell_area(i)
379      !--Cell dry air mass [kg]
380      !M_cell = rhodz * cell_area(i)
381
382
383      !-------------------------------------------------------------------
384      !--    SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD    --
385      !-------------------------------------------------------------------
386     
387      !--If there is a cloud
388      IF ( cldfra(i) .GT. eps ) THEN
389       
390        qvapincld = qvc(i) / cldfra(i)
391        qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
392       
393        !--If the ice water content is too low, the cloud is purely sublimated
394        !--Most probably, we advected a cloud with no ice water content (possible
395        !--if the entire cloud precipited for example)
396        IF ( qiceincld .LT. eps ) THEN
397          dcf_sub(i) = - cldfra(i)
398          dqvc_sub(i) = - qvc(i)
399          dqi_sub(i) = - ( qcld(i) - qvc(i) )
400
401          cldfra(i) = 0.
402          qcld(i) = 0.
403          qvc(i) = 0.
404
405        !--Else, the cloud is adjusted and sublimated
406        ELSE
407
408          !--The vapor in cloud cannot be higher than the
409          !--condensation threshold
410          qvapincld = MIN(qvapincld, gamma_cond(i) * qsat(i))
411          qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
412
413          ! AB WARM CLOUD
414          !IF ( ok_unadjusted_clouds ) THEN
415          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
416            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
417                                        pplay(i), dtime, qvapincld_new)
418            IF ( qvapincld_new .EQ. 0. ) THEN
419              !--If all the ice has been sublimated, we sublimate
420              !--completely the cloud and do not activate the sublimation
421              !--process
422              !--Tendencies and diagnostics
423              dcf_sub(i) = - cldfra(i)
424              dqvc_sub(i) = - qvc(i)
425              dqi_sub(i) = - ( qcld(i) - qvc(i) )
426
427              cldfra(i) = 0.
428              qcld(i) = 0.
429              qvc(i) = 0.
430            ENDIF
431          ELSE
432            !--We keep the saturation adjustment hypothesis, and the vapor in the
433            !--cloud is set equal to the saturation vapor
434            qvapincld_new = qsat(i)
435          ENDIF ! ok_unadjusted_clouds
436         
437          !--Adjustment of the IWC to the new vapor in cloud
438          !--(this can be either positive or negative)
439          dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
440          dqi_adj(i) = - dqvc_adj(i)
441
442          !--Add tendencies
443          !--The vapor in the cloud is updated, but not qcld as it is constant
444          !--through this process, as well as cldfra which is unmodified
445          qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
446         
447
448          !------------------------------------
449          !--    DISSIPATION OF THE CLOUD    --
450          !------------------------------------
451
452          !--If the vapor in cloud is below vapor needed for the cloud to survive
453          IF ( ( qvapincld .LT. qvapincld_new ) .OR. ( iflag_cloud_sublim_pdf .GE. 4 ) ) THEN
454            !--Sublimation of the subsaturated cloud
455            !--iflag_cloud_sublim_pdf selects the PDF of the ice water content
456            !--to use.
457            !--iflag = 1 --> uniform distribution
458            !--iflag = 2 --> exponential distribution
459            !--iflag = 3 --> gamma distribution (Karcher et al 2018)
460
461            IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN
462              !--Uniform distribution starting at qvapincld
463              pdf_e1 = 1. / ( 2. * qiceincld )
464
465              dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
466              dqt_sub    = - cldfra(i) * ( qvapincld_new**2 - qvapincld**2 ) &
467                                       * pdf_e1 / 2.
468
469            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN
470              !--Exponential distribution starting at qvapincld
471              pdf_alpha = 1. / qiceincld
472              pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) )
473
474              dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 )
475              dqt_sub    = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha &
476                                       + qvapincld - qvapincld_new * pdf_e1 )
477
478            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 3 ) THEN
479              !--Gamma distribution starting at qvapincld
480              pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
481              pdf_y = pdf_alpha * ( qvapincld_new - qvapincld )
482              pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
483              pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
484
485              dcf_sub(i) = - cldfra(i) * pdf_e1
486              dqt_sub    = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
487
488            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 4 ) THEN
489              !--Normal distribution around qvapincld + qiceincld with width
490              !--constant in the RHi space
491              pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
492                    / ( std_subl_pdf_lscp / 100. * qsat(i))
493              pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
494              pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )
495
496              dcf_sub(i) = - cldfra(i) * pdf_e1
497              dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
498                                         - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
499                                                                                       
500            ENDIF
501
502            !--Tendencies and diagnostics
503            dqvc_sub(i) = dqt_sub
504
505            !--Add tendencies
506            cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
507            qcld(i)   = MAX(0., qcld(i) + dqt_sub)
508            qvc(i)    = MAX(0., qvc(i) + dqvc_sub(i))
509
510          ENDIF ! qvapincld .LT. qvapincld_new
511
512        ENDIF   ! qiceincld .LT. eps
513      ENDIF ! cldfra(i) .GT. eps
514
515
516      !--------------------------------------------------------------------------
517      !--    CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS    --
518      !--------------------------------------------------------------------------
519      !--This section relies on a distribution of water in the clear-sky region of
520      !--the mesh.
521
522      !--If there is a clear-sky region
523      IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
524
525        !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
526        qclr = qtot(i) - qcld(i)
527
528        !--New PDF
529        rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.
530
531        !--Calculation of the properties of the PDF
532        !--Parameterization from IAGOS observations
533        !--pdf_e1 and pdf_e2 will be reused below
534
535        !--Coefficient for standard deviation:
536        !--  tuning coef * (clear sky area**0.25) * (function of temperature)
537        !pdf_e1 = beta_pdf_lscp &
538        !       * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
539        !       * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
540        !--  tuning coef * (cell area**0.25) * (function of temperature)
541        pdf_e1 = beta_pdf_lscp * ( ( 1. - cldfra(i) ) * cell_area(i) )**0.25 &
542                * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
543        IF ( rhl_clr .GT. 50. ) THEN
544          pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp
545        ELSE
546          pdf_std = pdf_e1 * rhl_clr / 50.
547        ENDIF
548        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
549        pdf_alpha = EXP( rhl_clr / 100. ) * pdf_e3
550        pdf_alpha = MIN(10., pdf_alpha)
551       
552        IF ( ok_warm_cloud ) THEN
553          !--If the statistical scheme is activated, the standard deviation is adapted
554          !--to depend on the pressure level. It is multiplied by ratqs, so that near the
555          !--surface std is almost 0, and upper than about 450 hPa the std is left untouched
556          pdf_std = pdf_std * ratqs(i)
557        ENDIF
558
559        pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
560        pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2 ))
561        pdf_loc = rhl_clr - pdf_scale * pdf_e2
562
563        !--Calculation of the newly condensed water and fraction (pronostic)
564        !--Integration of the clear sky PDF between gamma_cond*qsat and +inf
565        !--NB. the calculated values are clear-sky averaged
566
567        pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
568        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
569        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
570        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
571        cf_cond = EXP( - pdf_y )
572        qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.
573
574        ! AB WARM CLOUD
575        !IF ( ok_warm_cloud ) THEN
576        !  !--If the statistical scheme is activated, the calculated increase is equal
577        !  !--to the cloud fraction, we assume saturation adjustment, and the
578        !  !--condensation process stops
579        !  cldfra(i) = cf_cond
580        !  qcld(i) = qt_cond
581        !  qvc(i) = cldfra(i) * qsat(i)
582
583        !ELSEIF ( cf_cond .GT. eps ) THEN
584        IF ( cf_cond .GT. eps ) THEN
585
586          dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond
587          dqt_con = ( 1. - cldfra(i) ) * qt_cond
588         
589          !--Barriers
590          dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i))
591          dqt_con = MIN(dqt_con, qclr)
592
593
594          ! AB WARM CLOUD
595          !IF ( ok_unadjusted_clouds ) THEN
596          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
597            !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
598            !--the new vapor qvapincld. The timestep is divided by two because we do not
599            !--know when the condensation occurs
600            qvapincld = gamma_cond(i) * qsat(i)
601            qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
602            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
603                                        pplay(i), dtime / 2., qvapincld_new)
604            qvapincld = qvapincld_new
605          ELSE
606            !--We keep the saturation adjustment hypothesis, and the vapor in the
607            !--newly formed cloud is set equal to the saturation vapor.
608            qvapincld = qsat(i)
609          ENDIF
610
611          !--Tendency on cloud vapor and diagnostic
612          dqvc_con(i) = qvapincld * dcf_con(i)
613          dqi_con(i) = dqt_con - dqvc_con(i)
614
615          !--Add tendencies
616          cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
617          qcld(i)   = MIN(qtot(i), qcld(i) + dqt_con)
618          qvc(i)    = MIN(qcld(i), qvc(i) + dqvc_con(i))
619
620        ENDIF ! ok_warm_cloud, cf_cond .GT. eps
621       
622        !--We then calculate the part that is greater than qsat
623        !--and lower than gamma_cond * qsat, which is the ice supersaturated region
624        pdf_x = qsat(i) / qsatl(i) * 100.
625        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
626        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
627        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
628        issrfra(i) = EXP( - pdf_y ) * ( 1. - cldfra(i) )
629        qissr(i) = ( pdf_e3 * ( 1. - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
630
631        issrfra(i) = issrfra(i) - dcf_con(i)
632        qissr(i) = qissr(i) - dqvc_con(i) - dqi_con(i)
633
634      ENDIF ! ( 1. - cldfra(i) ) .GT. eps
635
636      !--Calculation of the subsaturated clear sky fraction and water
637      subfra(i) = 1. - cldfra(i) - issrfra(i)
638      qsub(i) = qtot(i) - qcld(i) - qissr(i)
639
640
641      !--------------------------------------
642      !--           CLOUD MIXING           --
643      !--------------------------------------
644      !--This process mixes the cloud with its surroundings: the subsaturated clear sky,
645      !--and the supersaturated clear sky. It is activated if the cloud is big enough,
646      !--but does not cover the entire mesh.
647      !
648      ! AB WARM CLOUD
649      !IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
650      !        .AND. .NOT. ok_warm_cloud ) THEN
651      IF ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) THEN
652
653        !--Initialisation
654        dcf_mix_sub   = 0.
655        dqt_mix_sub   = 0.
656        dqvc_mix_sub  = 0.
657        dcf_mix_issr  = 0.
658        dqt_mix_issr  = 0.
659        dqvc_mix_issr = 0.
660
661
662        !-- PART 1 - TURBULENT DIFFUSION
663
664        !--Clouds within the mesh are assumed to be ellipses. The length of the
665        !--semi-major axis is a and the length of the semi-minor axis is b.
666        !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer.
667        !--In particular, it is 0 if cldfra = 1.
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        !--P / a is a function of b / a only, that we can calculate
687        !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
688        Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
689        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
690        !--based on observations.
691        !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
692        !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
693        !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
694        !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
695        a_mix = Povera / coef_mixing_lscp / RPI / ( 1. - cldfra(i) ) / bovera
696        !--and finally,
697        !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
698        N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) ) * cell_area(i) &
699                  / Povera / a_mix
700
701        !--The time required for turbulent diffusion to homogenize a region of size
702        !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
703        !--We compute L_mix and assume that the cloud is mixed over this length
704        L_mix = SQRT( dtime**3 * pbl_eps(i) )
705        !--The mixing length cannot be greater than the semi-minor axis. In this case,
706        !--the entire cloud is mixed.
707        L_mix = MIN(L_mix, a_mix * bovera)
708
709        !--The fraction of clear sky mixed is
710        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
711        envfra_mix = N_cld_mix * RPI / cell_area(i) &
712                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
713        !--The fraction of cloudy sky mixed is
714        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
715        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
716                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
717
718
719        !-- PART 2 - SHEARING
720
721        !--The clouds are then sheared. We keep the shape and number
722        !--assumptions from before. The clouds are sheared along their
723        !--semi-major axis (a_mix), on the entire cell heigh dz.
724        !--The increase in size is
725        L_shear = coef_shear_lscp * shear(i) * dz * dtime
726        !--therefore, the fraction of clear sky mixed is
727        !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
728        !--and the fraction of cloud mixed is
729        !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area
730        !--(note that they are equal)
731        shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i)
732        !--and the environment and cloud mixed fractions are the same,
733        !--which we add to the previous calculated mixed fractions.
734        !--We therefore assume that the sheared clouds and the turbulent
735        !--mixed clouds are different.
736        envfra_mix = envfra_mix + shear_fra
737        cldfra_mix = cldfra_mix + shear_fra
738
739
740        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
741
742        !--The environment fraction is allocated to subsaturated sky or supersaturated sky,
743        !--according to the factor sigma_mix. This is computed as the ratio of the
744        !--subsaturated sky fraction to the environment fraction, corrected by a factor
745        !--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
746        !--supersaturated sky is favoured. Physically, this means that it is more likely
747        !--to have supersaturated sky around the cloud than subsaturated sky.
748        sigma_mix = subfra(i) / ( subfra(i) + chi_mixing_lscp * issrfra(i) )
749        subfra_mix = MIN( sigma_mix * envfra_mix, subfra(i) )
750        issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra(i) )
751        cldfra_mix = MIN( cldfra_mix, cldfra(i) )
752
753        !--First, we mix the subsaturated sky (subfra_mix) and the cloud close
754        !--to this fraction (sigma_mix * cldfra_mix).
755        IF ( subfra(i) .GT. eps ) THEN
756          !--We compute the total humidity in the mixed air, which
757          !--can be either sub- or supersaturated.
758          qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
759                    + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
760                    / ( subfra_mix + cldfra_mix * sigma_mix )
761
762          ! AB WARM CLOUD
763          !IF ( ok_unadjusted_clouds ) THEN
764          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
765            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * sigma_mix / cldfra(i) &
766                      / ( subfra_mix + cldfra_mix * sigma_mix )
767            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
768                                        pplay(i), dtime, qvapincld_new)
769            IF ( qvapincld_new .EQ. 0. ) THEN
770              !--If all the ice has been sublimated, we sublimate
771              !--completely the cloud and do not activate the sublimation
772              !--process
773              !--Tendencies and diagnostics
774              dcf_mix_sub = - sigma_mix * cldfra_mix
775              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
776              dqvc_mix_sub = dcf_mix_sub * qvc(i) / cldfra(i)
777            ELSE
778              dcf_mix_sub = subfra_mix
779              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
780              dqvc_mix_sub = dcf_mix_sub * qvapincld_new
781            ENDIF
782          ELSE
783            IF ( qvapinmix .GT. qsat(i) ) THEN
784              !--If the mixed air is supersaturated, we condense the subsaturated
785              !--region which was mixed.
786              dcf_mix_sub = subfra_mix
787              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
788              dqvc_mix_sub = dcf_mix_sub * qsat(i)
789            ELSE
790              !--Else, we sublimate the cloud which was mixed.
791              dcf_mix_sub = - sigma_mix * cldfra_mix
792              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
793              dqvc_mix_sub = dcf_mix_sub * qsat(i)
794            ENDIF
795          ENDIF ! ok_unadjusted_clouds
796        ENDIF ! subfra .GT. eps
797   
798        !--We then mix the supersaturated sky (issrfra_mix) and the cloud,
799        !--for which the mixed air is always supersatured, therefore
800        !--the cloud necessarily expands
801        IF ( issrfra(i) .GT. eps ) THEN
802
803          ! AB WARM CLOUD
804          !IF ( ok_unadjusted_clouds ) THEN
805          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
806            qvapinmix = ( qissr(i) * issrfra_mix / issrfra(i) &
807                      + qcld(i) * cldfra_mix * (1. - sigma_mix) / cldfra(i) ) &
808                      / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
809            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * (1. - sigma_mix) / cldfra(i) &
810                      / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
811            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
812                                        pplay(i), dtime, qvapincld_new)
813            dcf_mix_issr = issrfra_mix
814            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
815            dqvc_mix_issr = dcf_mix_issr * qvapincld_new
816          ELSE
817            !--In this case, the additionnal vapor condenses
818            dcf_mix_issr = issrfra_mix
819            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
820            dqvc_mix_issr = dcf_mix_issr * qsat(i)
821          ENDIF ! ok_unadjusted_clouds
822
823
824        ENDIF ! issrfra .GT. eps
825
826        !--Sum up the tendencies from subsaturated sky and supersaturated sky
827        dcf_mix(i)  = dcf_mix_sub  + dcf_mix_issr
828        dqt_mix     = dqt_mix_sub  + dqt_mix_issr
829        dqvc_mix(i) = dqvc_mix_sub + dqvc_mix_issr
830        dqi_mix(i)  = dqt_mix - dqvc_mix(i)
831
832        !--Add tendencies
833        issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
834        qissr(i)   = MAX(0., qissr(i) - dqt_mix_issr)
835        cldfra(i)  = MAX(0., MIN(1., cldfra(i) + dcf_mix(i)))
836        qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqt_mix))
837        qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))
838
839      ENDIF ! ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
840
841
842      !----------------------------------------
843      !--       CONTRAILS AND AVIATION       --
844      !----------------------------------------
845
846      !--Add a source of cirrus from aviation contrails
847      !IF ( ok_plane_contrail ) THEN
848      !  dcf_avi(i) = 0.
849      !  dqi_avi(i) = 0.
850      !  dqvc_avi(i) = 0.
851      !  ! TODO implement ok_unadjusted_clouds
852      !  IF ( issrfra(i) .GT. eps ) THEN
853      !    contrail_fra = MIN(1., flight_m(i,k) * dtime * contrail_cross_section / V_cell)
854      !    dcf_avi(i) = issrfra(i) * contrail_fra
855      !    dqt_avi = dcf_avi(i) * qissr(i) / issrfra(i)
856      !    dqvc_avi(i) = qsat(i) * dcf_avi(i)
857      !   
858      !    !--Add tendencies
859      !    cldfra(i) = cldfra(i) + dcf_avi(i)
860      !    issrfra(i) = issrfra(i) - dcf_avi(i)
861      !    qcld(i) = qcld(i) + dqt_avi
862      !    qvc(i) = qvc(i) + dqvc_avi(i)
863      !    qissr(i) = qissr(i) - dqt_avi
864
865      !    !--Diagnostics
866      !    dqi_avi(i) = dqt_avi - qsat(i) * dcf_avi(i)
867      !  ENDIF
868      !  dcf_avi(i)  = dcf_avi(i)  / dtime
869      !  dqi_avi(i)  = dqi_avi(i)  / dtime
870      !  dqvc_avi(i) = dqvc_avi(i) / dtime
871      !ENDIF
872
873
874
875      !-------------------------------------------
876      !--       FINAL BARRIERS AND OUTPUTS      --
877      !-------------------------------------------
878
879      IF ( cldfra(i) .LT. eps ) THEN
880        !--If the cloud is too small, it is sublimated.
881        cldfra(i) = 0.
882        qcld(i)   = 0.
883        qvc(i)    = 0.
884        qincld(i) = qsat(i)
885      ELSE
886        qincld(i) = qcld(i) / cldfra(i)
887      ENDIF ! cldfra .LT. eps
888
889      !--Diagnostics
890      dcf_sub(i)  = dcf_sub(i)  / dtime
891      dcf_con(i)  = dcf_con(i)  / dtime
892      dcf_mix(i)  = dcf_mix(i)  / dtime
893      dqi_adj(i)  = dqi_adj(i)  / dtime
894      dqi_sub(i)  = dqi_sub(i)  / dtime
895      dqi_con(i)  = dqi_con(i)  / dtime
896      dqi_mix(i)  = dqi_mix(i)  / dtime
897      dqvc_adj(i) = dqvc_adj(i) / dtime
898      dqvc_sub(i) = dqvc_sub(i) / dtime
899      dqvc_con(i) = dqvc_con(i) / dtime
900      dqvc_mix(i) = dqvc_mix(i) / dtime
901
902    ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
903
904  ENDIF   ! end keepgoing
905
906ENDDO     ! end loop on i
907
908END SUBROUTINE condensation_ice_supersat
909!**********************************************************************************
910
911SUBROUTINE deposition_sublimation( &
912    qvapincld, qiceincld, temp, qsat, pplay, dtime, &
913    qvapincld_new)
914
915USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, RD, EPS_W
916USE lmdz_lscp_ini,        ONLY: eps
917USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
918
919REAL,     INTENT(IN) :: qvapincld     !
920REAL,     INTENT(IN) :: qiceincld     !
921REAL,     INTENT(IN) :: temp          !
922REAL,     INTENT(IN) :: qsat          !
923REAL,     INTENT(IN) :: pplay         !
924REAL,     INTENT(IN) :: dtime         ! time step [s]
925REAL,     INTENT(OUT):: qvapincld_new !
926
927!
928! for unadjusted clouds
929REAL :: qice_ratio
930REAL :: pres_sat, rho, kappa
931REAL :: air_thermal_conduct, water_vapor_diff
932REAL :: iwc
933REAL :: iwc_log_inf100, iwc_log_sup100
934REAL :: iwc_inf100, alpha_inf100, coef_inf100
935REAL :: mu_sup100, sigma_sup100, coef_sup100
936REAL :: Dm_ice, rm_ice
937
938!--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
939!--hypothesis is lost, and the vapor in the cloud is purely prognostic.
940!
941!--The deposition equation is
942!-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
943!--from Lohmann et al. (2016), where
944!--alpha is the deposition coefficient [-]
945!--mi is the mass of one ice crystal [kg]
946!--C is the capacitance of an ice crystal [m]
947!--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
948!--R_v is the specific gas constant for humid air [J/kg/K]
949!--T is the temperature [K]
950!--esi is the saturation pressure w.r.t. ice [Pa]
951!--Dv is the diffusivity of water vapor [m2/s]
952!--Ls is the specific latent heat of sublimation [J/kg/K]
953!--ka is the thermal conductivity of dry air [J/m/s/K]
954!
955!--alpha is a coefficient to take into account the fact that during deposition, a water
956!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
957!--coherent (with the same structure). It has no impact for sublimation.
958!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
959!--during deposition, and alpha = 1. during sublimation.
960!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
961!-- C = capa_cond_cirrus * rm_ice
962!
963!--We have qice = Nice * mi, where Nice is the ice crystal
964!--number concentration per kg of moist air
965!--HYPOTHESIS 1: the ice crystals are spherical, therefore
966!-- mi = 4/3 * pi * rm_ice**3 * rho_ice
967!--HYPOTHESIS 2: the ice crystals are monodisperse with the
968!--initial radius rm_ice_0.
969!--NB. this is notably different than the assumption
970!--of a distributed qice in the cloud made in the sublimation process;
971!--should it be consistent?
972!
973!--As the deposition process does not create new ice crystals,
974!--and because we assume a same rm_ice value for all crystals
975!--therefore the sublimation process does not destroy ice crystals
976!--(or, in a limit case, it destroys all ice crystals), then
977!--Nice is a constant during the sublimation/deposition process.
978!-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
979!
980!--The deposition equation then reads:
981!-- 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
982!-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
983!--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
984!--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
985!-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
986!--  *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 )
987!--and we have
988!-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
989!-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
990!--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
991!
992!--This system of equations can be resolved with an exact
993!--explicit numerical integration, having one variable resolved
994!--explicitly, the other exactly. The exactly resolved variable is
995!--the one decreasing, so it is qvc if the process is
996!--condensation, qi if it is sublimation.
997!
998!--kappa is computed as an initialisation constant, as it depends only
999!--on temperature and other pre-computed values
1000pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
1001!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
1002air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
1003!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
1004water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
1005kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
1006      * ( RV * temp / water_vapor_diff / pres_sat &
1007        + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
1008!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
1009
1010!--Dry density [kg/m3]
1011rho = pplay / temp / RD
1012
1013!--Here, the initial vapor in the cloud is qvapincld, and we compute
1014!--the new vapor qvapincld_new
1015
1016!--rm_ice formula from McFarquhar and Heymsfield (1997)
1017iwc = qiceincld * rho * 1e3
1018iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
1019iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
1020iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
1021
1022alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
1023coef_inf100 = iwc_inf100 * alpha_inf100**3 / 120.
1024
1025mu_sup100 = ( 5.2 + 0.0013 * ( temp - RTT ) ) &
1026          + ( 0.026 - 1.2E-3 * ( temp - RTT ) ) * iwc_log_sup100
1027sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp - RTT ) ) &
1028             + ( 0.018 - 2.1E-4 * ( temp - RTT ) ) * iwc_log_sup100
1029coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2 )
1030
1031Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2 ) &
1032       * coef_sup100 ) / ( coef_inf100 + coef_sup100 )
1033rm_ice = Dm_ice / 2. * 1.E-6
1034
1035IF ( qvapincld .GE. qsat ) THEN
1036  !--If the cloud is initially supersaturated
1037  !--Exact explicit integration (qvc exact, qice explicit)
1038  qvapincld_new = qsat + ( qvapincld - qsat ) &
1039                * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
1040ELSE
1041  !--If the cloud is initially subsaturated
1042  !--Exact explicit integration (qice exact, qvc explicit)
1043  !--The barrier is set so that the resulting vapor in cloud
1044  !--cannot be greater than qsat
1045  !--qice_ratio is the ratio between the new ice content and
1046  !--the old one, it is comprised between 0 and 1
1047  qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) )
1048
1049  IF ( qice_ratio .LT. 0. ) THEN
1050    !--The new vapor in cloud is set to 0 so that the
1051    !--sublimation process does not activate
1052    qvapincld_new = 0.
1053  ELSE
1054    !--Else, the sublimation process is activated with the
1055    !--diagnosed new cloud water vapor
1056    !--The new vapor in the cloud is increased with the
1057    !--sublimated ice
1058    qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 )
1059    !--The new vapor in the cloud cannot be greater than qsat
1060    qvapincld_new = MIN(qvapincld_new, qsat)
1061  ENDIF ! qice_ratio .LT. 0.
1062ENDIF ! qvapincld .GT. qsat
1063
1064END SUBROUTINE deposition_sublimation
1065
1066END MODULE lmdz_lscp_condensation
Note: See TracBrowser for help on using the repository browser.