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

Last change on this file since 5453 was 5453, checked in by aborella, 37 hours ago

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