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

Last change on this file since 5456 was 5456, checked in by aborella, 24 hours ago

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