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

Last change on this file since 5452 was 5452, checked in by aborella, 5 days ago

First implementation of the contrails parameterisation
Lacks the emission of H2O + the impact on radiative transfer
File size: 44.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      rcont_seri, flight_dist, flight_h2o, contfra, &
102      Tcritcont, qcritcont, potcontfraP, potcontfraNP, 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: ok_plane_contrail
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, chi_mixing_lscp
127
128USE lmdz_aviation,        ONLY: contrails_mixing, contrails_formation_evolution
129
130IMPLICIT NONE
131
132!
133! Input
134!
135INTEGER,  INTENT(IN)                     :: klon          ! number of horizontal grid points
136REAL,     INTENT(IN)                     :: dtime         ! time step [s]
137REAL,     INTENT(IN)                     :: missing_val   ! missing value for output
138!
139REAL,     INTENT(IN)   , DIMENSION(klon) :: pplay         ! layer pressure [Pa]
140REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsdn       ! pressure at the lower interface [Pa]
141REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsup       ! pressure at the upper interface [Pa]
142REAL,     INTENT(IN)   , DIMENSION(klon) :: cf_seri       ! cloud fraction [-]
143REAL,     INTENT(IN)   , DIMENSION(klon) :: rvc_seri      ! gridbox-mean water vapor in cloud [kg/kg]
144REAL,     INTENT(IN)   , DIMENSION(klon) :: ql_seri       ! specific liquid water content [kg/kg]
145REAL,     INTENT(IN)   , DIMENSION(klon) :: qi_seri       ! specific ice water content [kg/kg]
146REAL,     INTENT(IN)   , DIMENSION(klon) :: shear         ! vertical shear [s-1]
147REAL,     INTENT(IN)   , DIMENSION(klon) :: pbl_eps       ! TKE dissipation [m2/s3]
148REAL,     INTENT(IN)   , DIMENSION(klon) :: cell_area     ! cell area [m2]
149REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
150REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot          ! total specific humidity (without precip) [kg/kg]
151REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
152REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
153REAL,     INTENT(IN)   , DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
154LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
155!
156! Input for aviation
157!
158REAL,     INTENT(INOUT), DIMENSION(klon) :: rcont_seri    ! ratio of contrails fraction to total cloud fraction [-]
159REAL,     INTENT(IN),    DIMENSION(klon) :: flight_dist   !
160REAL,     INTENT(IN),    DIMENSION(klon) :: flight_h2o    !
161!
162!  Output
163!  NB. cldfra and qincld should be outputed as cf_seri and qi_seri,
164!      or as tendencies (maybe in the future)
165!  NB. those are in INOUT because of the convergence loop on temperature
166!      (in some cases, the values are not re-computed) but the values
167!      are never used explicitely
168!
169REAL,     INTENT(INOUT), DIMENSION(klon) :: cldfra   ! cloud fraction [-]
170REAL,     INTENT(INOUT), DIMENSION(klon) :: qincld   ! cloud-mean in-cloud total specific water [kg/kg]
171REAL,     INTENT(INOUT), DIMENSION(klon) :: qvc      ! gridbox-mean vapor in the cloud [kg/kg]
172REAL,     INTENT(INOUT), DIMENSION(klon) :: issrfra  ! ISSR fraction [-]
173REAL,     INTENT(INOUT), DIMENSION(klon) :: qissr    ! gridbox-mean ISSR specific water [kg/kg]
174!
175!  Diagnostics for condensation and ice supersaturation
176!  NB. idem for the INOUT
177!
178REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_sub  ! cloud fraction tendency because of sublimation [s-1]
179REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_con  ! cloud fraction tendency because of condensation [s-1]
180REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_mix  ! cloud fraction tendency because of cloud mixing [s-1]
181REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_adj  ! specific ice content tendency because of temperature adjustment [kg/kg/s]
182REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_sub  ! specific ice content tendency because of sublimation [kg/kg/s]
183REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_con  ! specific ice content tendency because of condensation [kg/kg/s]
184REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_mix  ! specific ice content tendency because of cloud mixing [kg/kg/s]
185REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
186REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s]
187REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s]
188REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
189!
190!  Diagnostics for aviation
191!  NB. idem for the INOUT
192!
193REAL,     INTENT(INOUT), DIMENSION(klon) :: contfra      ! contrail fraction [-]
194REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
195REAL,     INTENT(INOUT), DIMENSION(klon) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
196REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraP  ! potential persistent contrail fraction [-]
197REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-]
198REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_avi      ! cloud fraction tendency because of aviation [s-1]
199REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_avi      ! specific ice content tendency because of aviation [kg/kg/s]
200REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_avi     ! specific cloud water vapor tendency because of aviation [kg/kg/s]
201!
202! Local
203!
204INTEGER :: i
205LOGICAL :: ok_warm_cloud
206REAL, DIMENSION(klon) :: qcld, qzero
207!
208! for lognormal
209REAL :: pdf_std, pdf_k, pdf_delta
210REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
211!
212! for unadjusted clouds
213REAL :: qvapincld, qvapincld_new
214REAL :: qiceincld
215!
216! for sublimation
217REAL :: pdf_alpha
218REAL :: dqt_sub
219!
220! for condensation
221REAL, DIMENSION(klon) :: qsatl, dqsatl
222REAL :: clrfra, qclr, sl_clr, rhl_clr
223REAL :: pdf_ratqs, pdf_skew, pdf_scale, pdf_loc
224REAL :: pdf_x, pdf_y, pdf_T
225REAL :: pdf_e3, pdf_e4
226REAL :: cf_cond, qt_cond, dqt_con
227!
228! for mixing
229REAL, DIMENSION(klon) :: subfra, qsub
230REAL :: dqt_mix_sub, dqt_mix_issr
231REAL :: dcf_mix_sub, dcf_mix_issr
232REAL :: dqvc_mix_sub, dqvc_mix_issr
233REAL :: dqt_mix
234REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix
235REAL :: envfra_mix, cldfra_mix
236REAL :: L_shear, shear_fra
237REAL :: sigma_mix, issrfra_mix, subfra_mix, qvapinmix
238!
239! for cell properties
240REAL :: rho, rhodz, dz
241REAL :: V_cell, M_cell
242
243qzero(:) = 0.
244
245!--Calculation of qsat w.r.t. liquid
246CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsatl)
247
248!
249!--Loop on klon
250!
251DO i = 1, klon
252
253  !--If a new calculation of the condensation is needed,
254  !--i.e., temperature has not yet converged (or the cloud is
255  !--formed elsewhere)
256  IF (keepgoing(i)) THEN
257
258    !--Initialisation
259    issrfra(i)  = 0.
260    qissr(i)    = 0.
261
262    !--If the temperature is higher than the threshold below which
263    !--there is no liquid in the gridbox, we activate the usual scheme
264    !--(generalised lognormal from Bony and Emanuel 2001)
265    !--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for
266    !--all clouds, and the lognormal scheme is not activated
267    IF ( ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ) THEN
268
269      pdf_std   = ratqs(i) * qtot(i)
270      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) )
271      pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
272      pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
273      pdf_b     = pdf_k / (2. * SQRT(2.))
274      pdf_e1    = pdf_a - pdf_b
275      pdf_e1    = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
276      pdf_e1    = 1. - ERF(pdf_e1)
277      pdf_e2    = pdf_a + pdf_b
278      pdf_e2    = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
279      pdf_e2    = 1. - ERF(pdf_e2)
280
281
282      IF ( pdf_e1 .LT. eps ) THEN
283        cldfra(i) = 0.
284        qincld(i) = qsat(i)
285        qvc(i) = 0.
286      ELSE
287        cldfra(i) = 0.5 * pdf_e1
288        qincld(i) = qtot(i) * pdf_e2 / pdf_e1
289        qvc(i) = qsat(i) * cldfra(i)
290      ENDIF
291
292    !--If the temperature is lower than temp_nowater, we use the new
293    !--condensation scheme that allows for ice supersaturation
294    ELSE
295
296      !--Initialisation
297      IF ( temp(i) .GT. temp_nowater ) THEN
298        !--If the air mass is warm (liquid water can exist),
299        !--all the memory is lost and the scheme becomes statistical,
300        !--i.e., the sublimation and mixing processes are deactivated,
301        !--and the condensation process is slightly adapted
302        !--This can happen only if ok_weibull_warm_clouds = .TRUE.
303        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
304        qcld(i)   = MAX(0., MIN(qtot(i), ql_seri(i) + qi_seri(i) + rvc_seri(i) * qtot(i)))
305        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
306        ok_warm_cloud = .TRUE.
307      ELSE
308        !--The following barriers ensure that the traced cloud properties
309        !--are consistent. In some rare cases, i.e. the cloud water vapor
310        !--can be greater than the total water in the gridbox
311        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
312        qcld(i)   = MAX(0., MIN(qtot(i), rvc_seri(i) * qtot(i) + qi_seri(i)))
313        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
314        ok_warm_cloud = .FALSE.
315      ENDIF
316
317      dcf_sub(i)  = 0.
318      dqi_sub(i)  = 0.
319      dqvc_sub(i) = 0.
320      dqi_adj(i)  = 0.
321      dqvc_adj(i) = 0.
322      dcf_con(i)  = 0.
323      dqi_con(i)  = 0.
324      dqvc_con(i) = 0.
325      dcf_mix(i)  = 0.
326      dqi_mix(i)  = 0.
327      dqvc_mix(i) = 0.
328
329      !--Initialisation of the cell properties
330      !--Dry density [kg/m3]
331      rho = pplay(i) / temp(i) / RD
332      !--Dry air mass [kg/m2]
333      rhodz = ( paprsdn(i) - paprsup(i) ) / RG
334      !--Cell thickness [m]
335      dz = rhodz / rho
336      !--Cell volume [m3]
337      V_cell = dz * cell_area(i)
338      !--Cell dry air mass [kg]
339      M_cell = rhodz * cell_area(i)
340
341
342      !-------------------------------------------------------------------
343      !--    SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD    --
344      !-------------------------------------------------------------------
345     
346      !--If there is a cloud
347      IF ( cldfra(i) .GT. eps ) THEN
348       
349        qvapincld = qvc(i) / cldfra(i)
350        qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
351       
352        !--If the ice water content is too low, the cloud is purely sublimated
353        !--Most probably, we advected a cloud with no ice water content (possible
354        !--if the entire cloud precipited for example)
355        IF ( qiceincld .LT. eps ) THEN
356          dcf_sub(i) = - cldfra(i)
357          dqvc_sub(i) = - qvc(i)
358          dqi_sub(i) = - ( qcld(i) - qvc(i) )
359
360          cldfra(i) = 0.
361          qcld(i) = 0.
362          qvc(i) = 0.
363
364        !--Else, the cloud is adjusted and sublimated
365        ELSE
366
367          !--The vapor in cloud cannot be higher than the
368          !--condensation threshold
369          qvapincld = MIN(qvapincld, gamma_cond(i) * qsat(i))
370          qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
371
372          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
373            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
374                                        pplay(i), dtime, qvapincld_new)
375            IF ( qvapincld_new .EQ. 0. ) THEN
376              !--If all the ice has been sublimated, we sublimate
377              !--completely the cloud and do not activate the sublimation
378              !--process
379              !--Tendencies and diagnostics
380              dcf_sub(i) = - cldfra(i)
381              dqvc_sub(i) = - qvc(i)
382              dqi_sub(i) = - ( qcld(i) - qvc(i) )
383
384              cldfra(i) = 0.
385              qcld(i) = 0.
386              qvc(i) = 0.
387            ENDIF
388          ELSE
389            !--We keep the saturation adjustment hypothesis, and the vapor in the
390            !--cloud is set equal to the saturation vapor
391            qvapincld_new = qsat(i)
392          ENDIF ! ok_unadjusted_clouds
393         
394          !--Adjustment of the IWC to the new vapor in cloud
395          !--(this can be either positive or negative)
396          dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
397          dqi_adj(i) = - dqvc_adj(i)
398
399          !--Add tendencies
400          !--The vapor in the cloud is updated, but not qcld as it is constant
401          !--through this process, as well as cldfra which is unmodified
402          qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
403         
404
405          !------------------------------------
406          !--    DISSIPATION OF THE CLOUD    --
407          !------------------------------------
408
409          !--If the vapor in cloud is below vapor needed for the cloud to survive
410          IF ( ( qvapincld .LT. qvapincld_new ) .OR. ( iflag_cloud_sublim_pdf .GE. 4 ) ) THEN
411            !--Sublimation of the subsaturated cloud
412            !--iflag_cloud_sublim_pdf selects the PDF of the ice water content
413            !--to use.
414            !--iflag = 1 --> uniform distribution
415            !--iflag = 2 --> exponential distribution
416            !--iflag = 3 --> gamma distribution (Karcher et al 2018)
417
418            IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN
419              !--Uniform distribution starting at qvapincld
420              pdf_e1 = 1. / ( 2. * qiceincld )
421
422              dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
423              dqt_sub    = - cldfra(i) * ( qvapincld_new**2 - qvapincld**2 ) &
424                                       * pdf_e1 / 2.
425
426            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN
427              !--Exponential distribution starting at qvapincld
428              pdf_alpha = 1. / qiceincld
429              pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) )
430
431              dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 )
432              dqt_sub    = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha &
433                                       + qvapincld - qvapincld_new * pdf_e1 )
434
435            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 3 ) THEN
436              !--Gamma distribution starting at qvapincld
437              pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
438              pdf_y = pdf_alpha * ( qvapincld_new - qvapincld )
439              pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
440              pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
441
442              dcf_sub(i) = - cldfra(i) * pdf_e1
443              dqt_sub    = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
444
445            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 4 ) THEN
446              !--Normal distribution around qvapincld + qiceincld with width
447              !--constant in the RHi space
448              pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
449                    / ( std_subl_pdf_lscp / 100. * qsat(i))
450              pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
451              pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )
452
453              dcf_sub(i) = - cldfra(i) * pdf_e1
454              dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
455                                         - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
456                                                                                       
457            ENDIF
458
459            !--Tendencies and diagnostics
460            dqvc_sub(i) = dqt_sub
461
462            !--Add tendencies
463            cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
464            qcld(i)   = MAX(0., qcld(i) + dqt_sub)
465            qvc(i)    = MAX(0., qvc(i) + dqvc_sub(i))
466
467          ENDIF ! qvapincld .LT. qvapincld_new
468
469        ENDIF   ! qiceincld .LT. eps
470      ENDIF ! cldfra(i) .GT. eps
471
472
473      !--------------------------------------------------------------------------
474      !--    CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS    --
475      !--------------------------------------------------------------------------
476      !--This section relies on a distribution of water in the clear-sky region of
477      !--the mesh.
478
479      !--If there is a clear-sky region
480      IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
481
482        !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
483        qclr = qtot(i) - qcld(i)
484
485        !--New PDF
486        rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.
487
488        !--Calculation of the properties of the PDF
489        !--Parameterization from IAGOS observations
490        !--pdf_e1 and pdf_e2 will be reused below
491
492        !--Coefficient for standard deviation:
493        !--  tuning coef * (clear sky area**0.25) * (function of temperature)
494        !pdf_e1 = beta_pdf_lscp &
495        !       * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
496        !       * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
497        !--  tuning coef * (cell area**0.25) * (function of temperature)
498        pdf_e1 = beta_pdf_lscp * ( ( 1. - cldfra(i) ) * cell_area(i) )**0.25 &
499                * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
500        IF ( rhl_clr .GT. 50. ) THEN
501          pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp
502        ELSE
503          pdf_std = pdf_e1 * rhl_clr / 50.
504        ENDIF
505        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
506        pdf_alpha = EXP( rhl_clr / 100. ) * pdf_e3
507        pdf_alpha = MIN(10., pdf_alpha)
508       
509        IF ( ok_warm_cloud ) THEN
510          !--If the statistical scheme is activated, the standard deviation is adapted
511          !--to depend on the pressure level. It is multiplied by ratqs, so that near the
512          !--surface std is almost 0, and upper than about 450 hPa the std is left untouched
513          pdf_std = pdf_std * ratqs(i)
514        ENDIF
515
516        pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
517        pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2 ))
518        pdf_loc = rhl_clr - pdf_scale * pdf_e2
519
520        !--Calculation of the newly condensed water and fraction (pronostic)
521        !--Integration of the clear sky PDF between gamma_cond*qsat and +inf
522        !--NB. the calculated values are clear-sky averaged
523
524        pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
525        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
526        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
527        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
528        cf_cond = EXP( - pdf_y )
529        qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.
530
531        IF ( cf_cond .GT. eps ) THEN
532
533          dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond
534          dqt_con = ( 1. - cldfra(i) ) * qt_cond
535         
536          !--Barriers
537          dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i))
538          dqt_con = MIN(dqt_con, qclr)
539
540          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
541            !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
542            !--the new vapor qvapincld. The timestep is divided by two because we do not
543            !--know when the condensation occurs
544            qvapincld = gamma_cond(i) * qsat(i)
545            qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
546            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
547                                        pplay(i), dtime / 2., qvapincld_new)
548            qvapincld = qvapincld_new
549          ELSE
550            !--We keep the saturation adjustment hypothesis, and the vapor in the
551            !--newly formed cloud is set equal to the saturation vapor.
552            qvapincld = qsat(i)
553          ENDIF
554
555          !--Tendency on cloud vapor and diagnostic
556          dqvc_con(i) = qvapincld * dcf_con(i)
557          dqi_con(i) = dqt_con - dqvc_con(i)
558
559          !--Add tendencies
560          cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
561          qcld(i)   = MIN(qtot(i), qcld(i) + dqt_con)
562          qvc(i)    = MIN(qcld(i), qvc(i) + dqvc_con(i))
563
564        ENDIF ! ok_warm_cloud, cf_cond .GT. eps
565       
566        !--We then calculate the part that is greater than qsat
567        !--and lower than gamma_cond * qsat, which is the ice supersaturated region
568        pdf_x = qsat(i) / qsatl(i) * 100.
569        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
570        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
571        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
572        issrfra(i) = EXP( - pdf_y ) * ( 1. - cldfra(i) )
573        qissr(i) = ( pdf_e3 * ( 1. - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
574
575        issrfra(i) = issrfra(i) - dcf_con(i)
576        qissr(i) = qissr(i) - dqvc_con(i) - dqi_con(i)
577
578      ENDIF ! ( 1. - cldfra(i) ) .GT. eps
579
580      !--Calculation of the subsaturated clear sky fraction and water
581      subfra(i) = 1. - cldfra(i) - issrfra(i)
582      qsub(i) = qtot(i) - qcld(i) - qissr(i)
583
584
585      !--------------------------------------
586      !--           CLOUD MIXING           --
587      !--------------------------------------
588      !--This process mixes the cloud with its surroundings: the subsaturated clear sky,
589      !--and the supersaturated clear sky. It is activated if the cloud is big enough,
590      !--but does not cover the entire mesh.
591      !
592      IF ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) THEN
593        !--Initialisation
594        dcf_mix_sub   = 0.
595        dqt_mix_sub   = 0.
596        dqvc_mix_sub  = 0.
597        dcf_mix_issr  = 0.
598        dqt_mix_issr  = 0.
599        dqvc_mix_issr = 0.
600
601        !-- PART 1 - TURBULENT DIFFUSION
602
603        !--Clouds within the mesh are assumed to be ellipses. The length of the
604        !--semi-major axis is a and the length of the semi-minor axis is b.
605        !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer.
606        !--In particular, it is 0 if cldfra = 1.
607        !--clouds_perim is the total perimeter of the clouds within the mesh,
608        !--not considering interfaces with other meshes (only the interfaces with clear
609        !--sky are taken into account).
610        !--
611        !--The area of each cloud is A = a * b * RPI,
612        !--and the perimeter of each cloud is
613        !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
614        !--
615        !--With cell_area the area of the cell, we have:
616        !-- cldfra = A * N_cld_mix / cell_area
617        !-- clouds_perim = P * N_cld_mix
618        !--
619        !--We assume that the ratio between b and a is a function of
620        !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
621        !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
622        !--are spherical.
623        !-- b / a = bovera = MAX(0.1, cldfra)
624        bovera = MAX(0.1, cldfra(i))
625        !--P / a is a function of b / a only, that we can calculate
626        !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
627        Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
628        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
629        !--based on observations.
630        !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
631        !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
632        !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
633        !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
634        a_mix = Povera / coef_mixing_lscp / RPI / ( 1. - cldfra(i) ) / bovera
635        !--and finally,
636        !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
637        N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) ) * cell_area(i) &
638                  / Povera / a_mix
639
640        !--The time required for turbulent diffusion to homogenize a region of size
641        !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
642        !--We compute L_mix and assume that the cloud is mixed over this length
643        L_mix = SQRT( dtime**3 * pbl_eps(i) )
644        !--The mixing length cannot be greater than the semi-minor axis. In this case,
645        !--the entire cloud is mixed.
646        L_mix = MIN(L_mix, a_mix * bovera)
647
648        !--The fraction of clear sky mixed is
649        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
650        envfra_mix = N_cld_mix * RPI / cell_area(i) &
651                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
652        !--The fraction of cloudy sky mixed is
653        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
654        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
655                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
656
657
658        !-- PART 2 - SHEARING
659
660        !--The clouds are then sheared. We keep the shape and number
661        !--assumptions from before. The clouds are sheared along their
662        !--semi-major axis (a_mix), on the entire cell heigh dz.
663        !--The increase in size is
664        L_shear = coef_shear_lscp * shear(i) * dz * dtime
665        !--therefore, the fraction of clear sky mixed is
666        !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
667        !--and the fraction of cloud mixed is
668        !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area
669        !--(note that they are equal)
670        shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i)
671        !--and the environment and cloud mixed fractions are the same,
672        !--which we add to the previous calculated mixed fractions.
673        !--We therefore assume that the sheared clouds and the turbulent
674        !--mixed clouds are different.
675        envfra_mix = envfra_mix + shear_fra
676        cldfra_mix = cldfra_mix + shear_fra
677
678
679        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
680
681        !--The environment fraction is allocated to subsaturated sky or supersaturated sky,
682        !--according to the factor sigma_mix. This is computed as the ratio of the
683        !--subsaturated sky fraction to the environment fraction, corrected by a factor
684        !--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
685        !--supersaturated sky is favoured. Physically, this means that it is more likely
686        !--to have supersaturated sky around the cloud than subsaturated sky.
687        sigma_mix = subfra(i) / ( subfra(i) + chi_mixing_lscp * issrfra(i) )
688        subfra_mix = MIN( sigma_mix * envfra_mix, subfra(i) )
689        issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra(i) )
690        cldfra_mix = MIN( cldfra_mix, cldfra(i) )
691
692        !--First, we mix the subsaturated sky (subfra_mix) and the cloud close
693        !--to this fraction (sigma_mix * cldfra_mix).
694        IF ( subfra(i) .GT. eps ) THEN
695          !--We compute the total humidity in the mixed air, which
696          !--can be either sub- or supersaturated.
697          qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
698                    + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
699                    / ( subfra_mix + cldfra_mix * sigma_mix )
700
701          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
702            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * sigma_mix / cldfra(i) &
703                      / ( subfra_mix + cldfra_mix * sigma_mix )
704            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
705                                        pplay(i), dtime, qvapincld_new)
706            IF ( qvapincld_new .EQ. 0. ) THEN
707              !--If all the ice has been sublimated, we sublimate
708              !--completely the cloud and do not activate the sublimation
709              !--process
710              !--Tendencies and diagnostics
711              dcf_mix_sub = - sigma_mix * cldfra_mix
712              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
713              dqvc_mix_sub = dcf_mix_sub * qvc(i) / cldfra(i)
714            ELSE
715              dcf_mix_sub = subfra_mix
716              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
717              dqvc_mix_sub = dcf_mix_sub * qvapincld_new
718            ENDIF
719          ELSE
720            IF ( qvapinmix .GT. qsat(i) ) THEN
721              !--If the mixed air is supersaturated, we condense the subsaturated
722              !--region which was mixed.
723              dcf_mix_sub = subfra_mix
724              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
725              dqvc_mix_sub = dcf_mix_sub * qsat(i)
726            ELSE
727              !--Else, we sublimate the cloud which was mixed.
728              dcf_mix_sub = - sigma_mix * cldfra_mix
729              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
730              dqvc_mix_sub = dcf_mix_sub * qsat(i)
731            ENDIF
732          ENDIF ! ok_unadjusted_clouds
733        ENDIF ! subfra .GT. eps
734   
735        !--We then mix the supersaturated sky (issrfra_mix) and the cloud,
736        !--for which the mixed air is always supersatured, therefore
737        !--the cloud necessarily expands
738        IF ( issrfra(i) .GT. eps ) THEN
739
740          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
741            qvapinmix = ( qissr(i) * issrfra_mix / issrfra(i) &
742                      + qcld(i) * cldfra_mix * (1. - sigma_mix) / cldfra(i) ) &
743                      / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
744            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * (1. - sigma_mix) / cldfra(i) &
745                      / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
746            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
747                                        pplay(i), dtime, qvapincld_new)
748            dcf_mix_issr = issrfra_mix
749            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
750            dqvc_mix_issr = dcf_mix_issr * qvapincld_new
751          ELSE
752            !--In this case, the additionnal vapor condenses
753            dcf_mix_issr = issrfra_mix
754            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
755            dqvc_mix_issr = dcf_mix_issr * qsat(i)
756          ENDIF ! ok_unadjusted_clouds
757
758        ENDIF ! issrfra .GT. eps
759
760        !--Sum up the tendencies from subsaturated sky and supersaturated sky
761        dcf_mix(i)  = dcf_mix_sub  + dcf_mix_issr
762        dqt_mix     = dqt_mix_sub  + dqt_mix_issr
763        dqvc_mix(i) = dqvc_mix_sub + dqvc_mix_issr
764        dqi_mix(i)  = dqt_mix - dqvc_mix(i)
765
766        IF ( ok_plane_contrail .AND. ( rcont_seri(i) .GT. eps ) ) THEN
767          CALL contrails_mixing( &
768              dtime, pplay(i), shear(i), pbl_eps(i), cell_area(i), dz, &
769              temp(i), qtot(i), qsat(i), subfra(i), qsub(i), issrfra(i), qissr(i), &
770              cldfra(i), qcld(i), qvc(i), rcont_seri(i), &
771              dcf_mix(i), dqvc_mix(i), dqi_mix(i), dqt_mix, dcf_mix_issr, dqt_mix_issr)
772        ENDIF
773
774        !--Add tendencies
775        issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
776        qissr(i)   = MAX(0., qissr(i) - dqt_mix_issr)
777        cldfra(i)  = MAX(0., MIN(1., cldfra(i) + dcf_mix(i)))
778        qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqt_mix))
779        qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))
780
781      ENDIF ! ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
782
783
784      !----------------------------------------
785      !--       CONTRAILS AND AVIATION       --
786      !----------------------------------------
787
788      IF ( ok_plane_contrail ) THEN
789        CALL contrails_formation_evolution( &
790            dtime, pplay(i), temp(i), qsat(i), qsatl(i), gamma_cond(i), &
791            rcont_seri(i), flight_dist(i), cldfra(i), qvc(i), &
792            V_cell, M_cell, pdf_loc, pdf_scale, pdf_alpha, &
793            Tcritcont(i), qcritcont(i), potcontfraP(i), potcontfraNP(i), contfra(i), &
794            dcf_avi(i), dqvc_avi(i), dqi_avi(i) &
795            )
796
797        !--Add tendencies
798        issrfra(i) = MAX(0., issrfra(i) - dcf_avi(i))
799        qissr(i)   = MAX(0., qissr(i) - dqvc_avi(i) - dqi_avi(i))
800        cldfra(i)  = MIN(1., cldfra(i) + dcf_avi(i))
801        qcld(i)    = MIN(qtot(i), qcld(i) + dqvc_avi(i) + dqi_avi(i))
802        qvc(i)     = MIN(qcld(i), qvc(i) + dqvc_avi(i))
803      ENDIF
804
805      !-------------------------------------------
806      !--       FINAL BARRIERS AND OUTPUTS      --
807      !-------------------------------------------
808
809      IF ( cldfra(i) .LT. eps ) THEN
810        !--If the cloud is too small, it is sublimated.
811        cldfra(i) = 0.
812        qcld(i)   = 0.
813        qvc(i)    = 0.
814        qincld(i) = qsat(i)
815      ELSE
816        qincld(i) = qcld(i) / cldfra(i)
817      ENDIF ! cldfra .LT. eps
818
819      !--Diagnostics
820      dcf_sub(i)  = dcf_sub(i)  / dtime
821      dcf_con(i)  = dcf_con(i)  / dtime
822      dcf_mix(i)  = dcf_mix(i)  / dtime
823      dqi_adj(i)  = dqi_adj(i)  / dtime
824      dqi_sub(i)  = dqi_sub(i)  / dtime
825      dqi_con(i)  = dqi_con(i)  / dtime
826      dqi_mix(i)  = dqi_mix(i)  / dtime
827      dqvc_adj(i) = dqvc_adj(i) / dtime
828      dqvc_sub(i) = dqvc_sub(i) / dtime
829      dqvc_con(i) = dqvc_con(i) / dtime
830      dqvc_mix(i) = dqvc_mix(i) / dtime
831      IF ( ok_plane_contrail ) THEN
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.