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

Last change on this file since 5748 was 5614, checked in by evignon, 4 months ago

Commission liée à un update majeur de la routine de condensation grande echelle suite au travail
de Lea, Audran et Etienne
Elle inclue une restructuration des routines pour clarifier le role "moniteur" de la routine lscp_main,
une mise à jour de la parametrisation de partitionnement de phase de Lea pour inclure les nuages de couche limite,
ainsi que des corrections des routines de precipitations "poprecip".
Convergence numerique verifiee en prod et debug pour les physiques NPv6.3 et 7.0.1c
File size: 113.5 KB
Line 
1MODULE lmdz_lscp_condensation
2!----------------------------------------------------------------
3! Module for condensation of clouds routines
4! that are called in LSCP
5
6
7IMPLICIT NONE
8
9CONTAINS
10
11!**********************************************************************************
12SUBROUTINE condensation_lognormal( &
13      klon, temp, qtot, qsat, gamma_cond, ratqs, &
14      keepgoing, cldfra, qincld, qvc)
15
16!----------------------------------------------------------------------
17! This subroutine calculates the formation of clouds, using a
18! statistical cloud scheme. It uses a generalised lognormal distribution
19! of total water in the gridbox
20! See Bony and Emanuel, 2001
21!----------------------------------------------------------------------
22
23USE lmdz_lscp_ini, ONLY: eps
24
25IMPLICIT NONE
26
27!
28! Input
29!
30INTEGER,  INTENT(IN)                     :: klon          ! number of horizontal grid points
31!
32REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
33REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot          ! total specific humidity (without precip) [kg/kg]
34REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
35REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
36REAL,     INTENT(IN)   , DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
37LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
38!
39!  Output
40!  NB. those are in INOUT because of the convergence loop on temperature
41!      (in some cases, the values are not re-computed) but the values
42!      are never used explicitely
43!
44REAL,     INTENT(INOUT), DIMENSION(klon) :: cldfra   ! cloud fraction [-]
45REAL,     INTENT(INOUT), DIMENSION(klon) :: qincld   ! cloud-mean in-cloud total specific water [kg/kg]
46REAL,     INTENT(INOUT), DIMENSION(klon) :: qvc      ! gridbox-mean vapor in the cloud [kg/kg]
47!
48! Local
49!
50INTEGER :: i
51REAL :: pdf_std, pdf_k, pdf_delta
52REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
53!
54!--Loop on klon
55!
56DO i = 1, klon
57
58  !--If a new calculation of the condensation is needed,
59  !--i.e., temperature has not yet converged (or the cloud is
60  !--formed elsewhere)
61  IF (keepgoing(i)) THEN
62
63    pdf_std   = ratqs(i) * qtot(i)
64    pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) )
65    pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
66    pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
67    pdf_b     = pdf_k / (2. * SQRT(2.))
68    pdf_e1    = pdf_a - pdf_b
69    pdf_e1    = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
70    pdf_e1    = 1. - ERF(pdf_e1)
71    pdf_e2    = pdf_a + pdf_b
72    pdf_e2    = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
73    pdf_e2    = 1. - ERF(pdf_e2)
74
75
76    IF ( pdf_e1 .LT. eps ) THEN
77      cldfra(i) = 0.
78      qincld(i) = qsat(i)
79      !--AB grid-mean vapor in the cloud - we assume saturation adjustment
80      qvc(i) = 0.
81    ELSE
82      cldfra(i) = 0.5 * pdf_e1
83      qincld(i) = qtot(i) * pdf_e2 / pdf_e1
84      !--AB grid-mean vapor in the cloud - we assume saturation adjustment
85      qvc(i) = qsat(i) * cldfra(i)
86    ENDIF
87
88  ENDIF   ! end keepgoing
89
90ENDDO     ! end loop on i
91END SUBROUTINE condensation_lognormal
92!**********************************************************************************
93
94!**********************************************************************************
95SUBROUTINE condensation_ice_supersat( &
96      klon, dtime, missing_val, pplay, paprsdn, paprsup, &
97      cf_seri, rvc_seri, ql_seri, qi_seri, shear, pbl_eps, cell_area, &
98      temp, qtot, qsat, gamma_cond, ratqs, keepgoing, &
99      cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, &
100      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
101      Tcontr, qcontr, qcontr2, fcontrN, fcontrP, flight_dist, flight_h2o, &
102      dcf_avi, dqi_avi, dqvc_avi)
103
104!----------------------------------------------------------------------
105! This subroutine calculates the formation, evolution and dissipation
106! of clouds, using a process-oriented treatment of the cloud properties
107! (cloud fraction, vapor in the cloud, condensed water in the cloud).
108! It allows for ice supersaturation in cold regions, in clear sky.
109! If ok_unadjusted_clouds, it also allows for sub- and supersaturated
110! cloud water vapors.
111! It also allows for the formation and evolution of condensation trails
112! (contrails) from aviation.
113! Authors: Audran Borella, Etienne Vignon, Olivier Boucher
114! April 2024
115!----------------------------------------------------------------------
116
117USE lmdz_lscp_tools,      ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC
118USE lmdz_lscp_ini,        ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
119USE lmdz_lscp_ini,        ONLY: eps, temp_nowater, ok_weibull_warm_clouds
120USE lmdz_lscp_ini,        ONLY: ok_unadjusted_clouds, iflag_cloud_sublim_pdf
121USE lmdz_lscp_ini,        ONLY: lunout
122
123USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp, &
124                                mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
125                                std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, &
126                                coef_mixing_lscp, coef_shear_lscp, &
127                                chi_mixing_lscp, rho_ice
128
129IMPLICIT NONE
130
131!
132! Input
133!
134INTEGER,  INTENT(IN)                     :: klon          ! number of horizontal grid points
135REAL,     INTENT(IN)                     :: dtime         ! time step [s]
136REAL,     INTENT(IN)                     :: missing_val   ! missing value for output
137!
138REAL,     INTENT(IN)   , DIMENSION(klon) :: pplay         ! layer pressure [Pa]
139REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsdn       ! pressure at the lower interface [Pa]
140REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsup       ! pressure at the upper interface [Pa]
141REAL,     INTENT(IN)   , DIMENSION(klon) :: cf_seri       ! cloud fraction [-]
142REAL,     INTENT(IN)   , DIMENSION(klon) :: rvc_seri      ! gridbox-mean water vapor in cloud [kg/kg]
143REAL,     INTENT(IN)   , DIMENSION(klon) :: ql_seri       ! specific liquid water content [kg/kg]
144REAL,     INTENT(IN)   , DIMENSION(klon) :: qi_seri       ! specific ice water content [kg/kg]
145REAL,     INTENT(IN)   , DIMENSION(klon) :: shear         ! vertical shear [s-1]
146REAL,     INTENT(IN)   , DIMENSION(klon) :: pbl_eps       !
147REAL,     INTENT(IN)   , DIMENSION(klon) :: cell_area     !
148REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
149REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot          ! total specific humidity (without precip) [kg/kg]
150REAL,     INTENT(IN)   , DIMENSION(klon) :: qsat          ! saturation specific humidity [kg/kg]
151REAL,     INTENT(IN)   , DIMENSION(klon) :: gamma_cond    ! condensation threshold w.r.t. qsat [-]
152REAL,     INTENT(IN)   , DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
153LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
154!
155! Input for aviation
156!
157REAL,     INTENT(IN),    DIMENSION(klon) :: flight_dist   !
158REAL,     INTENT(IN),    DIMENSION(klon) :: flight_h2o    !
159!
160!  Output
161!  NB. cldfra and qincld should be outputed as cf_seri and qi_seri,
162!      or as tendencies (maybe in the future)
163!  NB. those are in INOUT because of the convergence loop on temperature
164!      (in some cases, the values are not re-computed) but the values
165!      are never used explicitely
166!
167REAL,     INTENT(INOUT), DIMENSION(klon) :: cldfra   ! cloud fraction [-]
168REAL,     INTENT(INOUT), DIMENSION(klon) :: qincld   ! cloud-mean in-cloud total specific water [kg/kg]
169REAL,     INTENT(INOUT), DIMENSION(klon) :: qvc      ! gridbox-mean vapor in the cloud [kg/kg]
170REAL,     INTENT(INOUT), DIMENSION(klon) :: issrfra  ! ISSR fraction [-]
171REAL,     INTENT(INOUT), DIMENSION(klon) :: qissr    ! gridbox-mean ISSR specific water [kg/kg]
172!
173!  Diagnostics for condensation and ice supersaturation
174!  NB. idem for the INOUT
175!
176REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_sub  ! cloud fraction tendency because of sublimation [s-1]
177REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_con  ! cloud fraction tendency because of condensation [s-1]
178REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_mix  ! cloud fraction tendency because of cloud mixing [s-1]
179REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_adj  ! specific ice content tendency because of temperature adjustment [kg/kg/s]
180REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_sub  ! specific ice content tendency because of sublimation [kg/kg/s]
181REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_con  ! specific ice content tendency because of condensation [kg/kg/s]
182REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_mix  ! specific ice content tendency because of cloud mixing [kg/kg/s]
183REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
184REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s]
185REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s]
186REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
187!
188!  Diagnostics for aviation
189!  NB. idem for the INOUT
190!
191REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcontr   ! critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) [K]
192REAL,     INTENT(INOUT), DIMENSION(klon) :: qcontr   !
193REAL,     INTENT(INOUT), DIMENSION(klon) :: qcontr2  !
194REAL,     INTENT(INOUT), DIMENSION(klon) :: fcontrN  !
195REAL,     INTENT(INOUT), DIMENSION(klon) :: fcontrP  !
196REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_avi  ! cloud fraction tendency because of aviation [s-1]
197REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_avi  ! specific ice content tendency because of aviation [kg/kg/s]
198REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_avi ! specific cloud water vapor tendency because of aviation [kg/kg/s]
199!
200! Local
201!
202INTEGER :: i
203LOGICAL :: ok_warm_cloud
204REAL, DIMENSION(klon) :: qcld, qzero
205!
206! for lognormal
207REAL :: pdf_std, pdf_k, pdf_delta
208REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
209!
210! for unadjusted clouds
211REAL :: qvapincld, qvapincld_new
212REAL :: qiceincld
213!
214! for sublimation
215REAL :: pdf_alpha
216REAL :: dqt_sub
217!
218! for condensation
219REAL, DIMENSION(klon) :: qsatl, dqsatl
220REAL :: clrfra, qclr, sl_clr, rhl_clr
221REAL :: pdf_ratqs, pdf_skew, pdf_scale, pdf_loc
222REAL :: pdf_x, pdf_y, pdf_T
223REAL :: pdf_e3, pdf_e4
224REAL :: cf_cond, qt_cond, dqt_con
225!
226! for mixing
227REAL, DIMENSION(klon) :: subfra, qsub
228REAL :: dqt_mix_sub, dqt_mix_issr
229REAL :: dcf_mix_sub, dcf_mix_issr
230REAL :: dqvc_mix_sub, dqvc_mix_issr
231REAL :: dqt_mix
232REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix
233REAL :: envfra_mix, cldfra_mix
234REAL :: L_shear, shear_fra
235REAL :: sigma_mix, issrfra_mix, subfra_mix, qvapinmix
236!
237! for cell properties
238REAL :: rho, rhodz, dz
239!REAL :: V_cell, M_cell
240!
241! for aviation and cell properties
242!REAL :: dqt_avi
243!REAL :: contrail_fra
244!
245!
246!--more local variables for diagnostics
247!--imported from YOMCST.h
248!--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1)
249!--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air
250!--values from Schumann, Meteorol Zeitschrift, 1996
251!--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
252!--Qheat = 43.  /  50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
253!REAL, PARAMETER :: EiH2O=1.25  !--emission index of water vapour for kerosene (kg kg-1)
254!REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1)
255!REAL, PARAMETER :: eta=0.3     !--average propulsion efficiency of the aircraft
256!--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
257!--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
258!REAL :: Gcontr
259!--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2)
260!--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1)
261!REAL :: qsatliqcontr
262
263
264!-----------------------------------------------
265! Initialisations
266!-----------------------------------------------
267
268! Ajout des émissions de H2O dues à l'aviation
269! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
270! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) 
271!      = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
272! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
273! flight_h2O is in kg H2O / s / cell
274!
275!IF (ok_plane_h2o) THEN
276!  q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) )
277!ENDIF
278
279
280qzero(:) = 0.
281
282!--Calculation of qsat w.r.t. liquid
283CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsatl)
284
285!
286!--Loop on klon
287!
288DO i = 1, klon
289
290  !--If a new calculation of the condensation is needed,
291  !--i.e., temperature has not yet converged (or the cloud is
292  !--formed elsewhere)
293  IF (keepgoing(i)) THEN
294
295    !--Initialisation
296    issrfra(i)  = 0.
297    qissr(i)    = 0.
298
299    !--If the temperature is higher than the threshold below which
300    !--there is no liquid in the gridbox, we activate the usual scheme
301    !--(generalised lognormal from Bony and Emanuel 2001)
302    !--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for
303    !--all clouds, and the lognormal scheme is not activated
304    IF ( ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ) THEN
305
306      pdf_std   = ratqs(i) * qtot(i)
307      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) )
308      pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
309      pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
310      pdf_b     = pdf_k / (2. * SQRT(2.))
311      pdf_e1    = pdf_a - pdf_b
312      pdf_e1    = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
313      pdf_e1    = 1. - ERF(pdf_e1)
314      pdf_e2    = pdf_a + pdf_b
315      pdf_e2    = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
316      pdf_e2    = 1. - ERF(pdf_e2)
317
318
319      IF ( pdf_e1 .LT. eps ) THEN
320        cldfra(i) = 0.
321        qincld(i) = qsat(i)
322        qvc(i) = 0.
323      ELSE
324        cldfra(i) = 0.5 * pdf_e1
325        qincld(i) = qtot(i) * pdf_e2 / pdf_e1
326        qvc(i) = qsat(i) * cldfra(i)
327      ENDIF
328
329    !--If the temperature is lower than temp_nowater, we use the new
330    !--condensation scheme that allows for ice supersaturation
331    ELSE
332
333      !--Initialisation
334      IF ( temp(i) .GT. temp_nowater ) THEN
335        !--If the air mass is warm (liquid water can exist),
336        !--all the memory is lost and the scheme becomes statistical,
337        !--i.e., the sublimation and mixing processes are deactivated,
338        !--and the condensation process is slightly adapted
339        !--This can happen only if ok_weibull_warm_clouds = .TRUE.
340        ! AB WARM CLOUD
341        !cldfra(i) = 0.
342        !qcld(i)   = 0.
343        !qvc(i)    = 0.
344        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
345        qcld(i)   = MAX(0., MIN(qtot(i), ql_seri(i) + qi_seri(i) + rvc_seri(i) * qtot(i)))
346        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
347        ok_warm_cloud = .TRUE.
348      ELSE
349        !--The following barriers ensure that the traced cloud properties
350        !--are consistent. In some rare cases, i.e. the cloud water vapor
351        !--can be greater than the total water in the gridbox
352        cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
353        qcld(i)   = MAX(0., MIN(qtot(i), rvc_seri(i) * qtot(i) + qi_seri(i)))
354        qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
355        ok_warm_cloud = .FALSE.
356      ENDIF
357
358      dcf_sub(i)  = 0.
359      dqi_sub(i)  = 0.
360      dqvc_sub(i) = 0.
361      dqi_adj(i)  = 0.
362      dqvc_adj(i) = 0.
363      dcf_con(i)  = 0.
364      dqi_con(i)  = 0.
365      dqvc_con(i) = 0.
366      dcf_mix(i)  = 0.
367      dqi_mix(i)  = 0.
368      dqvc_mix(i) = 0.
369
370      !--Initialisation of the cell properties
371      !--Dry density [kg/m3]
372      rho = pplay(i) / temp(i) / RD
373      !--Dry air mass [kg/m2]
374      rhodz = ( paprsdn(i) - paprsup(i) ) / RG
375      !--Cell thickness [m]
376      dz = rhodz / rho
377      !--Cell volume [m3]
378      !V_cell = dz * cell_area(i)
379      !--Cell dry air mass [kg]
380      !M_cell = rhodz * cell_area(i)
381
382
383      !-------------------------------------------------------------------
384      !--    SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD    --
385      !-------------------------------------------------------------------
386     
387      !--If there is a cloud
388      IF ( cldfra(i) .GT. eps ) THEN
389       
390        qvapincld = qvc(i) / cldfra(i)
391        qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
392       
393        !--If the ice water content is too low, the cloud is purely sublimated
394        !--Most probably, we advected a cloud with no ice water content (possible
395        !--if the entire cloud precipited for example)
396        IF ( qiceincld .LT. eps ) THEN
397          dcf_sub(i) = - cldfra(i)
398          dqvc_sub(i) = - qvc(i)
399          dqi_sub(i) = - ( qcld(i) - qvc(i) )
400
401          cldfra(i) = 0.
402          qcld(i) = 0.
403          qvc(i) = 0.
404
405        !--Else, the cloud is adjusted and sublimated
406        ELSE
407
408          !--The vapor in cloud cannot be higher than the
409          !--condensation threshold
410          qvapincld = MIN(qvapincld, gamma_cond(i) * qsat(i))
411          qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
412
413          ! AB WARM CLOUD
414          !IF ( ok_unadjusted_clouds ) THEN
415          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
416            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
417                                        pplay(i), dtime, qvapincld_new)
418            IF ( qvapincld_new .EQ. 0. ) THEN
419              !--If all the ice has been sublimated, we sublimate
420              !--completely the cloud and do not activate the sublimation
421              !--process
422              !--Tendencies and diagnostics
423              dcf_sub(i) = - cldfra(i)
424              dqvc_sub(i) = - qvc(i)
425              dqi_sub(i) = - ( qcld(i) - qvc(i) )
426
427              cldfra(i) = 0.
428              qcld(i) = 0.
429              qvc(i) = 0.
430            ENDIF
431          ELSE
432            !--We keep the saturation adjustment hypothesis, and the vapor in the
433            !--cloud is set equal to the saturation vapor
434            qvapincld_new = qsat(i)
435          ENDIF ! ok_unadjusted_clouds
436         
437          !--Adjustment of the IWC to the new vapor in cloud
438          !--(this can be either positive or negative)
439          dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
440          dqi_adj(i) = - dqvc_adj(i)
441
442          !--Add tendencies
443          !--The vapor in the cloud is updated, but not qcld as it is constant
444          !--through this process, as well as cldfra which is unmodified
445          qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
446         
447
448          !------------------------------------
449          !--    DISSIPATION OF THE CLOUD    --
450          !------------------------------------
451
452          !--If the vapor in cloud is below vapor needed for the cloud to survive
453          IF ( ( qvapincld .LT. qvapincld_new ) .OR. ( iflag_cloud_sublim_pdf .GE. 4 ) ) THEN
454            !--Sublimation of the subsaturated cloud
455            !--iflag_cloud_sublim_pdf selects the PDF of the ice water content
456            !--to use.
457            !--iflag = 1 --> uniform distribution
458            !--iflag = 2 --> exponential distribution
459            !--iflag = 3 --> gamma distribution (Karcher et al 2018)
460
461            IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN
462              !--Uniform distribution starting at qvapincld
463              pdf_e1 = 1. / ( 2. * qiceincld )
464
465              dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
466              dqt_sub    = - cldfra(i) * ( qvapincld_new**2 - qvapincld**2 ) &
467                                       * pdf_e1 / 2.
468
469            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN
470              !--Exponential distribution starting at qvapincld
471              pdf_alpha = 1. / qiceincld
472              pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) )
473
474              dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 )
475              dqt_sub    = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha &
476                                       + qvapincld - qvapincld_new * pdf_e1 )
477
478            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 3 ) THEN
479              !--Gamma distribution starting at qvapincld
480              pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
481              pdf_y = pdf_alpha * ( qvapincld_new - qvapincld )
482              pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
483              pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
484
485              dcf_sub(i) = - cldfra(i) * pdf_e1
486              dqt_sub    = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
487
488            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 4 ) THEN
489              !--Normal distribution around qvapincld + qiceincld with width
490              !--constant in the RHi space
491              pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
492                    / ( std_subl_pdf_lscp / 100. * qsat(i))
493              pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
494              pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )
495
496              dcf_sub(i) = - cldfra(i) * pdf_e1
497              dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
498                                         - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
499                                                                                       
500            ENDIF
501
502            !--Tendencies and diagnostics
503            dqvc_sub(i) = dqt_sub
504
505            !--Add tendencies
506            cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
507            qcld(i)   = MAX(0., qcld(i) + dqt_sub)
508            qvc(i)    = MAX(0., qvc(i) + dqvc_sub(i))
509
510          ENDIF ! qvapincld .LT. qvapincld_new
511
512        ENDIF   ! qiceincld .LT. eps
513      ENDIF ! cldfra(i) .GT. eps
514
515
516      !--------------------------------------------------------------------------
517      !--    CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS    --
518      !--------------------------------------------------------------------------
519      !--This section relies on a distribution of water in the clear-sky region of
520      !--the mesh.
521
522      !--If there is a clear-sky region
523      IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
524
525        !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
526        qclr = qtot(i) - qcld(i)
527
528        !--New PDF
529        rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.
530
531        !--Calculation of the properties of the PDF
532        !--Parameterization from IAGOS observations
533        !--pdf_e1 and pdf_e2 will be reused below
534
535        !--Coefficient for standard deviation:
536        !--  tuning coef * (clear sky area**0.25) * (function of temperature)
537        !pdf_e1 = beta_pdf_lscp &
538        !       * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
539        !       * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
540        !--  tuning coef * (cell area**0.25) * (function of temperature)
541        pdf_e1 = beta_pdf_lscp * ( ( 1. - cldfra(i) ) * cell_area(i) )**0.25 &
542                * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
543        IF ( rhl_clr .GT. 50. ) THEN
544          pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp
545        ELSE
546          pdf_std = pdf_e1 * rhl_clr / 50.
547        ENDIF
548        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
549        pdf_alpha = EXP( rhl_clr / 100. ) * pdf_e3
550        pdf_alpha = MIN(10., pdf_alpha)
551       
552        IF ( ok_warm_cloud ) THEN
553          !--If the statistical scheme is activated, the standard deviation is adapted
554          !--to depend on the pressure level. It is multiplied by ratqs, so that near the
555          !--surface std is almost 0, and upper than about 450 hPa the std is left untouched
556          pdf_std = pdf_std * ratqs(i)
557        ENDIF
558
559        pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
560        pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2 ))
561        pdf_loc = rhl_clr - pdf_scale * pdf_e2
562
563        !--Calculation of the newly condensed water and fraction (pronostic)
564        !--Integration of the clear sky PDF between gamma_cond*qsat and +inf
565        !--NB. the calculated values are clear-sky averaged
566
567        pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
568        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
569        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
570        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
571        cf_cond = EXP( - pdf_y )
572        qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.
573
574        ! AB WARM CLOUD
575        !IF ( ok_warm_cloud ) THEN
576        !  !--If the statistical scheme is activated, the calculated increase is equal
577        !  !--to the cloud fraction, we assume saturation adjustment, and the
578        !  !--condensation process stops
579        !  cldfra(i) = cf_cond
580        !  qcld(i) = qt_cond
581        !  qvc(i) = cldfra(i) * qsat(i)
582
583        !ELSEIF ( cf_cond .GT. eps ) THEN
584        IF ( cf_cond .GT. eps ) THEN
585
586          dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond
587          dqt_con = ( 1. - cldfra(i) ) * qt_cond
588         
589          !--Barriers
590          dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i))
591          dqt_con = MIN(dqt_con, qclr)
592
593
594          ! AB WARM CLOUD
595          !IF ( ok_unadjusted_clouds ) THEN
596          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
597            !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
598            !--the new vapor qvapincld. The timestep is divided by two because we do not
599            !--know when the condensation occurs
600            qvapincld = gamma_cond(i) * qsat(i)
601            qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
602            CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
603                                        pplay(i), dtime / 2., qvapincld_new)
604            qvapincld = qvapincld_new
605          ELSE
606            !--We keep the saturation adjustment hypothesis, and the vapor in the
607            !--newly formed cloud is set equal to the saturation vapor.
608            qvapincld = qsat(i)
609          ENDIF
610
611          !--Tendency on cloud vapor and diagnostic
612          dqvc_con(i) = qvapincld * dcf_con(i)
613          dqi_con(i) = dqt_con - dqvc_con(i)
614
615          !--Add tendencies
616          cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
617          qcld(i)   = MIN(qtot(i), qcld(i) + dqt_con)
618          qvc(i)    = MIN(qcld(i), qvc(i) + dqvc_con(i))
619
620        ENDIF ! ok_warm_cloud, cf_cond .GT. eps
621       
622        !--We then calculate the part that is greater than qsat
623        !--and lower than gamma_cond * qsat, which is the ice supersaturated region
624        pdf_x = qsat(i) / qsatl(i) * 100.
625        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
626        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
627        pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
628        issrfra(i) = EXP( - pdf_y ) * ( 1. - cldfra(i) )
629        qissr(i) = ( pdf_e3 * ( 1. - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
630
631        issrfra(i) = issrfra(i) - dcf_con(i)
632        qissr(i) = qissr(i) - dqvc_con(i) - dqi_con(i)
633
634      ENDIF ! ( 1. - cldfra(i) ) .GT. eps
635
636      !--Calculation of the subsaturated clear sky fraction and water
637      subfra(i) = 1. - cldfra(i) - issrfra(i)
638      qsub(i) = qtot(i) - qcld(i) - qissr(i)
639
640
641      !--------------------------------------
642      !--           CLOUD MIXING           --
643      !--------------------------------------
644      !--This process mixes the cloud with its surroundings: the subsaturated clear sky,
645      !--and the supersaturated clear sky. It is activated if the cloud is big enough,
646      !--but does not cover the entire mesh.
647      !
648      ! AB WARM CLOUD
649      !IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
650      !        .AND. .NOT. ok_warm_cloud ) THEN
651      IF ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) THEN
652
653        !--Initialisation
654        dcf_mix_sub   = 0.
655        dqt_mix_sub   = 0.
656        dqvc_mix_sub  = 0.
657        dcf_mix_issr  = 0.
658        dqt_mix_issr  = 0.
659        dqvc_mix_issr = 0.
660
661
662        !-- PART 1 - TURBULENT DIFFUSION
663
664        !--Clouds within the mesh are assumed to be ellipses. The length of the
665        !--semi-major axis is a and the length of the semi-minor axis is b.
666        !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer.
667        !--In particular, it is 0 if cldfra = 1.
668        !--clouds_perim is the total perimeter of the clouds within the mesh,
669        !--not considering interfaces with other meshes (only the interfaces with clear
670        !--sky are taken into account).
671        !--
672        !--The area of each cloud is A = a * b * RPI,
673        !--and the perimeter of each cloud is
674        !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
675        !--
676        !--With cell_area the area of the cell, we have:
677        !-- cldfra = A * N_cld_mix / cell_area
678        !-- clouds_perim = P * N_cld_mix
679        !--
680        !--We assume that the ratio between b and a is a function of
681        !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
682        !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
683        !--are spherical.
684        !-- b / a = bovera = MAX(0.1, cldfra)
685        bovera = MAX(0.1, cldfra(i))
686        !--P / a is a function of b / a only, that we can calculate
687        !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
688        Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
689        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
690        !--based on observations.
691        !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
692        !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
693        !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
694        !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
695        a_mix = Povera / coef_mixing_lscp / RPI / ( 1. - cldfra(i) ) / bovera
696        !--and finally,
697        !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
698        N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) ) * cell_area(i) &
699                  / Povera / a_mix
700
701        !--The time required for turbulent diffusion to homogenize a region of size
702        !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
703        !--We compute L_mix and assume that the cloud is mixed over this length
704        L_mix = SQRT( dtime**3 * pbl_eps(i) )
705        !--The mixing length cannot be greater than the semi-minor axis. In this case,
706        !--the entire cloud is mixed.
707        L_mix = MIN(L_mix, a_mix * bovera)
708
709        !--The fraction of clear sky mixed is
710        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
711        envfra_mix = N_cld_mix * RPI / cell_area(i) &
712                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
713        !--The fraction of cloudy sky mixed is
714        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
715        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
716                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
717
718
719        !-- PART 2 - SHEARING
720
721        !--The clouds are then sheared. We keep the shape and number
722        !--assumptions from before. The clouds are sheared along their
723        !--semi-major axis (a_mix), on the entire cell heigh dz.
724        !--The increase in size is
725        L_shear = coef_shear_lscp * shear(i) * dz * dtime
726        !--therefore, the fraction of clear sky mixed is
727        !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
728        !--and the fraction of cloud mixed is
729        !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area
730        !--(note that they are equal)
731        shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i)
732        !--and the environment and cloud mixed fractions are the same,
733        !--which we add to the previous calculated mixed fractions.
734        !--We therefore assume that the sheared clouds and the turbulent
735        !--mixed clouds are different.
736        envfra_mix = envfra_mix + shear_fra
737        cldfra_mix = cldfra_mix + shear_fra
738
739
740        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
741
742        !--The environment fraction is allocated to subsaturated sky or supersaturated sky,
743        !--according to the factor sigma_mix. This is computed as the ratio of the
744        !--subsaturated sky fraction to the environment fraction, corrected by a factor
745        !--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
746        !--supersaturated sky is favoured. Physically, this means that it is more likely
747        !--to have supersaturated sky around the cloud than subsaturated sky.
748        sigma_mix = subfra(i) / ( subfra(i) + chi_mixing_lscp * issrfra(i) )
749        subfra_mix = MIN( sigma_mix * envfra_mix, subfra(i) )
750        issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra(i) )
751        cldfra_mix = MIN( cldfra_mix, cldfra(i) )
752
753        !--First, we mix the subsaturated sky (subfra_mix) and the cloud close
754        !--to this fraction (sigma_mix * cldfra_mix).
755        IF ( subfra(i) .GT. eps ) THEN
756          !--We compute the total humidity in the mixed air, which
757          !--can be either sub- or supersaturated.
758          qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
759                    + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
760                    / ( subfra_mix + cldfra_mix * sigma_mix )
761
762          ! AB WARM CLOUD
763          !IF ( ok_unadjusted_clouds ) THEN
764          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
765            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * sigma_mix / cldfra(i) &
766                      / ( subfra_mix + cldfra_mix * sigma_mix )
767            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
768                                        pplay(i), dtime, qvapincld_new)
769            IF ( qvapincld_new .EQ. 0. ) THEN
770              !--If all the ice has been sublimated, we sublimate
771              !--completely the cloud and do not activate the sublimation
772              !--process
773              !--Tendencies and diagnostics
774              dcf_mix_sub = - sigma_mix * cldfra_mix
775              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
776              dqvc_mix_sub = dcf_mix_sub * qvc(i) / cldfra(i)
777            ELSE
778              dcf_mix_sub = subfra_mix
779              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
780              dqvc_mix_sub = dcf_mix_sub * qvapincld_new
781            ENDIF
782          ELSE
783            IF ( qvapinmix .GT. qsat(i) ) THEN
784              !--If the mixed air is supersaturated, we condense the subsaturated
785              !--region which was mixed.
786              dcf_mix_sub = subfra_mix
787              dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
788              dqvc_mix_sub = dcf_mix_sub * qsat(i)
789            ELSE
790              !--Else, we sublimate the cloud which was mixed.
791              dcf_mix_sub = - sigma_mix * cldfra_mix
792              dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
793              dqvc_mix_sub = dcf_mix_sub * qsat(i)
794            ENDIF
795          ENDIF ! ok_unadjusted_clouds
796        ENDIF ! subfra .GT. eps
797   
798        !--We then mix the supersaturated sky (issrfra_mix) and the cloud,
799        !--for which the mixed air is always supersatured, therefore
800        !--the cloud necessarily expands
801        IF ( issrfra(i) .GT. eps ) THEN
802
803          ! AB WARM CLOUD
804          !IF ( ok_unadjusted_clouds ) THEN
805          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
806            qvapinmix = ( qissr(i) * issrfra_mix / issrfra(i) &
807                      + qcld(i) * cldfra_mix * (1. - sigma_mix) / cldfra(i) ) &
808                      / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
809            qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * (1. - sigma_mix) / cldfra(i) &
810                      / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
811            CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
812                                        pplay(i), dtime, qvapincld_new)
813            dcf_mix_issr = issrfra_mix
814            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
815            dqvc_mix_issr = dcf_mix_issr * qvapincld_new
816          ELSE
817            !--In this case, the additionnal vapor condenses
818            dcf_mix_issr = issrfra_mix
819            dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
820            dqvc_mix_issr = dcf_mix_issr * qsat(i)
821          ENDIF ! ok_unadjusted_clouds
822
823
824        ENDIF ! issrfra .GT. eps
825
826        !--Sum up the tendencies from subsaturated sky and supersaturated sky
827        dcf_mix(i)  = dcf_mix_sub  + dcf_mix_issr
828        dqt_mix     = dqt_mix_sub  + dqt_mix_issr
829        dqvc_mix(i) = dqvc_mix_sub + dqvc_mix_issr
830        dqi_mix(i)  = dqt_mix - dqvc_mix(i)
831
832        !--Add tendencies
833        issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
834        qissr(i)   = MAX(0., qissr(i) - dqt_mix_issr)
835        cldfra(i)  = MAX(0., MIN(1., cldfra(i) + dcf_mix(i)))
836        qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqt_mix))
837        qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))
838
839      ENDIF ! ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
840
841
842      !----------------------------------------
843      !--       CONTRAILS AND AVIATION       --
844      !----------------------------------------
845
846      !--Add a source of cirrus from aviation contrails
847      !IF ( ok_plane_contrail ) THEN
848      !  dcf_avi(i) = 0.
849      !  dqi_avi(i) = 0.
850      !  dqvc_avi(i) = 0.
851      !  ! TODO implement ok_unadjusted_clouds
852      !  IF ( issrfra(i) .GT. eps ) THEN
853      !    contrail_fra = MIN(1., flight_m(i,k) * dtime * contrail_cross_section / V_cell)
854      !    dcf_avi(i) = issrfra(i) * contrail_fra
855      !    dqt_avi = dcf_avi(i) * qissr(i) / issrfra(i)
856      !    dqvc_avi(i) = qsat(i) * dcf_avi(i)
857      !   
858      !    !--Add tendencies
859      !    cldfra(i) = cldfra(i) + dcf_avi(i)
860      !    issrfra(i) = issrfra(i) - dcf_avi(i)
861      !    qcld(i) = qcld(i) + dqt_avi
862      !    qvc(i) = qvc(i) + dqvc_avi(i)
863      !    qissr(i) = qissr(i) - dqt_avi
864
865      !    !--Diagnostics
866      !    dqi_avi(i) = dqt_avi - qsat(i) * dcf_avi(i)
867      !  ENDIF
868      !  dcf_avi(i)  = dcf_avi(i)  / dtime
869      !  dqi_avi(i)  = dqi_avi(i)  / dtime
870      !  dqvc_avi(i) = dqvc_avi(i) / dtime
871      !ENDIF
872
873
874
875      !-------------------------------------------
876      !--       FINAL BARRIERS AND OUTPUTS      --
877      !-------------------------------------------
878
879      IF ( cldfra(i) .LT. eps ) THEN
880        !--If the cloud is too small, it is sublimated.
881        cldfra(i) = 0.
882        qcld(i)   = 0.
883        qvc(i)    = 0.
884        qincld(i) = qsat(i)
885      ELSE
886        qincld(i) = qcld(i) / cldfra(i)
887      ENDIF ! cldfra .LT. eps
888
889      !--Diagnostics
890      dcf_sub(i)  = dcf_sub(i)  / dtime
891      dcf_con(i)  = dcf_con(i)  / dtime
892      dcf_mix(i)  = dcf_mix(i)  / dtime
893      dqi_adj(i)  = dqi_adj(i)  / dtime
894      dqi_sub(i)  = dqi_sub(i)  / dtime
895      dqi_con(i)  = dqi_con(i)  / dtime
896      dqi_mix(i)  = dqi_mix(i)  / dtime
897      dqvc_adj(i) = dqvc_adj(i) / dtime
898      dqvc_sub(i) = dqvc_sub(i) / dtime
899      dqvc_con(i) = dqvc_con(i) / dtime
900      dqvc_mix(i) = dqvc_mix(i) / dtime
901
902    ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
903
904  ENDIF   ! end keepgoing
905
906ENDDO     ! end loop on i
907
908END SUBROUTINE condensation_ice_supersat
909!**********************************************************************************
910
911!**********************************************************************************
912SUBROUTINE deposition_sublimation( &
913    qvapincld, qiceincld, temp, qsat, pplay, dtime, &
914    qvapincld_new)
915
916USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, RD, EPS_W
917USE lmdz_lscp_ini,        ONLY: eps
918USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
919
920REAL,     INTENT(IN) :: qvapincld     !
921REAL,     INTENT(IN) :: qiceincld     !
922REAL,     INTENT(IN) :: temp          !
923REAL,     INTENT(IN) :: qsat          !
924REAL,     INTENT(IN) :: pplay         !
925REAL,     INTENT(IN) :: dtime         ! time step [s]
926REAL,     INTENT(OUT):: qvapincld_new !
927
928!
929! for unadjusted clouds
930REAL :: qice_ratio
931REAL :: pres_sat, rho, kappa
932REAL :: air_thermal_conduct, water_vapor_diff
933REAL :: iwc
934REAL :: iwc_log_inf100, iwc_log_sup100
935REAL :: iwc_inf100, alpha_inf100, coef_inf100
936REAL :: mu_sup100, sigma_sup100, coef_sup100
937REAL :: Dm_ice, rm_ice
938
939!--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
940!--hypothesis is lost, and the vapor in the cloud is purely prognostic.
941!
942!--The deposition equation is
943!-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
944!--from Lohmann et al. (2016), where
945!--alpha is the deposition coefficient [-]
946!--mi is the mass of one ice crystal [kg]
947!--C is the capacitance of an ice crystal [m]
948!--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
949!--R_v is the specific gas constant for humid air [J/kg/K]
950!--T is the temperature [K]
951!--esi is the saturation pressure w.r.t. ice [Pa]
952!--Dv is the diffusivity of water vapor [m2/s]
953!--Ls is the specific latent heat of sublimation [J/kg/K]
954!--ka is the thermal conductivity of dry air [J/m/s/K]
955!
956!--alpha is a coefficient to take into account the fact that during deposition, a water
957!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
958!--coherent (with the same structure). It has no impact for sublimation.
959!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
960!--during deposition, and alpha = 1. during sublimation.
961!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
962!-- C = capa_cond_cirrus * rm_ice
963!
964!--We have qice = Nice * mi, where Nice is the ice crystal
965!--number concentration per kg of moist air
966!--HYPOTHESIS 1: the ice crystals are spherical, therefore
967!-- mi = 4/3 * pi * rm_ice**3 * rho_ice
968!--HYPOTHESIS 2: the ice crystals are monodisperse with the
969!--initial radius rm_ice_0.
970!--NB. this is notably different than the assumption
971!--of a distributed qice in the cloud made in the sublimation process;
972!--should it be consistent?
973!
974!--As the deposition process does not create new ice crystals,
975!--and because we assume a same rm_ice value for all crystals
976!--therefore the sublimation process does not destroy ice crystals
977!--(or, in a limit case, it destroys all ice crystals), then
978!--Nice is a constant during the sublimation/deposition process.
979!-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
980!
981!--The deposition equation then reads:
982!-- 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
983!-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
984!--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
985!--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
986!-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
987!--  *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 )
988!--and we have
989!-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
990!-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
991!--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
992!
993!--This system of equations can be resolved with an exact
994!--explicit numerical integration, having one variable resolved
995!--explicitly, the other exactly. The exactly resolved variable is
996!--the one decreasing, so it is qvc if the process is
997!--condensation, qi if it is sublimation.
998!
999!--kappa is computed as an initialisation constant, as it depends only
1000!--on temperature and other pre-computed values
1001pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
1002!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
1003air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
1004!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
1005water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
1006kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
1007      * ( RV * temp / water_vapor_diff / pres_sat &
1008        + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
1009!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
1010
1011!--Dry density [kg/m3]
1012rho = pplay / temp / RD
1013
1014!--Here, the initial vapor in the cloud is qvapincld, and we compute
1015!--the new vapor qvapincld_new
1016
1017!--rm_ice formula from McFarquhar and Heymsfield (1997)
1018iwc = qiceincld * rho * 1e3
1019iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
1020iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
1021iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
1022
1023alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
1024coef_inf100 = iwc_inf100 * alpha_inf100**3 / 120.
1025
1026mu_sup100 = ( 5.2 + 0.0013 * ( temp - RTT ) ) &
1027          + ( 0.026 - 1.2E-3 * ( temp - RTT ) ) * iwc_log_sup100
1028sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp - RTT ) ) &
1029             + ( 0.018 - 2.1E-4 * ( temp - RTT ) ) * iwc_log_sup100
1030coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2 )
1031
1032Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2 ) &
1033       * coef_sup100 ) / ( coef_inf100 + coef_sup100 )
1034rm_ice = Dm_ice / 2. * 1.E-6
1035
1036IF ( qvapincld .GE. qsat ) THEN
1037  !--If the cloud is initially supersaturated
1038  !--Exact explicit integration (qvc exact, qice explicit)
1039  qvapincld_new = qsat + ( qvapincld - qsat ) &
1040                * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
1041ELSE
1042  !--If the cloud is initially subsaturated
1043  !--Exact explicit integration (qice exact, qvc explicit)
1044  !--The barrier is set so that the resulting vapor in cloud
1045  !--cannot be greater than qsat
1046  !--qice_ratio is the ratio between the new ice content and
1047  !--the old one, it is comprised between 0 and 1
1048  qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) )
1049
1050  IF ( qice_ratio .LT. 0. ) THEN
1051    !--The new vapor in cloud is set to 0 so that the
1052    !--sublimation process does not activate
1053    qvapincld_new = 0.
1054  ELSE
1055    !--Else, the sublimation process is activated with the
1056    !--diagnosed new cloud water vapor
1057    !--The new vapor in the cloud is increased with the
1058    !--sublimated ice
1059    qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 )
1060    !--The new vapor in the cloud cannot be greater than qsat
1061    qvapincld_new = MIN(qvapincld_new, qsat)
1062  ENDIF ! qice_ratio .LT. 0.
1063ENDIF ! qvapincld .GT. qsat
1064
1065END SUBROUTINE deposition_sublimation
1066
1067
1068!**********************************************************************************
1069
1070!**********************************************************************************
1071SUBROUTINE  condensation_cloudth(klon,                     &
1072&           temp,qt,qt_th,frac_th,zpspsk,play,thetal_th,   &
1073&           ratqs,sigma_qtherm,qsth,qsenv,qcloud,ctot,ctotth,ctot_vol,  &
1074&           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1075! This routine computes the condensation of clouds in convective boundary layers
1076! with thermals assuming two separate distribution of the saturation deficit in
1077! the thermal plumes and in the environment
1078! It is based on the work of Arnaud Jam (Jam et al. 2013, BLM)
1079! Author : Etienne Vignon (LMDZ/CNRS)
1080! Date: February 2025
1081! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
1082! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
1083!-----------------------------------------------------------------------------------
1084
1085      use lmdz_lscp_ini,    only: iflag_cloudth_vert,iflag_ratqs,iflag_cloudth_vert_noratqs
1086      use lmdz_lscp_ini,    only: vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin
1087      use lmdz_lscp_ini,    only: RTT, RG, RPI, RD, RV, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th
1088
1089      IMPLICIT NONE
1090
1091
1092!------------------------------------------------------------------------------
1093! Declarations
1094!------------------------------------------------------------------------------
1095
1096! INPUT/OUTPUT
1097
1098      INTEGER, INTENT(IN)                         :: klon
1099     
1100
1101      REAL, DIMENSION(klon),      INTENT(IN)      ::  temp          ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip
1102      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt            ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip
1103      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt_th         ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip
1104      REAL, DIMENSION(klon),      INTENT(IN)      ::  thetal_th     ! Liquid potential temperature in thermals [K]: has not seen the evap of precip
1105      REAL, DIMENSION(klon),      INTENT(IN)      ::  frac_th       ! Fraction of the mesh covered by thermals [0-1]
1106      REAL, DIMENSION(klon),      INTENT(IN)      ::  zpspsk        ! Exner potential
1107      REAL, DIMENSION(klon),      INTENT(IN)      ::  play          ! Pressure of layers [Pa]
1108      REAL, DIMENSION(klon),      INTENT(IN)      ::  ratqs         ! Parameter that determines the width of the water distrib     [-]
1109      REAL, DIMENSION(klon),      INTENT(IN)      ::  sigma_qtherm  ! Parameter determining the width of the distrib in thermals   [-]
1110      REAL, DIMENSION(klon),      INTENT(IN)      ::  qsth          ! Saturation specific humidity in thermals
1111      REAL, DIMENSION(klon),      INTENT(IN)      ::  qsenv         ! Saturation specific humidity in environment
1112     
1113      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  ctot         ! Cloud fraction [0-1]
1114      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  ctotth       ! Cloud fraction [0-1] in thermals
1115      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
1116      REAL, DIMENSION(klon),      INTENT(INOUT)      ::  qcloud       ! In cloud total water content [kg/kg]
1117      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sth    ! mean saturation deficit in thermals
1118      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_senv   ! mean saturation deficit in environment
1119      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmath  ! std of saturation deficit in thermals
1120      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmaenv ! std of saturation deficit in environment
1121
1122
1123! LOCAL VARIABLES
1124
1125      INTEGER itap,ind1,l,ig,iter,k
1126      INTEGER iflag_topthermals, niter
1127
1128      REAL qcth(klon)
1129      REAL qcenv(klon)
1130      REAL qctot(klon)
1131      REAL cth(klon)
1132      REAL cenv(klon)   
1133      REAL cth_vol(klon)
1134      REAL cenv_vol(klon)
1135      REAL qt_env(klon), thetal_env(klon)
1136      REAL sqrtpi,sqrt2,sqrt2pi
1137      REAL alth,alenv,ath,aenv
1138      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
1139      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
1140      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
1141      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
1142      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
1143      REAL zdelta,qsatbef,zcor
1144      REAL Tbefth(klon), Tbefenv(klon)
1145      REAL qlbef
1146      REAL dqsatenv(klon), dqsatth(klon)
1147      REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
1148      REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
1149      REAL qincloud(klon)
1150      REAL alenvl, aenvl
1151      REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
1152
1153
1154!------------------------------------------------------------------------------
1155! Initialisation
1156!------------------------------------------------------------------------------
1157
1158
1159      sqrt2pi=sqrt(2.*rpi)
1160      sqrt2=sqrt(2.)
1161      sqrtpi=sqrt(rpi)
1162
1163!-------------------------------------------------------------------------------
1164! Thermal fraction calculation and standard deviation of the distribution
1165!------------------------------------------------------------------------------- 
1166
1167! initialisations and calculation of temperature, humidity and saturation specific humidity
1168
1169cloudth_senv(:) = 0.
1170cloudth_sth(:) = 0.
1171cloudth_sigmaenv(:) = 0.
1172cloudth_sigmath(:) = 0.
1173
1174
1175DO ind1=1,klon
1176 
1177   Tbefenv(ind1) = temp(ind1)
1178   thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1)
1179   Tbefth(ind1)  = thetal_th(ind1)*zpspsk(ind1)
1180   qt_env(ind1)  = (qt(ind1)-frac_th(ind1)*qt_th(ind1))/(1.-frac_th(ind1)) !qt = a*qtth + (1-a)*qtenv
1181
1182ENDDO
1183
1184
1185
1186DO ind1=1,klon
1187
1188
1189    IF (frac_th(ind1).GT.min_frac_th_cld) THEN !Thermal and environnement
1190
1191! Environment:
1192
1193
1194        alenv=(RD/RV*RLVTT*qsenv(ind1))/(rd*thetal_env(ind1)**2)     
1195        aenv=1./(1.+(alenv*RLVTT/rcpd))                             
1196        senv=aenv*(qt_env(ind1)-qsenv(ind1))                           
1197     
1198
1199! Thermals:
1200
1201
1202        alth=(RD/RV*RLVTT*qsth(ind1))/(rd*thetal_th(ind1)**2)       
1203        ath=1./(1.+(alth*RLVTT/rcpd))                                                         
1204        sth=ath*(qt_th(ind1)-qsth(ind1))                     
1205
1206
1207! Standard deviation of the distributions
1208
1209           ! environment
1210           sigma1s_fraca = (sigma1s_factor**0.5)*(frac_th(ind1)**sigma1s_power) / &
1211           &                (1-frac_th(ind1))*((sth-senv)**2)**0.5
1212
1213           IF (cloudth_ratqsmin>0.) THEN
1214             sigma1s_ratqs = cloudth_ratqsmin*qt(ind1)
1215           ELSE
1216             sigma1s_ratqs = ratqs(ind1)*qt(ind1)
1217           ENDIF
1218           sigma1s = sigma1s_fraca + sigma1s_ratqs
1219
1220           IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then
1221              sigma1s = ratqs(ind1)*qt(ind1)*aenv
1222           ENDIF
1223
1224           ! thermals
1225           sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((frac_th(ind1)+0.02)**sigma2s_power))+0.002*qt_th(ind1)
1226
1227          IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1).ne.0) then
1228             sigma2s = sigma_qtherm(ind1)*ath
1229          ENDIF
1230
1231 
1232! surface cloud fraction
1233
1234           deltasenv=aenv*vert_alpha*sigma1s
1235           deltasth=ath*vert_alpha_th*sigma2s
1236
1237           xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1238           xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1239           exp_xenv1 = exp(-1.*xenv1**2)
1240           exp_xenv2 = exp(-1.*xenv2**2)
1241           xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1242           xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1243           exp_xth1 = exp(-1.*xth1**2)
1244           exp_xth2 = exp(-1.*xth2**2)
1245           cth(ind1)=0.5*(1.-1.*erf(xth1))
1246           cenv(ind1)=0.5*(1.-1.*erf(xenv1))
1247           ctot(ind1)=frac_th(ind1)*cth(ind1)+(1.-1.*frac_th(ind1))*cenv(ind1)
1248           ctotth(ind1)=frac_th(ind1)*cth(ind1)
1249       
1250
1251!volume cloud fraction and condensed water
1252
1253            !environnement
1254
1255            IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1256            IntJ_CF=0.5*(1.-1.*erf(xenv2))
1257
1258            IF (deltasenv .LT. 1.e-10) THEN
1259              qcenv(ind1)=IntJ
1260              cenv_vol(ind1)=IntJ_CF
1261            ELSE
1262              IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
1263              IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
1264              IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
1265              IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
1266              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
1267              qcenv(ind1)=IntJ+IntI1+IntI2+IntI3
1268              cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
1269            ENDIF
1270             
1271
1272
1273            !thermals
1274
1275            IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1276            IntJ_CF=0.5*(1.-1.*erf(xth2))
1277     
1278            IF (deltasth .LT. 1.e-10) THEN
1279              qcth(ind1)=IntJ
1280              cth_vol(ind1)=IntJ_CF
1281            ELSE
1282              IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1283              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1284              IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
1285              IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
1286              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
1287              qcth(ind1)=IntJ+IntI1+IntI2+IntI3
1288              cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
1289            ENDIF
1290
1291            ! total
1292
1293            qctot(ind1)=frac_th(ind1)*qcth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1)
1294            ctot_vol(ind1)=frac_th(ind1)*cth_vol(ind1)+(1.-1.*frac_th(ind1))*cenv_vol(ind1)
1295
1296            IF (cenv(ind1).LT.min_neb_th.and.cth(ind1).LT.min_neb_th) THEN
1297                ctot(ind1)=0.
1298                ctot_vol(ind1)=0.
1299                qcloud(ind1)=qsenv(ind1)
1300                qincloud(ind1)=0.
1301            ELSE             
1302                qincloud(ind1)=qctot(ind1)/ctot(ind1)
1303                !to prevent situations with cloud condensed water greater than available total water
1304                qincloud(ind1)=min(qincloud(ind1),qt(ind1)/ctot(ind1))
1305                ! we assume that water vapor in cloud is qsenv
1306                qcloud(ind1)=qincloud(ind1)+qsenv(ind1)
1307            ENDIF
1308
1309
1310
1311           ! Outputs used to check the PDFs
1312           cloudth_senv(ind1) = senv
1313           cloudth_sth(ind1) = sth
1314           cloudth_sigmaenv(ind1) = sigma1s
1315           cloudth_sigmath(ind1) = sigma2s
1316
1317      ENDIF       ! selection of grid points concerned by thermals
1318
1319
1320    ENDDO       !loop on klon
1321
1322
1323RETURN
1324
1325
1326END SUBROUTINE condensation_cloudth
1327
1328
1329!*****************************************************************************************
1330!*****************************************************************************************
1331! pre-cmip7 routines are below and are becoming obsolete
1332!*****************************************************************************************
1333!*****************************************************************************************
1334
1335
1336       SUBROUTINE cloudth(ngrid,klev,ind2,  &
1337     &           ztv,po,zqta,fraca, &
1338     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
1339     &           ratqs,zqs,t, &
1340     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1341
1342
1343      use lmdz_lscp_ini, only: iflag_cloudth_vert,iflag_ratqs
1344
1345      USE yomcst_mod_h
1346      USE yoethf_mod_h
1347IMPLICIT NONE
1348
1349
1350!===========================================================================
1351! Auteur : Arnaud Octavio Jam (LMD/CNRS)
1352! Date : 25 Mai 2010
1353! Objet : calcule les valeurs de qc et rneb dans les thermiques
1354!===========================================================================
1355
1356      INCLUDE "FCTTRE.h"
1357
1358      INTEGER itap,ind1,ind2
1359      INTEGER ngrid,klev,klon,l,ig
1360      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
1361     
1362      REAL ztv(ngrid,klev)
1363      REAL po(ngrid)
1364      REAL zqenv(ngrid)   
1365      REAL zqta(ngrid,klev)
1366         
1367      REAL fraca(ngrid,klev+1)
1368      REAL zpspsk(ngrid,klev)
1369      REAL paprs(ngrid,klev+1)
1370      REAL pplay(ngrid,klev)
1371      REAL ztla(ngrid,klev)
1372      REAL zthl(ngrid,klev)
1373
1374      REAL zqsatth(ngrid,klev)
1375      REAL zqsatenv(ngrid,klev)
1376     
1377     
1378      REAL sigma1(ngrid,klev)
1379      REAL sigma2(ngrid,klev)
1380      REAL qlth(ngrid,klev)
1381      REAL qlenv(ngrid,klev)
1382      REAL qltot(ngrid,klev)
1383      REAL cth(ngrid,klev) 
1384      REAL cenv(ngrid,klev)   
1385      REAL ctot(ngrid,klev)
1386      REAL rneb(ngrid,klev)
1387      REAL t(ngrid,klev)
1388      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
1389      REAL rdd,cppd,Lv
1390      REAL alth,alenv,ath,aenv
1391      REAL sth,senv,sigma1s,sigma2s,xth,xenv
1392      REAL Tbef,zdelta,qsatbef,zcor
1393      REAL qlbef 
1394      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
1395     
1396      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
1397      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
1398      REAL zqs(ngrid), qcloud(ngrid)
1399
1400
1401
1402
1403!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1404! Gestion de deux versions de cloudth
1405!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1406
1407      IF (iflag_cloudth_vert.GE.1) THEN
1408      CALL cloudth_vert(ngrid,klev,ind2,  &
1409     &           ztv,po,zqta,fraca, &
1410     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
1411     &           ratqs,zqs,t)
1412      RETURN
1413      ENDIF
1414!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1415
1416
1417!-------------------------------------------------------------------------------
1418! Initialisation des variables r?elles
1419!-------------------------------------------------------------------------------
1420      sigma1(:,ind2)=0.
1421      sigma2(:,ind2)=0.
1422      qlth(:,ind2)=0.
1423      qlenv(:,ind2)=0. 
1424      qltot(:,ind2)=0.
1425      rneb(:,ind2)=0.
1426      qcloud(:)=0.
1427      cth(:,ind2)=0.
1428      cenv(:,ind2)=0.
1429      ctot(:,ind2)=0.
1430      qsatmmussig1=0.
1431      qsatmmussig2=0.
1432      rdd=287.04
1433      cppd=1005.7
1434      pi=3.14159
1435      Lv=2.5e6
1436      sqrt2pi=sqrt(2.*pi)
1437
1438
1439
1440!-------------------------------------------------------------------------------
1441! Calcul de la fraction du thermique et des ?cart-types des distributions
1442!-------------------------------------------------------------------------------                 
1443      do ind1=1,ngrid
1444
1445      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
1446
1447      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
1448
1449
1450!      zqenv(ind1)=po(ind1)
1451      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1452      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1453      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1454      qsatbef=MIN(0.5,qsatbef)
1455      zcor=1./(1.-retv*qsatbef)
1456      qsatbef=qsatbef*zcor
1457      zqsatenv(ind1,ind2)=qsatbef
1458
1459
1460
1461
1462      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
1463      aenv=1./(1.+(alenv*Lv/cppd))
1464      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1465
1466
1467
1468
1469      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1470      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1471      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1472      qsatbef=MIN(0.5,qsatbef)
1473      zcor=1./(1.-retv*qsatbef)
1474      qsatbef=qsatbef*zcor
1475      zqsatth(ind1,ind2)=qsatbef
1476           
1477      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
1478      ath=1./(1.+(alth*Lv/cppd))
1479      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
1480     
1481     
1482
1483!------------------------------------------------------------------------------
1484! Calcul des ?cart-types pour s
1485!------------------------------------------------------------------------------
1486
1487!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
1488!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
1489!       if (paprs(ind1,ind2).gt.90000) then
1490!       ratqs(ind1,ind2)=0.002
1491!       else
1492!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
1493!       endif
1494       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
1495       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
1496!       sigma1s=ratqs(ind1,ind2)*po(ind1)
1497!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
1498 
1499!------------------------------------------------------------------------------
1500! Calcul de l'eau condens?e et de la couverture nuageuse
1501!------------------------------------------------------------------------------
1502      sqrt2pi=sqrt(2.*pi)
1503      xth=sth/(sqrt(2.)*sigma2s)
1504      xenv=senv/(sqrt(2.)*sigma1s)
1505      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
1506      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1507      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
1508
1509      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
1510      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
1511      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1512
1513!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1514      if (ctot(ind1,ind2).lt.1.e-10) then
1515      ctot(ind1,ind2)=0.
1516      qcloud(ind1)=zqsatenv(ind1,ind2)
1517
1518      else   
1519               
1520      ctot(ind1,ind2)=ctot(ind1,ind2)
1521      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1522
1523      endif                           
1524     
1525         
1526
1527
1528      else  ! gaussienne environnement seule
1529     
1530      zqenv(ind1)=po(ind1)
1531      Tbef=t(ind1,ind2)
1532      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1533      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1534      qsatbef=MIN(0.5,qsatbef)
1535      zcor=1./(1.-retv*qsatbef)
1536      qsatbef=qsatbef*zcor
1537      zqsatenv(ind1,ind2)=qsatbef
1538     
1539
1540!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1541      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
1542      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
1543      aenv=1./(1.+(alenv*Lv/cppd))
1544      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1545     
1546
1547      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1548
1549      sqrt2pi=sqrt(2.*pi)
1550      xenv=senv/(sqrt(2.)*sigma1s)
1551      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1552      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
1553     
1554      if (ctot(ind1,ind2).lt.1.e-3) then
1555      ctot(ind1,ind2)=0.
1556      qcloud(ind1)=zqsatenv(ind1,ind2)
1557
1558      else   
1559               
1560      ctot(ind1,ind2)=ctot(ind1,ind2)
1561      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1562
1563      endif   
1564 
1565 
1566 
1567 
1568 
1569 
1570      endif   
1571      enddo
1572     
1573      return
1574!     end
1575END SUBROUTINE cloudth
1576
1577
1578
1579!===========================================================================
1580     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
1581     &           ztv,po,zqta,fraca, &
1582     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
1583     &           ratqs,zqs,t)
1584
1585!===========================================================================
1586! Auteur : Arnaud Octavio Jam (LMD/CNRS)
1587! Date : 25 Mai 2010
1588! Objet : calcule les valeurs de qc et rneb dans les thermiques
1589!===========================================================================
1590
1591
1592USE yoethf_mod_h
1593            use lmdz_lscp_ini, only: iflag_cloudth_vert, vert_alpha
1594
1595      USE yomcst_mod_h
1596IMPLICIT NONE
1597
1598
1599      INCLUDE "FCTTRE.h"
1600     
1601      INTEGER itap,ind1,ind2
1602      INTEGER ngrid,klev,klon,l,ig
1603     
1604      REAL ztv(ngrid,klev)
1605      REAL po(ngrid)
1606      REAL zqenv(ngrid)   
1607      REAL zqta(ngrid,klev)
1608         
1609      REAL fraca(ngrid,klev+1)
1610      REAL zpspsk(ngrid,klev)
1611      REAL paprs(ngrid,klev+1)
1612      REAL pplay(ngrid,klev)
1613      REAL ztla(ngrid,klev)
1614      REAL zthl(ngrid,klev)
1615
1616      REAL zqsatth(ngrid,klev)
1617      REAL zqsatenv(ngrid,klev)
1618     
1619     
1620      REAL sigma1(ngrid,klev)                                                         
1621      REAL sigma2(ngrid,klev)
1622      REAL qlth(ngrid,klev)
1623      REAL qlenv(ngrid,klev)
1624      REAL qltot(ngrid,klev)
1625      REAL cth(ngrid,klev) 
1626      REAL cenv(ngrid,klev)   
1627      REAL ctot(ngrid,klev)
1628      REAL rneb(ngrid,klev)
1629      REAL t(ngrid,klev)                                                                 
1630      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
1631      REAL rdd,cppd,Lv,sqrt2,sqrtpi
1632      REAL alth,alenv,ath,aenv
1633      REAL sth,senv,sigma1s,sigma2s,xth,xenv
1634      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
1635      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
1636      REAL Tbef,zdelta,qsatbef,zcor
1637      REAL qlbef 
1638      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
1639      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
1640      ! (J Jouhaud, JL Dufresne, JB Madeleine)
1641     
1642      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
1643      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
1644      REAL zqs(ngrid), qcloud(ngrid)
1645
1646!------------------------------------------------------------------------------
1647! Initialisation des variables r?elles
1648!------------------------------------------------------------------------------
1649      sigma1(:,ind2)=0.
1650      sigma2(:,ind2)=0.
1651      qlth(:,ind2)=0.
1652      qlenv(:,ind2)=0. 
1653      qltot(:,ind2)=0.
1654      rneb(:,ind2)=0.
1655      qcloud(:)=0.
1656      cth(:,ind2)=0.
1657      cenv(:,ind2)=0.
1658      ctot(:,ind2)=0.
1659      qsatmmussig1=0.
1660      qsatmmussig2=0.
1661      rdd=287.04
1662      cppd=1005.7
1663      pi=3.14159
1664      Lv=2.5e6
1665      sqrt2pi=sqrt(2.*pi)
1666      sqrt2=sqrt(2.)
1667      sqrtpi=sqrt(pi)
1668
1669!-------------------------------------------------------------------------------
1670! Calcul de la fraction du thermique et des ?cart-types des distributions
1671!-------------------------------------------------------------------------------                 
1672      do ind1=1,ngrid
1673
1674      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
1675
1676      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
1677
1678
1679!      zqenv(ind1)=po(ind1)
1680      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1681      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1682      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1683      qsatbef=MIN(0.5,qsatbef)
1684      zcor=1./(1.-retv*qsatbef)
1685      qsatbef=qsatbef*zcor
1686      zqsatenv(ind1,ind2)=qsatbef
1687
1688
1689
1690
1691      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
1692      aenv=1./(1.+(alenv*Lv/cppd))
1693      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1694
1695
1696
1697
1698      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1699      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1700      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1701      qsatbef=MIN(0.5,qsatbef)
1702      zcor=1./(1.-retv*qsatbef)
1703      qsatbef=qsatbef*zcor
1704      zqsatth(ind1,ind2)=qsatbef
1705           
1706      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
1707      ath=1./(1.+(alth*Lv/cppd))
1708      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
1709     
1710     
1711
1712!------------------------------------------------------------------------------
1713! Calcul des ?cart-types pour s
1714!------------------------------------------------------------------------------
1715
1716      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
1717      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
1718!       if (paprs(ind1,ind2).gt.90000) then
1719!       ratqs(ind1,ind2)=0.002
1720!       else
1721!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
1722!       endif
1723!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
1724!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
1725!       sigma1s=ratqs(ind1,ind2)*po(ind1)
1726!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
1727 
1728!------------------------------------------------------------------------------
1729! Calcul de l'eau condens?e et de la couverture nuageuse
1730!------------------------------------------------------------------------------
1731      sqrt2pi=sqrt(2.*pi)
1732      xth=sth/(sqrt(2.)*sigma2s)
1733      xenv=senv/(sqrt(2.)*sigma1s)
1734      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
1735      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1736      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
1737
1738      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
1739      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
1740      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1741     
1742       IF (iflag_cloudth_vert == 1) THEN
1743!-------------------------------------------------------------------------------
1744!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
1745!-------------------------------------------------------------------------------
1746!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
1747!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
1748      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1749      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1750!      deltasenv=aenv*0.01*po(ind1)
1751!     deltasth=ath*0.01*zqta(ind1,ind2)   
1752      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
1753      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
1754      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
1755      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
1756      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
1757      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
1758     
1759      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
1760      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
1761      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
1762
1763      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
1764      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
1765      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
1766      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
1767
1768      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1769!      qlenv(ind1,ind2)=IntJ
1770!      print*, qlenv(ind1,ind2),'VERIF EAU'
1771
1772
1773      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
1774!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
1775!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
1776      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
1777      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
1778      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
1779      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1780!      qlth(ind1,ind2)=IntJ
1781!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
1782      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1783
1784      ELSE IF (iflag_cloudth_vert == 2) THEN
1785
1786!-------------------------------------------------------------------------------
1787!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
1788!-------------------------------------------------------------------------------
1789!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
1790!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
1791!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1792!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1793      deltasenv=aenv*vert_alpha*sigma1s
1794      deltasth=ath*vert_alpha*sigma2s
1795     
1796      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1797      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1798      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1799      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1800!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
1801!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
1802     
1803      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
1804      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
1805      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1806
1807      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
1808      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
1809      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
1810      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
1811
1812!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
1813!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
1814!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
1815
1816      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1817!      qlenv(ind1,ind2)=IntJ
1818!      print*, qlenv(ind1,ind2),'VERIF EAU'
1819
1820      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
1821      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1822      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
1823      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
1824     
1825      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1826!      qlth(ind1,ind2)=IntJ
1827!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
1828      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1829     
1830
1831
1832
1833      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
1834
1835!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1836
1837      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1838      ctot(ind1,ind2)=0.
1839      qcloud(ind1)=zqsatenv(ind1,ind2)
1840
1841      else
1842               
1843      ctot(ind1,ind2)=ctot(ind1,ind2)
1844      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1845!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1846!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1847
1848      endif 
1849                       
1850     
1851         
1852!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
1853
1854
1855      else  ! gaussienne environnement seule
1856     
1857      zqenv(ind1)=po(ind1)
1858      Tbef=t(ind1,ind2)
1859      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1860      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1861      qsatbef=MIN(0.5,qsatbef)
1862      zcor=1./(1.-retv*qsatbef)
1863      qsatbef=qsatbef*zcor
1864      zqsatenv(ind1,ind2)=qsatbef
1865     
1866
1867!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1868      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
1869      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
1870      aenv=1./(1.+(alenv*Lv/cppd))
1871      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1872     
1873
1874      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1875
1876      sqrt2pi=sqrt(2.*pi)
1877      xenv=senv/(sqrt(2.)*sigma1s)
1878      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1879      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
1880     
1881      if (ctot(ind1,ind2).lt.1.e-3) then
1882      ctot(ind1,ind2)=0.
1883      qcloud(ind1)=zqsatenv(ind1,ind2)
1884
1885      else   
1886               
1887      ctot(ind1,ind2)=ctot(ind1,ind2)
1888      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1889
1890      endif   
1891 
1892 
1893 
1894 
1895 
1896 
1897      endif   
1898      enddo
1899     
1900      return
1901!     end
1902END SUBROUTINE cloudth_vert
1903
1904
1905
1906
1907       SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
1908     &           ztv,po,zqta,fraca, &
1909     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1910     &           ratqs,sigma_qtherm,zqs,t, &
1911     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1912
1913      use lmdz_lscp_ini, only: iflag_cloudth_vert
1914
1915      USE yomcst_mod_h
1916      USE yoethf_mod_h
1917IMPLICIT NONE
1918
1919
1920!===========================================================================
1921! Author : Arnaud Octavio Jam (LMD/CNRS)
1922! Date : 25 Mai 2010
1923! Objet : calcule les valeurs de qc et rneb dans les thermiques
1924!===========================================================================
1925      INCLUDE "FCTTRE.h"
1926
1927      integer, intent(in) :: ind2
1928      integer, intent(in) :: ngrid,klev
1929     
1930      real, dimension(ngrid,klev), intent(in) :: ztv
1931      real, dimension(ngrid), intent(in) :: po
1932      real, dimension(ngrid,klev), intent(in) :: zqta
1933      real, dimension(ngrid,klev+1), intent(in) :: fraca
1934      real, dimension(ngrid), intent(out) :: qcloud
1935      real, dimension(ngrid,klev), intent(out) :: ctot
1936      real, dimension(ngrid,klev), intent(out) :: ctot_vol
1937      real, dimension(ngrid,klev), intent(in) :: zpspsk
1938      real, dimension(ngrid,klev+1), intent(in) :: paprs
1939      real, dimension(ngrid,klev), intent(in) :: pplay
1940      real, dimension(ngrid,klev), intent(in) :: ztla
1941      real, dimension(ngrid,klev), intent(inout) :: zthl
1942      real, dimension(ngrid,klev), intent(in) :: ratqs,sigma_qtherm
1943      real, dimension(ngrid), intent(in) :: zqs
1944      real, dimension(ngrid,klev), intent(in) :: t
1945      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
1946
1947
1948      REAL zqenv(ngrid)   
1949      REAL zqsatth(ngrid,klev)
1950      REAL zqsatenv(ngrid,klev)
1951     
1952      REAL sigma1(ngrid,klev)                                                         
1953      REAL sigma2(ngrid,klev)
1954      REAL qlth(ngrid,klev)
1955      REAL qlenv(ngrid,klev)
1956      REAL qltot(ngrid,klev)
1957      REAL cth(ngrid,klev)
1958      REAL cenv(ngrid,klev)   
1959      REAL cth_vol(ngrid,klev)
1960      REAL cenv_vol(ngrid,klev)
1961      REAL rneb(ngrid,klev)     
1962      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
1963      REAL rdd,cppd,Lv
1964      REAL alth,alenv,ath,aenv
1965      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
1966      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
1967      REAL Tbef,zdelta,qsatbef,zcor
1968      REAL qlbef 
1969      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
1970      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
1971
1972
1973      INTEGER :: ind1,l, ig
1974
1975      IF (iflag_cloudth_vert.GE.1) THEN
1976      CALL cloudth_vert_v3(ngrid,klev,ind2,  &
1977     &           ztv,po,zqta,fraca, &
1978     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1979     &           ratqs,sigma_qtherm,zqs,t, &
1980     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1981      RETURN
1982      ENDIF
1983!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1984
1985
1986!-------------------------------------------------------------------------------
1987! Initialisation des variables r?elles
1988!-------------------------------------------------------------------------------
1989      sigma1(:,ind2)=0.
1990      sigma2(:,ind2)=0.
1991      qlth(:,ind2)=0.
1992      qlenv(:,ind2)=0. 
1993      qltot(:,ind2)=0.
1994      rneb(:,ind2)=0.
1995      qcloud(:)=0.
1996      cth(:,ind2)=0.
1997      cenv(:,ind2)=0.
1998      ctot(:,ind2)=0.
1999      cth_vol(:,ind2)=0.
2000      cenv_vol(:,ind2)=0.
2001      ctot_vol(:,ind2)=0.
2002      qsatmmussig1=0.
2003      qsatmmussig2=0.
2004      rdd=287.04
2005      cppd=1005.7
2006      pi=3.14159
2007      Lv=2.5e6
2008      sqrt2pi=sqrt(2.*pi)
2009      sqrt2=sqrt(2.)
2010      sqrtpi=sqrt(pi)
2011
2012
2013!-------------------------------------------------------------------------------
2014! Cloud fraction in the thermals and standard deviation of the PDFs
2015!-------------------------------------------------------------------------------                 
2016      do ind1=1,ngrid
2017
2018      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
2019
2020      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
2021
2022      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
2023      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2024      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2025      qsatbef=MIN(0.5,qsatbef)
2026      zcor=1./(1.-retv*qsatbef)
2027      qsatbef=qsatbef*zcor
2028      zqsatenv(ind1,ind2)=qsatbef
2029
2030
2031      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
2032      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
2033      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
2034
2035!po = qt de l'environnement ET des thermique
2036!zqenv = qt environnement
2037!zqsatenv = qsat environnement
2038!zthl = Tl environnement
2039
2040
2041      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
2042      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2043      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2044      qsatbef=MIN(0.5,qsatbef)
2045      zcor=1./(1.-retv*qsatbef)
2046      qsatbef=qsatbef*zcor
2047      zqsatth(ind1,ind2)=qsatbef
2048           
2049      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
2050      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
2051      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
2052     
2053!zqta = qt thermals
2054!zqsatth = qsat thermals
2055!ztla = Tl thermals
2056
2057!------------------------------------------------------------------------------
2058! s standard deviations
2059!------------------------------------------------------------------------------
2060
2061!     tests
2062!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
2063!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
2064!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
2065!     final option
2066      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
2067      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
2068 
2069!------------------------------------------------------------------------------
2070! Condensed water and cloud cover
2071!------------------------------------------------------------------------------
2072      xth=sth/(sqrt2*sigma2s)
2073      xenv=senv/(sqrt2*sigma1s)
2074      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
2075      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
2076      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
2077      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
2078
2079      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
2080      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
2081      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
2082
2083      if (ctot(ind1,ind2).lt.1.e-10) then
2084      ctot(ind1,ind2)=0.
2085      qcloud(ind1)=zqsatenv(ind1,ind2)
2086      else
2087      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
2088      endif
2089
2090      else  ! Environnement only, follow the if l.110
2091     
2092      zqenv(ind1)=po(ind1)
2093      Tbef=t(ind1,ind2)
2094      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2095      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2096      qsatbef=MIN(0.5,qsatbef)
2097      zcor=1./(1.-retv*qsatbef)
2098      qsatbef=qsatbef*zcor
2099      zqsatenv(ind1,ind2)=qsatbef
2100
2101!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
2102      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
2103      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
2104      aenv=1./(1.+(alenv*Lv/cppd))
2105      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
2106     
2107      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
2108
2109      xenv=senv/(sqrt2*sigma1s)
2110      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
2111      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
2112      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
2113
2114      if (ctot(ind1,ind2).lt.1.e-3) then
2115      ctot(ind1,ind2)=0.
2116      qcloud(ind1)=zqsatenv(ind1,ind2)
2117      else   
2118      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
2119      endif
2120
2121
2122      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
2123      enddo       ! from the loop on ngrid l.108
2124      return
2125!     end
2126END SUBROUTINE cloudth_v3
2127
2128
2129
2130!===========================================================================
2131     SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
2132     &           ztv,po,zqta,fraca, &
2133     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
2134     &           ratqs,sigma_qtherm,zqs,t, &
2135     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
2136
2137!===========================================================================
2138! Auteur : Arnaud Octavio Jam (LMD/CNRS)
2139! Date : 25 Mai 2010
2140! Objet : calcule les valeurs de qc et rneb dans les thermiques
2141!===========================================================================
2142
2143      use yoethf_mod_h
2144      use lmdz_lscp_ini, only : iflag_cloudth_vert,iflag_ratqs
2145      use lmdz_lscp_ini, only : vert_alpha,vert_alpha_th, sigma1s_factor, sigma1s_power , sigma2s_factor , sigma2s_power , cloudth_ratqsmin , iflag_cloudth_vert_noratqs
2146
2147      USE yomcst_mod_h
2148IMPLICIT NONE
2149
2150
2151
2152
2153      INCLUDE "FCTTRE.h"
2154     
2155      INTEGER itap,ind1,ind2
2156      INTEGER ngrid,klev,klon,l,ig
2157      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
2158     
2159      REAL ztv(ngrid,klev)
2160      REAL po(ngrid)
2161      REAL zqenv(ngrid)   
2162      REAL zqta(ngrid,klev)
2163         
2164      REAL fraca(ngrid,klev+1)
2165      REAL zpspsk(ngrid,klev)
2166      REAL paprs(ngrid,klev+1)
2167      REAL pplay(ngrid,klev)
2168      REAL ztla(ngrid,klev)
2169      REAL zthl(ngrid,klev)
2170
2171      REAL zqsatth(ngrid,klev)
2172      REAL zqsatenv(ngrid,klev)
2173     
2174      REAL sigma1(ngrid,klev)                                                         
2175      REAL sigma2(ngrid,klev)
2176      REAL qlth(ngrid,klev)
2177      REAL qlenv(ngrid,klev)
2178      REAL qltot(ngrid,klev)
2179      REAL cth(ngrid,klev)
2180      REAL cenv(ngrid,klev)   
2181      REAL ctot(ngrid,klev)
2182      REAL cth_vol(ngrid,klev)
2183      REAL cenv_vol(ngrid,klev)
2184      REAL ctot_vol(ngrid,klev)
2185      REAL rneb(ngrid,klev)
2186      REAL t(ngrid,klev)                                                                 
2187      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
2188      REAL rdd,cppd,Lv
2189      REAL alth,alenv,ath,aenv
2190      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
2191      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
2192      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
2193      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
2194      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
2195      REAL Tbef,zdelta,qsatbef,zcor
2196      REAL qlbef 
2197      REAL ratqs(ngrid,klev),sigma_qtherm(ngrid,klev) ! determine la largeur de distribution de vapeur
2198      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
2199      ! (J Jouhaud, JL Dufresne, JB Madeleine)
2200
2201      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
2202      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
2203      REAL zqs(ngrid), qcloud(ngrid)
2204
2205      REAL rhodz(ngrid,klev)
2206      REAL zrho(ngrid,klev)
2207      REAL dz(ngrid,klev)
2208
2209      DO ind1 = 1, ngrid
2210        !Layer calculation
2211        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
2212        zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
2213        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
2214      END DO
2215
2216!------------------------------------------------------------------------------
2217! Initialize
2218!------------------------------------------------------------------------------
2219
2220      sigma1(:,ind2)=0.
2221      sigma2(:,ind2)=0.
2222      qlth(:,ind2)=0.
2223      qlenv(:,ind2)=0. 
2224      qltot(:,ind2)=0.
2225      rneb(:,ind2)=0.
2226      qcloud(:)=0.
2227      cth(:,ind2)=0.
2228      cenv(:,ind2)=0.
2229      ctot(:,ind2)=0.
2230      cth_vol(:,ind2)=0.
2231      cenv_vol(:,ind2)=0.
2232      ctot_vol(:,ind2)=0.
2233      qsatmmussig1=0.
2234      qsatmmussig2=0.
2235      rdd=287.04
2236      cppd=1005.7
2237      pi=3.14159
2238      Lv=2.5e6
2239      sqrt2pi=sqrt(2.*pi)
2240      sqrt2=sqrt(2.)
2241      sqrtpi=sqrt(pi)
2242
2243
2244
2245!-------------------------------------------------------------------------------
2246! Calcul de la fraction du thermique et des ecart-types des distributions
2247!-------------------------------------------------------------------------------                 
2248      do ind1=1,ngrid
2249
2250      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
2251
2252      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
2253
2254
2255      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
2256      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2257      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2258      qsatbef=MIN(0.5,qsatbef)
2259      zcor=1./(1.-retv*qsatbef)
2260      qsatbef=qsatbef*zcor
2261      zqsatenv(ind1,ind2)=qsatbef
2262
2263
2264      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
2265      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
2266      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
2267
2268!zqenv = qt environnement
2269!zqsatenv = qsat environnement
2270!zthl = Tl environnement
2271
2272
2273      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
2274      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2275      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2276      qsatbef=MIN(0.5,qsatbef)
2277      zcor=1./(1.-retv*qsatbef)
2278      qsatbef=qsatbef*zcor
2279      zqsatth(ind1,ind2)=qsatbef
2280           
2281      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
2282      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
2283      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
2284     
2285     
2286!zqta = qt thermals
2287!zqsatth = qsat thermals
2288!ztla = Tl thermals
2289!------------------------------------------------------------------------------
2290! s standard deviation
2291!------------------------------------------------------------------------------
2292
2293      sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
2294     &                (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
2295!     sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
2296      IF (cloudth_ratqsmin>0.) THEN
2297         sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
2298      ELSE
2299         sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
2300      ENDIF
2301      sigma1s = sigma1s_fraca + sigma1s_ratqs
2302      sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
2303      IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then
2304         sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
2305         IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1,ind2).ne.0) then
2306            sigma2s = sigma_qtherm(ind1,ind2)*ath
2307         ENDIF
2308      ENDIF
2309     
2310!      tests
2311!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
2312!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
2313!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
2314!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
2315!       if (paprs(ind1,ind2).gt.90000) then
2316!       ratqs(ind1,ind2)=0.002
2317!       else
2318!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
2319!       endif
2320!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
2321!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
2322!       sigma1s=ratqs(ind1,ind2)*po(ind1)
2323!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
2324 
2325       IF (iflag_cloudth_vert == 1) THEN
2326!-------------------------------------------------------------------------------
2327!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
2328!-------------------------------------------------------------------------------
2329
2330      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
2331      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
2332
2333      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
2334      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
2335      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
2336      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
2337      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
2338      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
2339     
2340      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
2341      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
2342      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
2343
2344      ! Environment
2345      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
2346      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
2347      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
2348      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
2349
2350      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
2351
2352      ! Thermal
2353      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
2354      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
2355      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
2356      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
2357      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
2358      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
2359
2360      ELSE IF (iflag_cloudth_vert >= 3) THEN
2361      IF (iflag_cloudth_vert < 5) THEN
2362!-------------------------------------------------------------------------------
2363!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
2364!-------------------------------------------------------------------------------
2365!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
2366!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
2367!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
2368!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
2369      IF (iflag_cloudth_vert == 3) THEN
2370        deltasenv=aenv*vert_alpha*sigma1s
2371        deltasth=ath*vert_alpha_th*sigma2s
2372      ELSE IF (iflag_cloudth_vert == 4) THEN
2373        IF (iflag_cloudth_vert_noratqs == 1) THEN
2374          deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
2375          deltasth=vert_alpha_th*sigma2s
2376        ELSE
2377          deltasenv=vert_alpha*sigma1s
2378          deltasth=vert_alpha_th*sigma2s
2379        ENDIF
2380      ENDIF
2381     
2382      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
2383      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
2384      exp_xenv1 = exp(-1.*xenv1**2)
2385      exp_xenv2 = exp(-1.*xenv2**2)
2386      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
2387      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
2388      exp_xth1 = exp(-1.*xth1**2)
2389      exp_xth2 = exp(-1.*xth2**2)
2390     
2391      !CF_surfacique
2392      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
2393      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
2394      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
2395
2396
2397      !CF_volumique & eau condense
2398      !environnement
2399      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
2400      IntJ_CF=0.5*(1.-1.*erf(xenv2))
2401      if (deltasenv .lt. 1.e-10) then
2402      qlenv(ind1,ind2)=IntJ
2403      cenv_vol(ind1,ind2)=IntJ_CF
2404      else
2405      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
2406      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
2407      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
2408      IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
2409      IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
2410      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
2411      cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
2412      endif
2413
2414      !thermique
2415      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
2416      IntJ_CF=0.5*(1.-1.*erf(xth2))
2417      if (deltasth .lt. 1.e-10) then
2418      qlth(ind1,ind2)=IntJ
2419      cth_vol(ind1,ind2)=IntJ_CF
2420      else
2421      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
2422      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
2423      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
2424      IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
2425      IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
2426      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
2427      cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
2428      endif
2429
2430      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
2431      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
2432
2433      ELSE IF (iflag_cloudth_vert == 5) THEN
2434         sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) &
2435              /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) &
2436              +ratqs(ind1,ind2)*po(ind1) !Environment
2437      sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)                   !Thermals
2438      !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
2439      !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
2440      xth=sth/(sqrt(2.)*sigma2s)
2441      xenv=senv/(sqrt(2.)*sigma1s)
2442
2443      !Volumique
2444      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
2445      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
2446      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
2447      !print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
2448
2449      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2))
2450      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2)) 
2451      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
2452
2453      !Surfacique
2454      !Neggers
2455      !beta=0.0044
2456      !inverse_rho=1.+beta*dz(ind1,ind2)
2457      !print *,'jeanjean : beta=',beta
2458      !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
2459      !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
2460      !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
2461
2462      !Brooks
2463      a_Brooks=0.6694
2464      b_Brooks=0.1882
2465      A_Maj_Brooks=0.1635 !-- sans shear
2466      !A_Maj_Brooks=0.17   !-- ARM LES
2467      !A_Maj_Brooks=0.18   !-- RICO LES
2468      !A_Maj_Brooks=0.19   !-- BOMEX LES
2469      Dx_Brooks=200000.
2470      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
2471      !print *,'jeanjean_f=',f_Brooks
2472
2473      cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
2474      cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
2475      ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
2476      !print *,'JJ_ctot_1',ctot(ind1,ind2)
2477
2478
2479
2480
2481
2482      ENDIF ! of if (iflag_cloudth_vert<5)
2483      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
2484
2485!      if (ctot(ind1,ind2).lt.1.e-10) then
2486      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
2487      ctot(ind1,ind2)=0.
2488      ctot_vol(ind1,ind2)=0.
2489      qcloud(ind1)=zqsatenv(ind1,ind2)
2490
2491      else
2492               
2493      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
2494!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
2495!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
2496
2497      endif 
2498
2499      else  ! gaussienne environnement seule
2500     
2501
2502      zqenv(ind1)=po(ind1)
2503      Tbef=t(ind1,ind2)
2504      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2505      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2506      qsatbef=MIN(0.5,qsatbef)
2507      zcor=1./(1.-retv*qsatbef)
2508      qsatbef=qsatbef*zcor
2509      zqsatenv(ind1,ind2)=qsatbef
2510     
2511
2512!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
2513      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
2514      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
2515      aenv=1./(1.+(alenv*Lv/cppd))
2516      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
2517      sth=0.
2518     
2519
2520      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
2521      sigma2s=0.
2522
2523      sqrt2pi=sqrt(2.*pi)
2524      xenv=senv/(sqrt(2.)*sigma1s)
2525      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
2526      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
2527      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
2528     
2529      if (ctot(ind1,ind2).lt.1.e-3) then
2530      ctot(ind1,ind2)=0.
2531      qcloud(ind1)=zqsatenv(ind1,ind2)
2532
2533      else   
2534               
2535!      ctot(ind1,ind2)=ctot(ind1,ind2)
2536      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
2537
2538      endif 
2539 
2540
2541
2542
2543      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
2544      ! Outputs used to check the PDFs
2545      cloudth_senv(ind1,ind2) = senv
2546      cloudth_sth(ind1,ind2) = sth
2547      cloudth_sigmaenv(ind1,ind2) = sigma1s
2548      cloudth_sigmath(ind1,ind2) = sigma2s
2549
2550      enddo       ! from the loop on ngrid l.333
2551      return
2552!     end
2553END SUBROUTINE cloudth_vert_v3
2554!
2555
2556
2557
2558
2559
2560
2561
2562
2563
2564
2565
2566       SUBROUTINE cloudth_v6(ngrid,klev,ind2,  &
2567     &           ztv,po,zqta,fraca, &
2568     &           qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
2569     &           ratqs,zqs,T, &
2570     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
2571
2572      USE yoethf_mod_h
2573      USE lmdz_lscp_ini, only: iflag_cloudth_vert
2574
2575      USE yomcst_mod_h
2576IMPLICIT NONE
2577
2578
2579
2580      INCLUDE "FCTTRE.h"
2581
2582
2583        !Domain variables
2584      INTEGER ngrid !indice Max lat-lon
2585      INTEGER klev  !indice Max alt
2586      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
2587      INTEGER ind1  !indice in [1:ngrid]
2588      INTEGER ind2  !indice in [1:klev]
2589        !thermal plume fraction
2590      REAL fraca(ngrid,klev+1)   !thermal plumes fraction in the gridbox
2591        !temperatures
2592      REAL T(ngrid,klev)       !temperature
2593      REAL zpspsk(ngrid,klev)  !factor (p/p0)**kappa (used for potential variables)
2594      REAL ztv(ngrid,klev)     !potential temperature (voir thermcell_env.F90)
2595      REAL ztla(ngrid,klev)    !liquid temperature in the thermals (Tl_th)
2596      REAL zthl(ngrid,klev)    !liquid temperature in the environment (Tl_env)
2597        !pressure
2598      REAL paprs(ngrid,klev+1)   !pressure at the interface of levels
2599      REAL pplay(ngrid,klev)     !pressure at the middle of the level
2600        !humidity
2601      REAL ratqs(ngrid,klev)   !width of the total water subgrid-scale distribution
2602      REAL po(ngrid)           !total water (qt)
2603      REAL zqenv(ngrid)        !total water in the environment (qt_env)
2604      REAL zqta(ngrid,klev)    !total water in the thermals (qt_th)
2605      REAL zqsatth(ngrid,klev)   !water saturation level in the thermals (q_sat_th)
2606      REAL zqsatenv(ngrid,klev)  !water saturation level in the environment (q_sat_env)
2607      REAL qlth(ngrid,klev)    !condensed water in the thermals
2608      REAL qlenv(ngrid,klev)   !condensed water in the environment
2609      REAL qltot(ngrid,klev)   !condensed water in the gridbox
2610        !cloud fractions
2611      REAL cth_vol(ngrid,klev)   !cloud fraction by volume in the thermals
2612      REAL cenv_vol(ngrid,klev)  !cloud fraction by volume in the environment
2613      REAL ctot_vol(ngrid,klev)  !cloud fraction by volume in the gridbox
2614      REAL cth_surf(ngrid,klev)  !cloud fraction by surface in the thermals
2615      REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment 
2616      REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
2617        !PDF of saturation deficit variables
2618      REAL rdd,cppd,Lv
2619      REAL Tbef,zdelta,qsatbef,zcor
2620      REAL alth,alenv,ath,aenv
2621      REAL sth,senv              !saturation deficits in the thermals and environment
2622      REAL sigma_env,sigma_th    !standard deviations of the biGaussian PDF
2623        !cloud fraction variables
2624      REAL xth,xenv
2625      REAL inverse_rho,beta                                  !Neggers et al. (2011) method
2626      REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
2627        !Incloud total water variables
2628      REAL zqs(ngrid)    !q_sat
2629      REAL qcloud(ngrid) !eau totale dans le nuage
2630        !Some arithmetic variables
2631      REAL  pi,sqrt2,sqrt2pi
2632        !Depth of the layer
2633      REAL dz(ngrid,klev)    !epaisseur de la couche en metre
2634      REAL rhodz(ngrid,klev)
2635      REAL zrho(ngrid,klev)
2636      DO ind1 = 1, ngrid
2637        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
2638        zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd          ![kg/m3]
2639        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2)            ![m]
2640      END DO
2641
2642!------------------------------------------------------------------------------
2643! Initialization
2644!------------------------------------------------------------------------------
2645      qlth(:,ind2)=0.
2646      qlenv(:,ind2)=0. 
2647      qltot(:,ind2)=0.
2648      cth_vol(:,ind2)=0.
2649      cenv_vol(:,ind2)=0.
2650      ctot_vol(:,ind2)=0.
2651      cth_surf(:,ind2)=0.
2652      cenv_surf(:,ind2)=0.
2653      ctot_surf(:,ind2)=0.
2654      qcloud(:)=0.
2655      rdd=287.04
2656      cppd=1005.7
2657      pi=3.14159
2658      Lv=2.5e6
2659      sqrt2=sqrt(2.)
2660      sqrt2pi=sqrt(2.*pi)
2661
2662
2663      DO ind1=1,ngrid
2664!-------------------------------------------------------------------------------
2665!Both thermal and environment in the gridbox
2666!-------------------------------------------------------------------------------
2667      IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
2668        !--------------------------------------------
2669        !calcul de qsat_env
2670        !--------------------------------------------
2671      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
2672      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2673      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2674      qsatbef=MIN(0.5,qsatbef)
2675      zcor=1./(1.-retv*qsatbef)
2676      qsatbef=qsatbef*zcor
2677      zqsatenv(ind1,ind2)=qsatbef
2678        !--------------------------------------------
2679        !calcul de s_env
2680        !--------------------------------------------
2681      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
2682      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
2683      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
2684        !--------------------------------------------
2685        !calcul de qsat_th
2686        !--------------------------------------------
2687      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
2688      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2689      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2690      qsatbef=MIN(0.5,qsatbef)
2691      zcor=1./(1.-retv*qsatbef)
2692      qsatbef=qsatbef*zcor
2693      zqsatth(ind1,ind2)=qsatbef
2694        !--------------------------------------------
2695        !calcul de s_th 
2696        !--------------------------------------------
2697      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84 these Arnaud Jam
2698      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84 these Arnaud Jam
2699      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84 these Arnaud Jam
2700        !--------------------------------------------
2701        !calcul standard deviations bi-Gaussian PDF
2702        !--------------------------------------------
2703      sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)
2704      sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) &
2705           /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) &
2706           +ratqs(ind1,ind2)*po(ind1)
2707      xth=sth/(sqrt2*sigma_th)
2708      xenv=senv/(sqrt2*sigma_env)
2709        !--------------------------------------------
2710        !Cloud fraction by volume CF_vol
2711        !--------------------------------------------
2712      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
2713      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
2714      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
2715        !--------------------------------------------
2716        !Condensed water qc
2717        !--------------------------------------------
2718      qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2))
2719      qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2)) 
2720      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
2721        !--------------------------------------------
2722        !Cloud fraction by surface CF_surf
2723        !--------------------------------------------
2724        !Method Neggers et al. (2011) : ok for cumulus clouds only
2725      !beta=0.0044 (Jouhaud et al.2018)
2726      !inverse_rho=1.+beta*dz(ind1,ind2)
2727      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
2728        !Method Brooks et al. (2005) : ok for all types of clouds
2729      a_Brooks=0.6694
2730      b_Brooks=0.1882
2731      A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
2732      Dx_Brooks=200000.   !-- si l'on considere des mailles de 200km de cote
2733      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
2734      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
2735        !--------------------------------------------
2736        !Incloud Condensed water qcloud
2737        !--------------------------------------------
2738      if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
2739      ctot_vol(ind1,ind2)=0.
2740      ctot_surf(ind1,ind2)=0.
2741      qcloud(ind1)=zqsatenv(ind1,ind2)
2742      else
2743      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
2744      endif
2745
2746
2747
2748!-------------------------------------------------------------------------------
2749!Environment only in the gridbox
2750!-------------------------------------------------------------------------------
2751      ELSE
2752        !--------------------------------------------
2753        !calcul de qsat_env
2754        !--------------------------------------------
2755      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
2756      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2757      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2758      qsatbef=MIN(0.5,qsatbef)
2759      zcor=1./(1.-retv*qsatbef)
2760      qsatbef=qsatbef*zcor
2761      zqsatenv(ind1,ind2)=qsatbef
2762        !--------------------------------------------
2763        !calcul de s_env
2764        !--------------------------------------------
2765      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
2766      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
2767      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
2768        !--------------------------------------------
2769        !calcul standard deviations Gaussian PDF
2770        !--------------------------------------------
2771      zqenv(ind1)=po(ind1)
2772      sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
2773      xenv=senv/(sqrt2*sigma_env)
2774        !--------------------------------------------
2775        !Cloud fraction by volume CF_vol
2776        !--------------------------------------------
2777      ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
2778        !--------------------------------------------
2779        !Condensed water qc
2780        !--------------------------------------------
2781      qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2))
2782        !--------------------------------------------
2783        !Cloud fraction by surface CF_surf
2784        !--------------------------------------------
2785        !Method Neggers et al. (2011) : ok for cumulus clouds only
2786      !beta=0.0044 (Jouhaud et al.2018)
2787      !inverse_rho=1.+beta*dz(ind1,ind2)
2788      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
2789        !Method Brooks et al. (2005) : ok for all types of clouds
2790      a_Brooks=0.6694
2791      b_Brooks=0.1882
2792      A_Maj_Brooks=0.1635 !-- sans dependence au shear
2793      Dx_Brooks=200000.
2794      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
2795      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
2796        !--------------------------------------------
2797        !Incloud Condensed water qcloud
2798        !--------------------------------------------
2799      if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
2800      ctot_vol(ind1,ind2)=0.
2801      ctot_surf(ind1,ind2)=0.
2802      qcloud(ind1)=zqsatenv(ind1,ind2)
2803      else
2804      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
2805      endif
2806
2807
2808      END IF  ! From the separation (thermal/envrionnement) et (environnement only)
2809
2810      ! Outputs used to check the PDFs
2811      cloudth_senv(ind1,ind2) = senv
2812      cloudth_sth(ind1,ind2) = sth
2813      cloudth_sigmaenv(ind1,ind2) = sigma_env
2814      cloudth_sigmath(ind1,ind2) = sigma_th
2815
2816      END DO  ! From the loop on ngrid
2817      return
2818
2819END SUBROUTINE cloudth_v6
2820
2821
2822
2823END MODULE lmdz_lscp_condensation
Note: See TracBrowser for help on using the repository browser.