source: LMDZ6/trunk/libf/phylmd/lmdz_cloudth.f90 @ 5441

Last change on this file since 5441 was 5285, checked in by abarral, 8 weeks ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

File size: 87.3 KB
Line 
1MODULE lmdz_cloudth
2
3
4  IMPLICIT NONE
5
6CONTAINS
7
8       SUBROUTINE cloudth(ngrid,klev,ind2,  &
9     &           ztv,po,zqta,fraca, &
10     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
11     &           ratqs,zqs,t, &
12     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
13
14
15      use lmdz_cloudth_ini, only: iflag_cloudth_vert,iflag_ratqs
16
17      USE yomcst_mod_h
18      USE yoethf_mod_h
19IMPLICIT NONE
20
21
22!===========================================================================
23! Auteur : Arnaud Octavio Jam (LMD/CNRS)
24! Date : 25 Mai 2010
25! Objet : calcule les valeurs de qc et rneb dans les thermiques
26!===========================================================================
27
28      INCLUDE "FCTTRE.h"
29
30      INTEGER itap,ind1,ind2
31      INTEGER ngrid,klev,klon,l,ig
32      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
33     
34      REAL ztv(ngrid,klev)
35      REAL po(ngrid)
36      REAL zqenv(ngrid)   
37      REAL zqta(ngrid,klev)
38         
39      REAL fraca(ngrid,klev+1)
40      REAL zpspsk(ngrid,klev)
41      REAL paprs(ngrid,klev+1)
42      REAL pplay(ngrid,klev)
43      REAL ztla(ngrid,klev)
44      REAL zthl(ngrid,klev)
45
46      REAL zqsatth(ngrid,klev)
47      REAL zqsatenv(ngrid,klev)
48     
49     
50      REAL sigma1(ngrid,klev)
51      REAL sigma2(ngrid,klev)
52      REAL qlth(ngrid,klev)
53      REAL qlenv(ngrid,klev)
54      REAL qltot(ngrid,klev)
55      REAL cth(ngrid,klev) 
56      REAL cenv(ngrid,klev)   
57      REAL ctot(ngrid,klev)
58      REAL rneb(ngrid,klev)
59      REAL t(ngrid,klev)
60      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
61      REAL rdd,cppd,Lv
62      REAL alth,alenv,ath,aenv
63      REAL sth,senv,sigma1s,sigma2s,xth,xenv
64      REAL Tbef,zdelta,qsatbef,zcor
65      REAL qlbef 
66      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
67     
68      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
69      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
70      REAL zqs(ngrid), qcloud(ngrid)
71      REAL erf
72
73
74
75
76!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
77! Gestion de deux versions de cloudth
78!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
79
80      IF (iflag_cloudth_vert.GE.1) THEN
81      CALL cloudth_vert(ngrid,klev,ind2,  &
82     &           ztv,po,zqta,fraca, &
83     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
84     &           ratqs,zqs,t)
85      RETURN
86      ENDIF
87!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
88
89
90!-------------------------------------------------------------------------------
91! Initialisation des variables r?elles
92!-------------------------------------------------------------------------------
93      sigma1(:,:)=0.
94      sigma2(:,:)=0.
95      qlth(:,:)=0.
96      qlenv(:,:)=0. 
97      qltot(:,:)=0.
98      rneb(:,:)=0.
99      qcloud(:)=0.
100      cth(:,:)=0.
101      cenv(:,:)=0.
102      ctot(:,:)=0.
103      qsatmmussig1=0.
104      qsatmmussig2=0.
105      rdd=287.04
106      cppd=1005.7
107      pi=3.14159
108      Lv=2.5e6
109      sqrt2pi=sqrt(2.*pi)
110
111
112
113!-------------------------------------------------------------------------------
114! Calcul de la fraction du thermique et des ?cart-types des distributions
115!-------------------------------------------------------------------------------                 
116      do ind1=1,ngrid
117
118      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
119
120      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
121
122
123!      zqenv(ind1)=po(ind1)
124      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
125      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
126      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
127      qsatbef=MIN(0.5,qsatbef)
128      zcor=1./(1.-retv*qsatbef)
129      qsatbef=qsatbef*zcor
130      zqsatenv(ind1,ind2)=qsatbef
131
132
133
134
135      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
136      aenv=1./(1.+(alenv*Lv/cppd))
137      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
138
139
140
141
142      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
143      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
144      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
145      qsatbef=MIN(0.5,qsatbef)
146      zcor=1./(1.-retv*qsatbef)
147      qsatbef=qsatbef*zcor
148      zqsatth(ind1,ind2)=qsatbef
149           
150      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
151      ath=1./(1.+(alth*Lv/cppd))
152      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
153     
154     
155
156!------------------------------------------------------------------------------
157! Calcul des ?cart-types pour s
158!------------------------------------------------------------------------------
159
160!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
161!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
162!       if (paprs(ind1,ind2).gt.90000) then
163!       ratqs(ind1,ind2)=0.002
164!       else
165!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
166!       endif
167       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
168       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
169!       sigma1s=ratqs(ind1,ind2)*po(ind1)
170!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
171 
172!------------------------------------------------------------------------------
173! Calcul de l'eau condens?e et de la couverture nuageuse
174!------------------------------------------------------------------------------
175      sqrt2pi=sqrt(2.*pi)
176      xth=sth/(sqrt(2.)*sigma2s)
177      xenv=senv/(sqrt(2.)*sigma1s)
178      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
179      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
180      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
181
182      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
183      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
184      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
185
186!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
187      if (ctot(ind1,ind2).lt.1.e-10) then
188      ctot(ind1,ind2)=0.
189      qcloud(ind1)=zqsatenv(ind1,ind2)
190
191      else   
192               
193      ctot(ind1,ind2)=ctot(ind1,ind2)
194      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
195
196      endif                           
197     
198         
199!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
200
201
202      else  ! gaussienne environnement seule
203     
204      zqenv(ind1)=po(ind1)
205      Tbef=t(ind1,ind2)
206      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
207      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
208      qsatbef=MIN(0.5,qsatbef)
209      zcor=1./(1.-retv*qsatbef)
210      qsatbef=qsatbef*zcor
211      zqsatenv(ind1,ind2)=qsatbef
212     
213
214!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
215      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
216      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
217      aenv=1./(1.+(alenv*Lv/cppd))
218      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
219     
220
221      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
222
223      sqrt2pi=sqrt(2.*pi)
224      xenv=senv/(sqrt(2.)*sigma1s)
225      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
226      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
227     
228      if (ctot(ind1,ind2).lt.1.e-3) then
229      ctot(ind1,ind2)=0.
230      qcloud(ind1)=zqsatenv(ind1,ind2)
231
232      else   
233               
234      ctot(ind1,ind2)=ctot(ind1,ind2)
235      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
236
237      endif   
238 
239 
240 
241 
242 
243 
244      endif   
245      enddo
246     
247      return
248!     end
249END SUBROUTINE cloudth
250
251
252
253!===========================================================================
254     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
255     &           ztv,po,zqta,fraca, &
256     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
257     &           ratqs,zqs,t)
258
259!===========================================================================
260! Auteur : Arnaud Octavio Jam (LMD/CNRS)
261! Date : 25 Mai 2010
262! Objet : calcule les valeurs de qc et rneb dans les thermiques
263!===========================================================================
264
265
266USE yoethf_mod_h
267            use lmdz_cloudth_ini, only: iflag_cloudth_vert, vert_alpha
268
269      USE yomcst_mod_h
270IMPLICIT NONE
271
272
273      INCLUDE "FCTTRE.h"
274     
275      INTEGER itap,ind1,ind2
276      INTEGER ngrid,klev,klon,l,ig
277     
278      REAL ztv(ngrid,klev)
279      REAL po(ngrid)
280      REAL zqenv(ngrid)   
281      REAL zqta(ngrid,klev)
282         
283      REAL fraca(ngrid,klev+1)
284      REAL zpspsk(ngrid,klev)
285      REAL paprs(ngrid,klev+1)
286      REAL pplay(ngrid,klev)
287      REAL ztla(ngrid,klev)
288      REAL zthl(ngrid,klev)
289
290      REAL zqsatth(ngrid,klev)
291      REAL zqsatenv(ngrid,klev)
292     
293     
294      REAL sigma1(ngrid,klev)                                                         
295      REAL sigma2(ngrid,klev)
296      REAL qlth(ngrid,klev)
297      REAL qlenv(ngrid,klev)
298      REAL qltot(ngrid,klev)
299      REAL cth(ngrid,klev) 
300      REAL cenv(ngrid,klev)   
301      REAL ctot(ngrid,klev)
302      REAL rneb(ngrid,klev)
303      REAL t(ngrid,klev)                                                                 
304      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
305      REAL rdd,cppd,Lv,sqrt2,sqrtpi
306      REAL alth,alenv,ath,aenv
307      REAL sth,senv,sigma1s,sigma2s,xth,xenv
308      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
309      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
310      REAL Tbef,zdelta,qsatbef,zcor
311      REAL qlbef 
312      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
313      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
314      ! (J Jouhaud, JL Dufresne, JB Madeleine)
315     
316      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
317      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
318      REAL zqs(ngrid), qcloud(ngrid)
319      REAL erf
320
321!------------------------------------------------------------------------------
322! Initialisation des variables r?elles
323!------------------------------------------------------------------------------
324      sigma1(:,:)=0.
325      sigma2(:,:)=0.
326      qlth(:,:)=0.
327      qlenv(:,:)=0. 
328      qltot(:,:)=0.
329      rneb(:,:)=0.
330      qcloud(:)=0.
331      cth(:,:)=0.
332      cenv(:,:)=0.
333      ctot(:,:)=0.
334      qsatmmussig1=0.
335      qsatmmussig2=0.
336      rdd=287.04
337      cppd=1005.7
338      pi=3.14159
339      Lv=2.5e6
340      sqrt2pi=sqrt(2.*pi)
341      sqrt2=sqrt(2.)
342      sqrtpi=sqrt(pi)
343
344!-------------------------------------------------------------------------------
345! Calcul de la fraction du thermique et des ?cart-types des distributions
346!-------------------------------------------------------------------------------                 
347      do ind1=1,ngrid
348
349      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
350
351      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
352
353
354!      zqenv(ind1)=po(ind1)
355      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
356      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
357      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
358      qsatbef=MIN(0.5,qsatbef)
359      zcor=1./(1.-retv*qsatbef)
360      qsatbef=qsatbef*zcor
361      zqsatenv(ind1,ind2)=qsatbef
362
363
364
365
366      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
367      aenv=1./(1.+(alenv*Lv/cppd))
368      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
369
370
371
372
373      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
374      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
375      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
376      qsatbef=MIN(0.5,qsatbef)
377      zcor=1./(1.-retv*qsatbef)
378      qsatbef=qsatbef*zcor
379      zqsatth(ind1,ind2)=qsatbef
380           
381      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
382      ath=1./(1.+(alth*Lv/cppd))
383      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
384     
385     
386
387!------------------------------------------------------------------------------
388! Calcul des ?cart-types pour s
389!------------------------------------------------------------------------------
390
391      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
392      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
393!       if (paprs(ind1,ind2).gt.90000) then
394!       ratqs(ind1,ind2)=0.002
395!       else
396!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
397!       endif
398!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
399!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
400!       sigma1s=ratqs(ind1,ind2)*po(ind1)
401!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
402 
403!------------------------------------------------------------------------------
404! Calcul de l'eau condens?e et de la couverture nuageuse
405!------------------------------------------------------------------------------
406      sqrt2pi=sqrt(2.*pi)
407      xth=sth/(sqrt(2.)*sigma2s)
408      xenv=senv/(sqrt(2.)*sigma1s)
409      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
410      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
411      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
412
413      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
414      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
415      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
416     
417       IF (iflag_cloudth_vert == 1) THEN
418!-------------------------------------------------------------------------------
419!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
420!-------------------------------------------------------------------------------
421!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
422!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
423      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
424      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
425!      deltasenv=aenv*0.01*po(ind1)
426!     deltasth=ath*0.01*zqta(ind1,ind2)   
427      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
428      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
429      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
430      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
431      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
432      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
433     
434      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
435      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
436      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
437
438      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
439      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
440      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
441      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
442
443      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
444!      qlenv(ind1,ind2)=IntJ
445!      print*, qlenv(ind1,ind2),'VERIF EAU'
446
447
448      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
449!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
450!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
451      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
452      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
453      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
454      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
455!      qlth(ind1,ind2)=IntJ
456!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
457      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
458
459      ELSE IF (iflag_cloudth_vert == 2) THEN
460
461!-------------------------------------------------------------------------------
462!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
463!-------------------------------------------------------------------------------
464!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
465!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
466!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
467!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
468      deltasenv=aenv*vert_alpha*sigma1s
469      deltasth=ath*vert_alpha*sigma2s
470     
471      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
472      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
473      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
474      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
475!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
476!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
477     
478      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
479      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
480      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
481
482      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
483      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
484      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
485      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
486
487!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
488!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
489!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
490
491      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
492!      qlenv(ind1,ind2)=IntJ
493!      print*, qlenv(ind1,ind2),'VERIF EAU'
494
495      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
496      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
497      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
498      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
499     
500      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
501!      qlth(ind1,ind2)=IntJ
502!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
503      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
504     
505
506
507
508      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
509
510!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
511
512      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
513      ctot(ind1,ind2)=0.
514      qcloud(ind1)=zqsatenv(ind1,ind2)
515
516      else
517               
518      ctot(ind1,ind2)=ctot(ind1,ind2)
519      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
520!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
521!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
522
523      endif 
524                       
525     
526         
527!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
528
529
530      else  ! gaussienne environnement seule
531     
532      zqenv(ind1)=po(ind1)
533      Tbef=t(ind1,ind2)
534      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
535      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
536      qsatbef=MIN(0.5,qsatbef)
537      zcor=1./(1.-retv*qsatbef)
538      qsatbef=qsatbef*zcor
539      zqsatenv(ind1,ind2)=qsatbef
540     
541
542!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
543      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
544      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
545      aenv=1./(1.+(alenv*Lv/cppd))
546      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
547     
548
549      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
550
551      sqrt2pi=sqrt(2.*pi)
552      xenv=senv/(sqrt(2.)*sigma1s)
553      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
554      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
555     
556      if (ctot(ind1,ind2).lt.1.e-3) then
557      ctot(ind1,ind2)=0.
558      qcloud(ind1)=zqsatenv(ind1,ind2)
559
560      else   
561               
562      ctot(ind1,ind2)=ctot(ind1,ind2)
563      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
564
565      endif   
566 
567 
568 
569 
570 
571 
572      endif   
573      enddo
574     
575      return
576!     end
577END SUBROUTINE cloudth_vert
578
579
580
581
582       SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
583     &           ztv,po,zqta,fraca, &
584     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
585     &           ratqs,sigma_qtherm,zqs,t, &
586     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
587
588      use lmdz_cloudth_ini, only: iflag_cloudth_vert
589
590      USE yomcst_mod_h
591      USE yoethf_mod_h
592IMPLICIT NONE
593
594
595!===========================================================================
596! Author : Arnaud Octavio Jam (LMD/CNRS)
597! Date : 25 Mai 2010
598! Objet : calcule les valeurs de qc et rneb dans les thermiques
599!===========================================================================
600      INCLUDE "FCTTRE.h"
601
602      integer, intent(in) :: ind2
603      integer, intent(in) :: ngrid,klev
604     
605      real, dimension(ngrid,klev), intent(in) :: ztv
606      real, dimension(ngrid), intent(in) :: po
607      real, dimension(ngrid,klev), intent(in) :: zqta
608      real, dimension(ngrid,klev+1), intent(in) :: fraca
609      real, dimension(ngrid), intent(out) :: qcloud
610      real, dimension(ngrid,klev), intent(out) :: ctot
611      real, dimension(ngrid,klev), intent(out) :: ctot_vol
612      real, dimension(ngrid,klev), intent(in) :: zpspsk
613      real, dimension(ngrid,klev+1), intent(in) :: paprs
614      real, dimension(ngrid,klev), intent(in) :: pplay
615      real, dimension(ngrid,klev), intent(in) :: ztla
616      real, dimension(ngrid,klev), intent(inout) :: zthl
617      real, dimension(ngrid,klev), intent(in) :: ratqs,sigma_qtherm
618      real, dimension(ngrid), intent(in) :: zqs
619      real, dimension(ngrid,klev), intent(in) :: t
620      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
621
622
623      REAL zqenv(ngrid)   
624      REAL zqsatth(ngrid,klev)
625      REAL zqsatenv(ngrid,klev)
626     
627      REAL sigma1(ngrid,klev)                                                         
628      REAL sigma2(ngrid,klev)
629      REAL qlth(ngrid,klev)
630      REAL qlenv(ngrid,klev)
631      REAL qltot(ngrid,klev)
632      REAL cth(ngrid,klev)
633      REAL cenv(ngrid,klev)   
634      REAL cth_vol(ngrid,klev)
635      REAL cenv_vol(ngrid,klev)
636      REAL rneb(ngrid,klev)     
637      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
638      REAL rdd,cppd,Lv
639      REAL alth,alenv,ath,aenv
640      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
641      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
642      REAL Tbef,zdelta,qsatbef,zcor
643      REAL qlbef 
644      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
645      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
646      REAL erf
647
648
649      INTEGER :: ind1,l, ig
650
651      IF (iflag_cloudth_vert.GE.1) THEN
652      CALL cloudth_vert_v3(ngrid,klev,ind2,  &
653     &           ztv,po,zqta,fraca, &
654     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
655     &           ratqs,sigma_qtherm,zqs,t, &
656     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
657      RETURN
658      ENDIF
659!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
660
661
662!-------------------------------------------------------------------------------
663! Initialisation des variables r?elles
664!-------------------------------------------------------------------------------
665      sigma1(:,:)=0.
666      sigma2(:,:)=0.
667      qlth(:,:)=0.
668      qlenv(:,:)=0. 
669      qltot(:,:)=0.
670      rneb(:,:)=0.
671      qcloud(:)=0.
672      cth(:,:)=0.
673      cenv(:,:)=0.
674      ctot(:,:)=0.
675      cth_vol(:,:)=0.
676      cenv_vol(:,:)=0.
677      ctot_vol(:,:)=0.
678      qsatmmussig1=0.
679      qsatmmussig2=0.
680      rdd=287.04
681      cppd=1005.7
682      pi=3.14159
683      Lv=2.5e6
684      sqrt2pi=sqrt(2.*pi)
685      sqrt2=sqrt(2.)
686      sqrtpi=sqrt(pi)
687
688
689!-------------------------------------------------------------------------------
690! Cloud fraction in the thermals and standard deviation of the PDFs
691!-------------------------------------------------------------------------------                 
692      do ind1=1,ngrid
693
694      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
695
696      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
697
698      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
699      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
700      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
701      qsatbef=MIN(0.5,qsatbef)
702      zcor=1./(1.-retv*qsatbef)
703      qsatbef=qsatbef*zcor
704      zqsatenv(ind1,ind2)=qsatbef
705
706
707      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
708      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
709      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
710
711!po = qt de l'environnement ET des thermique
712!zqenv = qt environnement
713!zqsatenv = qsat environnement
714!zthl = Tl environnement
715
716
717      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
718      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
719      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
720      qsatbef=MIN(0.5,qsatbef)
721      zcor=1./(1.-retv*qsatbef)
722      qsatbef=qsatbef*zcor
723      zqsatth(ind1,ind2)=qsatbef
724           
725      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
726      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
727      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
728     
729!zqta = qt thermals
730!zqsatth = qsat thermals
731!ztla = Tl thermals
732
733!------------------------------------------------------------------------------
734! s standard deviations
735!------------------------------------------------------------------------------
736
737!     tests
738!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
739!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
740!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
741!     final option
742      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
743      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
744 
745!------------------------------------------------------------------------------
746! Condensed water and cloud cover
747!------------------------------------------------------------------------------
748      xth=sth/(sqrt2*sigma2s)
749      xenv=senv/(sqrt2*sigma1s)
750      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
751      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
752      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
753      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
754
755      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
756      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
757      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
758
759      if (ctot(ind1,ind2).lt.1.e-10) then
760      ctot(ind1,ind2)=0.
761      qcloud(ind1)=zqsatenv(ind1,ind2)
762      else
763      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
764      endif
765
766      else  ! Environnement only, follow the if l.110
767     
768      zqenv(ind1)=po(ind1)
769      Tbef=t(ind1,ind2)
770      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
771      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
772      qsatbef=MIN(0.5,qsatbef)
773      zcor=1./(1.-retv*qsatbef)
774      qsatbef=qsatbef*zcor
775      zqsatenv(ind1,ind2)=qsatbef
776
777!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
778      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
779      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
780      aenv=1./(1.+(alenv*Lv/cppd))
781      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
782     
783      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
784
785      xenv=senv/(sqrt2*sigma1s)
786      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
787      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
788      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
789
790      if (ctot(ind1,ind2).lt.1.e-3) then
791      ctot(ind1,ind2)=0.
792      qcloud(ind1)=zqsatenv(ind1,ind2)
793      else   
794      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
795      endif
796
797
798      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
799      enddo       ! from the loop on ngrid l.108
800      return
801!     end
802END SUBROUTINE cloudth_v3
803
804
805
806!===========================================================================
807     SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
808     &           ztv,po,zqta,fraca, &
809     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
810     &           ratqs,sigma_qtherm,zqs,t, &
811     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
812
813!===========================================================================
814! Auteur : Arnaud Octavio Jam (LMD/CNRS)
815! Date : 25 Mai 2010
816! Objet : calcule les valeurs de qc et rneb dans les thermiques
817!===========================================================================
818
819USE yoethf_mod_h
820            use lmdz_cloudth_ini, only : iflag_cloudth_vert,iflag_ratqs
821      use lmdz_cloudth_ini, only : vert_alpha,vert_alpha_th, sigma1s_factor, sigma1s_power , sigma2s_factor , sigma2s_power , cloudth_ratqsmin , iflag_cloudth_vert_noratqs
822
823      USE yomcst_mod_h
824IMPLICIT NONE
825
826
827
828
829      INCLUDE "FCTTRE.h"
830     
831      INTEGER itap,ind1,ind2
832      INTEGER ngrid,klev,klon,l,ig
833      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
834     
835      REAL ztv(ngrid,klev)
836      REAL po(ngrid)
837      REAL zqenv(ngrid)   
838      REAL zqta(ngrid,klev)
839         
840      REAL fraca(ngrid,klev+1)
841      REAL zpspsk(ngrid,klev)
842      REAL paprs(ngrid,klev+1)
843      REAL pplay(ngrid,klev)
844      REAL ztla(ngrid,klev)
845      REAL zthl(ngrid,klev)
846
847      REAL zqsatth(ngrid,klev)
848      REAL zqsatenv(ngrid,klev)
849     
850      REAL sigma1(ngrid,klev)                                                         
851      REAL sigma2(ngrid,klev)
852      REAL qlth(ngrid,klev)
853      REAL qlenv(ngrid,klev)
854      REAL qltot(ngrid,klev)
855      REAL cth(ngrid,klev)
856      REAL cenv(ngrid,klev)   
857      REAL ctot(ngrid,klev)
858      REAL cth_vol(ngrid,klev)
859      REAL cenv_vol(ngrid,klev)
860      REAL ctot_vol(ngrid,klev)
861      REAL rneb(ngrid,klev)
862      REAL t(ngrid,klev)                                                                 
863      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
864      REAL rdd,cppd,Lv
865      REAL alth,alenv,ath,aenv
866      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
867      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
868      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
869      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
870      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
871      REAL Tbef,zdelta,qsatbef,zcor
872      REAL qlbef 
873      REAL ratqs(ngrid,klev),sigma_qtherm(ngrid,klev) ! determine la largeur de distribution de vapeur
874      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
875      ! (J Jouhaud, JL Dufresne, JB Madeleine)
876
877      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
878      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
879      REAL zqs(ngrid), qcloud(ngrid)
880      REAL erf
881
882      REAL rhodz(ngrid,klev)
883      REAL zrho(ngrid,klev)
884      REAL dz(ngrid,klev)
885
886      DO ind1 = 1, ngrid
887        !Layer calculation
888        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
889        zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
890        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
891      END DO
892
893!------------------------------------------------------------------------------
894! Initialize
895!------------------------------------------------------------------------------
896
897      sigma1(:,:)=0.
898      sigma2(:,:)=0.
899      qlth(:,:)=0.
900      qlenv(:,:)=0. 
901      qltot(:,:)=0.
902      rneb(:,:)=0.
903      qcloud(:)=0.
904      cth(:,:)=0.
905      cenv(:,:)=0.
906      ctot(:,:)=0.
907      cth_vol(:,:)=0.
908      cenv_vol(:,:)=0.
909      ctot_vol(:,:)=0.
910      qsatmmussig1=0.
911      qsatmmussig2=0.
912      rdd=287.04
913      cppd=1005.7
914      pi=3.14159
915      Lv=2.5e6
916      sqrt2pi=sqrt(2.*pi)
917      sqrt2=sqrt(2.)
918      sqrtpi=sqrt(pi)
919
920
921
922!-------------------------------------------------------------------------------
923! Calcul de la fraction du thermique et des ecart-types des distributions
924!-------------------------------------------------------------------------------                 
925      do ind1=1,ngrid
926
927      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
928
929      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
930
931
932      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
933      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
934      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
935      qsatbef=MIN(0.5,qsatbef)
936      zcor=1./(1.-retv*qsatbef)
937      qsatbef=qsatbef*zcor
938      zqsatenv(ind1,ind2)=qsatbef
939
940
941      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
942      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
943      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
944
945!zqenv = qt environnement
946!zqsatenv = qsat environnement
947!zthl = Tl environnement
948
949
950      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
951      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
952      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
953      qsatbef=MIN(0.5,qsatbef)
954      zcor=1./(1.-retv*qsatbef)
955      qsatbef=qsatbef*zcor
956      zqsatth(ind1,ind2)=qsatbef
957           
958      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
959      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
960      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
961     
962     
963!zqta = qt thermals
964!zqsatth = qsat thermals
965!ztla = Tl thermals
966!------------------------------------------------------------------------------
967! s standard deviation
968!------------------------------------------------------------------------------
969
970      sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
971     &                (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
972!     sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
973      IF (cloudth_ratqsmin>0.) THEN
974         sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
975      ELSE
976         sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
977      ENDIF
978      sigma1s = sigma1s_fraca + sigma1s_ratqs
979      sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
980      IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then
981         sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
982         IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1,ind2).ne.0) then
983            sigma2s = sigma_qtherm(ind1,ind2)*ath
984         ENDIF
985      ENDIF
986     
987!      tests
988!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
989!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
990!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
991!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
992!       if (paprs(ind1,ind2).gt.90000) then
993!       ratqs(ind1,ind2)=0.002
994!       else
995!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
996!       endif
997!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
998!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
999!       sigma1s=ratqs(ind1,ind2)*po(ind1)
1000!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
1001 
1002       IF (iflag_cloudth_vert == 1) THEN
1003!-------------------------------------------------------------------------------
1004!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
1005!-------------------------------------------------------------------------------
1006
1007      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1008      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1009
1010      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
1011      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
1012      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
1013      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
1014      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
1015      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
1016     
1017      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
1018      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
1019      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
1020
1021      ! Environment
1022      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
1023      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
1024      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
1025      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
1026
1027      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1028
1029      ! Thermal
1030      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
1031      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
1032      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
1033      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
1034      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1035      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1036
1037      ELSE IF (iflag_cloudth_vert >= 3) THEN
1038      IF (iflag_cloudth_vert < 5) THEN
1039!-------------------------------------------------------------------------------
1040!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1041!-------------------------------------------------------------------------------
1042!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
1043!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
1044!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1045!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1046      IF (iflag_cloudth_vert == 3) THEN
1047        deltasenv=aenv*vert_alpha*sigma1s
1048        deltasth=ath*vert_alpha_th*sigma2s
1049      ELSE IF (iflag_cloudth_vert == 4) THEN
1050        IF (iflag_cloudth_vert_noratqs == 1) THEN
1051          deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
1052          deltasth=vert_alpha_th*sigma2s
1053        ELSE
1054          deltasenv=vert_alpha*sigma1s
1055          deltasth=vert_alpha_th*sigma2s
1056        ENDIF
1057      ENDIF
1058     
1059      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1060      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1061      exp_xenv1 = exp(-1.*xenv1**2)
1062      exp_xenv2 = exp(-1.*xenv2**2)
1063      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1064      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1065      exp_xth1 = exp(-1.*xth1**2)
1066      exp_xth2 = exp(-1.*xth2**2)
1067     
1068      !CF_surfacique
1069      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
1070      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
1071      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1072
1073
1074      !CF_volumique & eau condense
1075      !environnement
1076      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1077      IntJ_CF=0.5*(1.-1.*erf(xenv2))
1078      if (deltasenv .lt. 1.e-10) then
1079      qlenv(ind1,ind2)=IntJ
1080      cenv_vol(ind1,ind2)=IntJ_CF
1081      else
1082      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
1083      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
1084      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
1085      IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
1086      IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
1087      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1088      cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1089      endif
1090
1091      !thermique
1092      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1093      IntJ_CF=0.5*(1.-1.*erf(xth2))
1094      if (deltasth .lt. 1.e-10) then
1095      qlth(ind1,ind2)=IntJ
1096      cth_vol(ind1,ind2)=IntJ_CF
1097      else
1098      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1099      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1100      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
1101      IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
1102      IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
1103      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1104      cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1105      endif
1106
1107      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1108      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1109
1110      ELSE IF (iflag_cloudth_vert == 5) THEN
1111         sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) &
1112              /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) &
1113              +ratqs(ind1,ind2)*po(ind1) !Environment
1114      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
1115      !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
1116      !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
1117      xth=sth/(sqrt(2.)*sigma2s)
1118      xenv=senv/(sqrt(2.)*sigma1s)
1119
1120      !Volumique
1121      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
1122      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1123      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1124      !print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
1125
1126      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2))
1127      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2)) 
1128      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1129
1130      !Surfacique
1131      !Neggers
1132      !beta=0.0044
1133      !inverse_rho=1.+beta*dz(ind1,ind2)
1134      !print *,'jeanjean : beta=',beta
1135      !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
1136      !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
1137      !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1138
1139      !Brooks
1140      a_Brooks=0.6694
1141      b_Brooks=0.1882
1142      A_Maj_Brooks=0.1635 !-- sans shear
1143      !A_Maj_Brooks=0.17   !-- ARM LES
1144      !A_Maj_Brooks=0.18   !-- RICO LES
1145      !A_Maj_Brooks=0.19   !-- BOMEX LES
1146      Dx_Brooks=200000.
1147      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
1148      !print *,'jeanjean_f=',f_Brooks
1149
1150      cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
1151      cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
1152      ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1153      !print *,'JJ_ctot_1',ctot(ind1,ind2)
1154
1155
1156
1157
1158
1159      ENDIF ! of if (iflag_cloudth_vert<5)
1160      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1161
1162!      if (ctot(ind1,ind2).lt.1.e-10) then
1163      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1164      ctot(ind1,ind2)=0.
1165      ctot_vol(ind1,ind2)=0.
1166      qcloud(ind1)=zqsatenv(ind1,ind2)
1167
1168      else
1169               
1170      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1171!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1172!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1173
1174      endif 
1175
1176      else  ! gaussienne environnement seule
1177     
1178
1179      zqenv(ind1)=po(ind1)
1180      Tbef=t(ind1,ind2)
1181      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1182      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1183      qsatbef=MIN(0.5,qsatbef)
1184      zcor=1./(1.-retv*qsatbef)
1185      qsatbef=qsatbef*zcor
1186      zqsatenv(ind1,ind2)=qsatbef
1187     
1188
1189!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1190      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
1191      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
1192      aenv=1./(1.+(alenv*Lv/cppd))
1193      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1194      sth=0.
1195     
1196
1197      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1198      sigma2s=0.
1199
1200      sqrt2pi=sqrt(2.*pi)
1201      xenv=senv/(sqrt(2.)*sigma1s)
1202      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1203      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1204      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
1205     
1206      if (ctot(ind1,ind2).lt.1.e-3) then
1207      ctot(ind1,ind2)=0.
1208      qcloud(ind1)=zqsatenv(ind1,ind2)
1209
1210      else   
1211               
1212!      ctot(ind1,ind2)=ctot(ind1,ind2)
1213      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1214
1215      endif 
1216 
1217
1218
1219
1220      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1221      ! Outputs used to check the PDFs
1222      cloudth_senv(ind1,ind2) = senv
1223      cloudth_sth(ind1,ind2) = sth
1224      cloudth_sigmaenv(ind1,ind2) = sigma1s
1225      cloudth_sigmath(ind1,ind2) = sigma2s
1226
1227      enddo       ! from the loop on ngrid l.333
1228      return
1229!     end
1230END SUBROUTINE cloudth_vert_v3
1231!
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243       SUBROUTINE cloudth_v6(ngrid,klev,ind2,  &
1244     &           ztv,po,zqta,fraca, &
1245     &           qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1246     &           ratqs,zqs,T, &
1247     &           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1248
1249USE yoethf_mod_h
1250            use lmdz_cloudth_ini, only: iflag_cloudth_vert
1251
1252      USE yomcst_mod_h
1253IMPLICIT NONE
1254
1255
1256
1257      INCLUDE "FCTTRE.h"
1258
1259
1260        !Domain variables
1261      INTEGER ngrid !indice Max lat-lon
1262      INTEGER klev  !indice Max alt
1263      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
1264      INTEGER ind1  !indice in [1:ngrid]
1265      INTEGER ind2  !indice in [1:klev]
1266        !thermal plume fraction
1267      REAL fraca(ngrid,klev+1)   !thermal plumes fraction in the gridbox
1268        !temperatures
1269      REAL T(ngrid,klev)       !temperature
1270      REAL zpspsk(ngrid,klev)  !factor (p/p0)**kappa (used for potential variables)
1271      REAL ztv(ngrid,klev)     !potential temperature (voir thermcell_env.F90)
1272      REAL ztla(ngrid,klev)    !liquid temperature in the thermals (Tl_th)
1273      REAL zthl(ngrid,klev)    !liquid temperature in the environment (Tl_env)
1274        !pressure
1275      REAL paprs(ngrid,klev+1)   !pressure at the interface of levels
1276      REAL pplay(ngrid,klev)     !pressure at the middle of the level
1277        !humidity
1278      REAL ratqs(ngrid,klev)   !width of the total water subgrid-scale distribution
1279      REAL po(ngrid)           !total water (qt)
1280      REAL zqenv(ngrid)        !total water in the environment (qt_env)
1281      REAL zqta(ngrid,klev)    !total water in the thermals (qt_th)
1282      REAL zqsatth(ngrid,klev)   !water saturation level in the thermals (q_sat_th)
1283      REAL zqsatenv(ngrid,klev)  !water saturation level in the environment (q_sat_env)
1284      REAL qlth(ngrid,klev)    !condensed water in the thermals
1285      REAL qlenv(ngrid,klev)   !condensed water in the environment
1286      REAL qltot(ngrid,klev)   !condensed water in the gridbox
1287        !cloud fractions
1288      REAL cth_vol(ngrid,klev)   !cloud fraction by volume in the thermals
1289      REAL cenv_vol(ngrid,klev)  !cloud fraction by volume in the environment
1290      REAL ctot_vol(ngrid,klev)  !cloud fraction by volume in the gridbox
1291      REAL cth_surf(ngrid,klev)  !cloud fraction by surface in the thermals
1292      REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment 
1293      REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
1294        !PDF of saturation deficit variables
1295      REAL rdd,cppd,Lv
1296      REAL Tbef,zdelta,qsatbef,zcor
1297      REAL alth,alenv,ath,aenv
1298      REAL sth,senv              !saturation deficits in the thermals and environment
1299      REAL sigma_env,sigma_th    !standard deviations of the biGaussian PDF
1300        !cloud fraction variables
1301      REAL xth,xenv
1302      REAL inverse_rho,beta                                  !Neggers et al. (2011) method
1303      REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
1304        !Incloud total water variables
1305      REAL zqs(ngrid)    !q_sat
1306      REAL qcloud(ngrid) !eau totale dans le nuage
1307        !Some arithmetic variables
1308      REAL erf,pi,sqrt2,sqrt2pi
1309        !Depth of the layer
1310      REAL dz(ngrid,klev)    !epaisseur de la couche en metre
1311      REAL rhodz(ngrid,klev)
1312      REAL zrho(ngrid,klev)
1313      DO ind1 = 1, ngrid
1314        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
1315        zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd          ![kg/m3]
1316        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2)            ![m]
1317      END DO
1318
1319!------------------------------------------------------------------------------
1320! Initialization
1321!------------------------------------------------------------------------------
1322      qlth(:,:)=0.
1323      qlenv(:,:)=0. 
1324      qltot(:,:)=0.
1325      cth_vol(:,:)=0.
1326      cenv_vol(:,:)=0.
1327      ctot_vol(:,:)=0.
1328      cth_surf(:,:)=0.
1329      cenv_surf(:,:)=0.
1330      ctot_surf(:,:)=0.
1331      qcloud(:)=0.
1332      rdd=287.04
1333      cppd=1005.7
1334      pi=3.14159
1335      Lv=2.5e6
1336      sqrt2=sqrt(2.)
1337      sqrt2pi=sqrt(2.*pi)
1338
1339
1340      DO ind1=1,ngrid
1341!-------------------------------------------------------------------------------
1342!Both thermal and environment in the gridbox
1343!-------------------------------------------------------------------------------
1344      IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
1345        !--------------------------------------------
1346        !calcul de qsat_env
1347        !--------------------------------------------
1348      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1349      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1350      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1351      qsatbef=MIN(0.5,qsatbef)
1352      zcor=1./(1.-retv*qsatbef)
1353      qsatbef=qsatbef*zcor
1354      zqsatenv(ind1,ind2)=qsatbef
1355        !--------------------------------------------
1356        !calcul de s_env
1357        !--------------------------------------------
1358      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
1359      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
1360      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
1361        !--------------------------------------------
1362        !calcul de qsat_th
1363        !--------------------------------------------
1364      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1365      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1366      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1367      qsatbef=MIN(0.5,qsatbef)
1368      zcor=1./(1.-retv*qsatbef)
1369      qsatbef=qsatbef*zcor
1370      zqsatth(ind1,ind2)=qsatbef
1371        !--------------------------------------------
1372        !calcul de s_th 
1373        !--------------------------------------------
1374      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84 these Arnaud Jam
1375      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84 these Arnaud Jam
1376      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84 these Arnaud Jam
1377        !--------------------------------------------
1378        !calcul standard deviations bi-Gaussian PDF
1379        !--------------------------------------------
1380      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)
1381      sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) &
1382           /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) &
1383           +ratqs(ind1,ind2)*po(ind1)
1384      xth=sth/(sqrt2*sigma_th)
1385      xenv=senv/(sqrt2*sigma_env)
1386        !--------------------------------------------
1387        !Cloud fraction by volume CF_vol
1388        !--------------------------------------------
1389      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
1390      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1391      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1392        !--------------------------------------------
1393        !Condensed water qc
1394        !--------------------------------------------
1395      qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2))
1396      qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2)) 
1397      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1398        !--------------------------------------------
1399        !Cloud fraction by surface CF_surf
1400        !--------------------------------------------
1401        !Method Neggers et al. (2011) : ok for cumulus clouds only
1402      !beta=0.0044 (Jouhaud et al.2018)
1403      !inverse_rho=1.+beta*dz(ind1,ind2)
1404      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1405        !Method Brooks et al. (2005) : ok for all types of clouds
1406      a_Brooks=0.6694
1407      b_Brooks=0.1882
1408      A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
1409      Dx_Brooks=200000.   !-- si l'on considere des mailles de 200km de cote
1410      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
1411      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1412        !--------------------------------------------
1413        !Incloud Condensed water qcloud
1414        !--------------------------------------------
1415      if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
1416      ctot_vol(ind1,ind2)=0.
1417      ctot_surf(ind1,ind2)=0.
1418      qcloud(ind1)=zqsatenv(ind1,ind2)
1419      else
1420      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
1421      endif
1422
1423
1424
1425!-------------------------------------------------------------------------------
1426!Environment only in the gridbox
1427!-------------------------------------------------------------------------------
1428      ELSE
1429        !--------------------------------------------
1430        !calcul de qsat_env
1431        !--------------------------------------------
1432      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1433      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1434      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1435      qsatbef=MIN(0.5,qsatbef)
1436      zcor=1./(1.-retv*qsatbef)
1437      qsatbef=qsatbef*zcor
1438      zqsatenv(ind1,ind2)=qsatbef
1439        !--------------------------------------------
1440        !calcul de s_env
1441        !--------------------------------------------
1442      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
1443      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
1444      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
1445        !--------------------------------------------
1446        !calcul standard deviations Gaussian PDF
1447        !--------------------------------------------
1448      zqenv(ind1)=po(ind1)
1449      sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
1450      xenv=senv/(sqrt2*sigma_env)
1451        !--------------------------------------------
1452        !Cloud fraction by volume CF_vol
1453        !--------------------------------------------
1454      ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1455        !--------------------------------------------
1456        !Condensed water qc
1457        !--------------------------------------------
1458      qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2))
1459        !--------------------------------------------
1460        !Cloud fraction by surface CF_surf
1461        !--------------------------------------------
1462        !Method Neggers et al. (2011) : ok for cumulus clouds only
1463      !beta=0.0044 (Jouhaud et al.2018)
1464      !inverse_rho=1.+beta*dz(ind1,ind2)
1465      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1466        !Method Brooks et al. (2005) : ok for all types of clouds
1467      a_Brooks=0.6694
1468      b_Brooks=0.1882
1469      A_Maj_Brooks=0.1635 !-- sans dependence au shear
1470      Dx_Brooks=200000.
1471      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
1472      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1473        !--------------------------------------------
1474        !Incloud Condensed water qcloud
1475        !--------------------------------------------
1476      if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
1477      ctot_vol(ind1,ind2)=0.
1478      ctot_surf(ind1,ind2)=0.
1479      qcloud(ind1)=zqsatenv(ind1,ind2)
1480      else
1481      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
1482      endif
1483
1484
1485      END IF  ! From the separation (thermal/envrionnement) et (environnement only)
1486
1487      ! Outputs used to check the PDFs
1488      cloudth_senv(ind1,ind2) = senv
1489      cloudth_sth(ind1,ind2) = sth
1490      cloudth_sigmaenv(ind1,ind2) = sigma_env
1491      cloudth_sigmath(ind1,ind2) = sigma_th
1492
1493      END DO  ! From the loop on ngrid
1494      return
1495
1496END SUBROUTINE cloudth_v6
1497
1498
1499
1500
1501
1502!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1503SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points,                        &
1504&           temp,qt,qt_th,frac_th,zpspsk,paprsup, paprsdn,play,thetal_th,   &
1505&           ratqs,qcloud,qincloud,icefrac,ctot,ctot_vol,                    &
1506&           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1507!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1508! Author : Etienne Vignon (LMDZ/CNRS)
1509! Date: April 2024
1510! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
1511! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
1512!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1513
1514      use lmdz_cloudth_ini, only: iflag_cloudth_vert,iflag_ratqs
1515      use lmdz_cloudth_ini, only: C_mpc ,Ni,C_cap,Ei,d_top ,vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin ,iflag_cloudth_vert_noratqs
1516      use lmdz_lscp_tools,  only: CALC_QSAT_ECMWF
1517      use lmdz_lscp_ini,    only: RTT, RG, RPI, RD, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th
1518
1519      IMPLICIT NONE
1520
1521
1522!------------------------------------------------------------------------------
1523! Declaration
1524!------------------------------------------------------------------------------
1525
1526! INPUT/OUTPUT
1527
1528      INTEGER, INTENT(IN)                         :: klon,klev,ind2
1529     
1530
1531      REAL, DIMENSION(klon),      INTENT(IN)      ::  temp          ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip
1532      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt            ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip
1533      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt_th         ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip
1534      REAL, DIMENSION(klon),      INTENT(IN)      ::  thetal_th     ! Liquid potential temperature in thermals [K]: has not seen the evap of precip
1535      REAL, DIMENSION(klon),      INTENT(IN)      ::  frac_th       ! Fraction of the mesh covered by thermals [0-1]
1536      REAL, DIMENSION(klon),      INTENT(IN)      ::  zpspsk        ! Exner potential
1537      REAL, DIMENSION(klon),      INTENT(IN)      ::  paprsup       ! Pressure at top layer interface [Pa]
1538      REAL, DIMENSION(klon),      INTENT(IN)      ::  paprsdn       ! Pressure at bottom layer interface [Pa]
1539      REAL, DIMENSION(klon),      INTENT(IN)      ::  play          ! Pressure of layers [Pa]
1540      REAL, DIMENSION(klon),      INTENT(IN)      ::  ratqs         ! Parameter that determines the width of the total water distrib.
1541
1542
1543      INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points  ! grid points with convective (thermals) mixed phase clouds
1544     
1545      REAL, DIMENSION(klon),      INTENT(OUT)      ::  ctot         ! Cloud fraction [0-1]
1546      REAL, DIMENSION(klon),      INTENT(OUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
1547      REAL, DIMENSION(klon),      INTENT(OUT)      ::  qcloud       ! In cloud total water content [kg/kg]
1548      REAL, DIMENSION(klon),      INTENT(OUT)      ::  qincloud     ! In cloud condensed water content [kg/kg]
1549      REAL, DIMENSION(klon),      INTENT(OUT)      ::  icefrac      ! Fraction of ice in clouds [0-1]
1550      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sth  ! mean saturation deficit in thermals
1551      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_senv ! mean saturation deficit in environment
1552      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmath ! std of saturation deficit in thermals
1553      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmaenv ! std of saturation deficit in environment
1554
1555
1556! LOCAL VARIABLES
1557
1558      INTEGER itap,ind1,l,ig,iter,k
1559      INTEGER iflag_topthermals, niter
1560
1561      REAL qcth(klon)
1562      REAL qcenv(klon)
1563      REAL qctot(klon)
1564      REAL cth(klon)
1565      REAL cenv(klon)   
1566      REAL cth_vol(klon)
1567      REAL cenv_vol(klon)
1568      REAL qt_env(klon), thetal_env(klon)
1569      REAL icefracenv(klon), icefracth(klon)
1570      REAL sqrtpi,sqrt2,sqrt2pi
1571      REAL alth,alenv,ath,aenv
1572      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
1573      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
1574      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
1575      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
1576      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
1577      REAL zdelta,qsatbef,zcor
1578      REAL Tbefth(klon), Tbefenv(klon), zrho(klon), rhoth(klon)
1579      REAL qlbef
1580      REAL dqsatenv(klon), dqsatth(klon)
1581      REAL erf
1582      REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
1583      REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
1584      REAL rhodz(klon)
1585      REAL rho(klon)
1586      REAL dz(klon)
1587      REAL qslenv(klon), qslth(klon), qsienv(klon), qsith(klon) 
1588      REAL alenvl, aenvl
1589      REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
1590      REAL senvi, senvl, qbase, sbase, qliqth, qiceth
1591      REAL qmax, ttarget, stmp, cout, coutref, fraci
1592      REAL maxi, mini, pas
1593
1594      ! Modifty the saturation deficit PDF in thermals
1595      ! in the presence of ice crystals
1596      CHARACTER (len = 80) :: abort_message
1597      CHARACTER (len = 20) :: routname = 'cloudth_mpc'
1598
1599
1600!------------------------------------------------------------------------------
1601! Initialisation
1602!------------------------------------------------------------------------------
1603
1604
1605! Few initial checksS
1606
1607      DO k = 1,klev
1608      DO ind1 = 1, klon
1609        rhodz(ind1) = (paprsdn(ind1)-paprsup(ind1))/rg !kg/m2
1610        zrho(ind1)  = play(ind1)/temp(ind1)/rd !kg/m3
1611        dz(ind1)    = rhodz(ind1)/zrho(ind1) !m : epaisseur de la couche en metre
1612      END DO
1613      END DO
1614
1615
1616      icefracth(:)=0.
1617      icefracenv(:)=0.
1618      sqrt2pi=sqrt(2.*rpi)
1619      sqrt2=sqrt(2.)
1620      sqrtpi=sqrt(rpi)
1621      icefrac(:)=0.
1622
1623!-------------------------------------------------------------------------------
1624! Identify grid points with potential mixed-phase conditions
1625!-------------------------------------------------------------------------------
1626
1627
1628      DO ind1=1,klon
1629            IF ((temp(ind1) .LT. RTT) .AND. (temp(ind1) .GT. temp_nowater) &
1630            .AND. (ind2<=klev-2)  &
1631            .AND. (frac_th(ind1).GT.min_frac_th_cld)) THEN
1632                mpc_bl_points(ind1,ind2)=1
1633            ELSE
1634                mpc_bl_points(ind1,ind2)=0
1635            ENDIF
1636      ENDDO
1637
1638
1639!-------------------------------------------------------------------------------
1640! Thermal fraction calculation and standard deviation of the distribution
1641!------------------------------------------------------------------------------- 
1642
1643! initialisations and calculation of temperature, humidity and saturation specific humidity
1644
1645DO ind1=1,klon
1646 
1647   rhodz(ind1)   = (paprsdn(ind1)-paprsup(ind1))/rg  ! kg/m2
1648   rho(ind1)     = play(ind1)/temp(ind1)/rd          ! kg/m3
1649   dz(ind1)      = rhodz(ind1)/rho(ind1)             ! m : epaisseur de la couche en metre
1650   Tbefenv(ind1) = temp(ind1)
1651   thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1)
1652   Tbefth(ind1)  = thetal_th(ind1)*zpspsk(ind1)
1653   rhoth(ind1)   = play(ind1)/Tbefth(ind1)/RD
1654   qt_env(ind1)  = (qt(ind1)-frac_th(ind1)*qt_th(ind1))/(1.-frac_th(ind1)) !qt = a*qtth + (1-a)*qtenv
1655
1656ENDDO
1657
1658! Calculation of saturation specific humidity
1659CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,1,.false.,qslenv,dqsatenv)
1660CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,2,.false.,qsienv,dqsatenv)
1661CALL CALC_QSAT_ECMWF(klon,Tbefth,qt_th,play,RTT,1,.false.,qslth,dqsatth)
1662
1663
1664DO ind1=1,klon
1665
1666
1667    IF (frac_th(ind1).GT.min_frac_th_cld) THEN !Thermal and environnement
1668
1669! unlike in the other cloudth routine,
1670! We consider distributions of the saturation deficit WRT liquid
1671! at positive AND negative celcius temperature
1672! subsequently, cloud fraction corresponds to the part of the pdf corresponding
1673! to superstauration with respect to liquid WHATEVER the temperature
1674
1675! Environment:
1676
1677
1678        alenv=(0.622*RLVTT*qslenv(ind1))/(rd*thetal_env(ind1)**2)     
1679        aenv=1./(1.+(alenv*RLVTT/rcpd))                             
1680        senv=aenv*(qt_env(ind1)-qslenv(ind1))                           
1681     
1682
1683! Thermals:
1684
1685
1686        alth=(0.622*RLVTT*qslth(ind1))/(rd*thetal_th(ind1)**2)       
1687        ath=1./(1.+(alth*RLVTT/rcpd))                                                         
1688        sth=ath*(qt_th(ind1)-qslth(ind1))                     
1689
1690
1691!       IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
1692
1693
1694! Standard deviation of the distributions
1695
1696           sigma1s_fraca = (sigma1s_factor**0.5)*(frac_th(ind1)**sigma1s_power) / &
1697           &                (1-frac_th(ind1))*((sth-senv)**2)**0.5
1698
1699           IF (cloudth_ratqsmin>0.) THEN
1700             sigma1s_ratqs = cloudth_ratqsmin*qt(ind1)
1701           ELSE
1702             sigma1s_ratqs = ratqs(ind1)*qt(ind1)
1703           ENDIF
1704 
1705           sigma1s = sigma1s_fraca + sigma1s_ratqs
1706           IF (iflag_ratqs.eq.11) then
1707              sigma1s = ratqs(ind1)*qt(ind1)*aenv
1708           ENDIF
1709           sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((frac_th(ind1)+0.02)**sigma2s_power))+0.002*qt_th(ind1)
1710
1711
1712           deltasenv=aenv*vert_alpha*sigma1s
1713           deltasth=ath*vert_alpha_th*sigma2s
1714
1715           xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1716           xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1717           exp_xenv1 = exp(-1.*xenv1**2)
1718           exp_xenv2 = exp(-1.*xenv2**2)
1719           xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1720           xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1721           exp_xth1 = exp(-1.*xth1**2)
1722           exp_xth2 = exp(-1.*xth2**2)
1723     
1724!surface CF
1725
1726           cth(ind1)=0.5*(1.-1.*erf(xth1))
1727           cenv(ind1)=0.5*(1.-1.*erf(xenv1))
1728           ctot(ind1)=frac_th(ind1)*cth(ind1)+(1.-1.*frac_th(ind1))*cenv(ind1)
1729
1730
1731!volume CF, condensed water and ice fraction
1732
1733            !environnement
1734
1735            IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1736            IntJ_CF=0.5*(1.-1.*erf(xenv2))
1737
1738            IF (deltasenv .LT. 1.e-10) THEN
1739              qcenv(ind1)=IntJ
1740              cenv_vol(ind1)=IntJ_CF
1741            ELSE
1742              IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
1743              IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
1744              IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
1745              IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
1746              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
1747              qcenv(ind1)=IntJ+IntI1+IntI2+IntI3
1748              cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
1749              IF (Tbefenv(ind1) .LT. temp_nowater) THEN
1750                  ! freeze all droplets in cirrus temperature regime
1751                  icefracenv(ind1)=1.
1752              ENDIF
1753            ENDIF
1754             
1755
1756
1757            !thermals
1758
1759            IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1760            IntJ_CF=0.5*(1.-1.*erf(xth2))
1761     
1762            IF (deltasth .LT. 1.e-10) THEN
1763              qcth(ind1)=IntJ
1764              cth_vol(ind1)=IntJ_CF
1765            ELSE
1766              IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1767              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1768              IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
1769              IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
1770              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
1771              qcth(ind1)=IntJ+IntI1+IntI2+IntI3
1772              cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
1773              IF (Tbefth(ind1) .LT. temp_nowater) THEN
1774                  ! freeze all droplets in cirrus temperature regime
1775                  icefracth(ind1)=1.
1776              ENDIF
1777
1778            ENDIF
1779
1780            qctot(ind1)=frac_th(ind1)*qcth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1)
1781            ctot_vol(ind1)=frac_th(ind1)*cth_vol(ind1)+(1.-1.*frac_th(ind1))*cenv_vol(ind1)
1782
1783            IF (cenv(ind1).LT.min_neb_th.and.cth(ind1).LT.min_neb_th) THEN
1784                ctot(ind1)=0.
1785                ctot_vol(ind1)=0.
1786                qcloud(ind1)=qslenv(ind1)
1787                qincloud(ind1)=0.
1788                icefrac(ind1)=0.
1789            ELSE               
1790                qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
1791                qincloud(ind1)=qctot(ind1)/ctot(ind1)
1792                IF (qctot(ind1) .GT. 0) THEN
1793                  icefrac(ind1)=(frac_th(ind1)*qcth(ind1)*icefracth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1)*icefracenv(ind1))/qctot(ind1)
1794                  icefrac(ind1)=max(min(1.,icefrac(ind1)), 0.)
1795                ENDIF
1796            ENDIF
1797
1798
1799!        ELSE ! mpc_bl_points>0
1800
1801    ELSE  ! gaussian for environment only
1802
1803     
1804        alenv=(0.622*RLVTT*qslenv(ind1))/(rd*thetal_env(ind1)**2)
1805        aenv=1./(1.+(alenv*RLVTT/rcpd))
1806        senv=aenv*(qt_env(ind1)-qslenv(ind1))
1807        sth=0.
1808        sigma1s=ratqs(ind1)*qt_env(ind1)
1809        sigma2s=0.
1810
1811        xenv=senv/(sqrt(2.)*sigma1s)
1812        cenv(ind1)=0.5*(1.-1.*erf(xenv))
1813        ctot(ind1)=0.5*(1.+1.*erf(xenv))
1814        ctot_vol(ind1)=ctot(ind1)
1815        qctot(ind1)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1))
1816
1817
1818        IF (ctot(ind1).LT.min_neb_th) THEN
1819          ctot(ind1)=0.
1820          qcloud(ind1)=qslenv(ind1)
1821          qincloud(ind1)=0.
1822        ELSE
1823          qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
1824          qincloud(ind1)=MAX(qctot(ind1)/ctot(ind1),0.)
1825        ENDIF
1826 
1827
1828    ENDIF       ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492
1829
1830    ! Outputs used to check the PDFs
1831    cloudth_senv(ind1) = senv
1832    cloudth_sth(ind1) = sth
1833    cloudth_sigmaenv(ind1) = sigma1s
1834    cloudth_sigmath(ind1) = sigma2s
1835
1836
1837    ENDDO       !loop on klon
1838
1839
1840RETURN
1841
1842
1843END SUBROUTINE cloudth_mpc
1844
1845
1846
1847!        ELSE ! mpc_bl_points>0
1848!
1849!            ! Treat boundary layer mixed phase clouds
1850!           
1851!            ! thermals
1852!            !=========
1853!
1854!           ! ice phase
1855!           !...........
1856!
1857!            qiceth=0.
1858!            deltazlev_mpc=dz(ind1,:)
1859!            temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
1860!            pres_mpc=pplay(ind1,:)
1861!            fraca_mpc=fraca(ind1,:)
1862!            snowf_mpc=snowflux(ind1,:)
1863!            iflag_topthermals=0
1864!            IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0))  THEN
1865!                iflag_topthermals = 1
1866!            ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) &
1867!                    .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
1868!                iflag_topthermals = 2
1869!            ELSE
1870!                iflag_topthermals = 0
1871!            ENDIF
1872!
1873!            CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
1874!                                   qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
1875!
1876!            ! qmax calculation
1877!            sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
1878!            deltasth=athl*vert_alpha_th*sigma2s
1879!            xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
1880!            xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
1881!            exp_xth1 = exp(-1.*xth1**2)
1882!            exp_xth2 = exp(-1.*xth2**2)
1883!            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1884!            IntJ_CF=0.5*(1.-1.*erf(xth2))
1885!            IntI1=(((sthl+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1886!            IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1887!            IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
1888!            IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
1889!            IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
1890!            qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
1891!
1892!
1893!            ! Liquid phase
1894!            !................
1895!            ! We account for the effect of ice crystals in thermals on sthl
1896!            ! and on the width of the distribution
1897!
1898!            sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2))  &
1899!                + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1)) 
1900!
1901!            sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5) &
1902!                 /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
1903!                 +0.002*zqta(ind1,ind2)
1904!            deltasthc=athl*vert_alpha_th*sigma2sc
1905!     
1906!           
1907!            xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
1908!            xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)           
1909!            exp_xth1 = exp(-1.*xth1**2)
1910!            exp_xth2 = exp(-1.*xth2**2)
1911!            IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
1912!            IntJ_CF=0.5*(1.-1.*erf(xth2))
1913!            IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1))
1914!            IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1915!            IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2)
1916!            IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc)
1917!            IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc)
1918!            qliqth=IntJ+IntI1+IntI2+IntI3
1919!           
1920!            qlth(ind1,ind2)=MAX(0.,qliqth)
1921!
1922!            ! Condensed water
1923!           
1924!            qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
1925!
1926!
1927!            ! consistency with subgrid distribution
1928!           
1929!             IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
1930!                 fraci=qith(ind1,ind2)/qcth(ind1,ind2)
1931!                 qcth(ind1,ind2)=qmax
1932!                 qith(ind1,ind2)=fraci*qmax
1933!                 qlth(ind1,ind2)=(1.-fraci)*qmax
1934!             ENDIF
1935!
1936!            ! Cloud Fraction
1937!            !...............
1938!            ! calculation of qbase which is the value of the water vapor within mixed phase clouds
1939!            ! such that the total water in cloud = qbase+qliqth+qiceth
1940!            ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
1941!            ! sbase and qbase calculation (note that sbase is wrt liq so negative)
1942!            ! look for an approximate solution with iteration
1943!           
1944!            ttarget=qcth(ind1,ind2)
1945!            mini= athl*(qsith(ind1,ind2)-qslth(ind1))
1946!            maxi= 0. !athl*(3.*sqrt(sigma2s))
1947!            niter=20
1948!            pas=(maxi-mini)/niter
1949!            stmp=mini
1950!            sbase=stmp
1951!            coutref=1.E6
1952!            DO iter=1,niter
1953!                cout=ABS(sigma2s/SQRT(2.*RPI)*EXP(-((sthl-stmp)/sigma2s)**2)+(sthl-stmp)/SQRT(2.)*(1.-erf(-(sthl-stmp)/sigma2s)) &
1954!                     + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
1955!               IF (cout .LT. coutref) THEN
1956!                     sbase=stmp
1957!                     coutref=cout
1958!                ELSE
1959!                     stmp=stmp+pas
1960!                ENDIF
1961!            ENDDO
1962!            qbase=MAX(0., sbase/athl+qslth(ind1))
1963!
1964!            ! surface cloud fraction in thermals
1965!            cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
1966!            cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
1967!
1968!
1969!            !volume cloud fraction in thermals
1970!            !to be checked
1971!            xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
1972!            xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)           
1973!            exp_xth1 = exp(-1.*xth1**2)
1974!            exp_xth2 = exp(-1.*xth2**2)
1975!
1976!            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1977!            IntJ_CF=0.5*(1.-1.*erf(xth2))
1978!     
1979!            IF (deltasth .LT. 1.e-10) THEN
1980!              cth_vol(ind1,ind2)=IntJ_CF
1981!            ELSE
1982!              IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1983!              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1984!              IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
1985!              IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
1986!              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
1987!              cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1988!            ENDIF
1989!              cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
1990!
1991!
1992!
1993!            ! Environment
1994!            !=============
1995!            ! In the environment/downdrafts, only liquid clouds
1996!            ! See Shupe et al. 2008, JAS
1997!
1998!            ! standard deviation of the distribution in the environment
1999!            sth=sthl
2000!            senv=senvl
2001!            sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
2002!                          &                (1-fraca(ind1,ind2))*(MAX((sth-senv),0.)**2)**0.5
2003!            ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
2004!            ! in the environement
2005!
2006!            sigma1s_ratqs=1E-10
2007!            IF (cloudth_ratqsmin>0.) THEN
2008!                sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
2009!            ENDIF
2010!
2011!            sigma1s = sigma1s_fraca + sigma1s_ratqs
2012!            IF (iflag_ratqs.eq.11) then
2013!               sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
2014!            ENDIF
2015!            IF (iflag_ratqs.eq.11) then
2016!              sigma1s = ratqs(ind1,ind2)*po(ind1)*aenvl
2017!            ENDIF
2018!            deltasenv=aenvl*vert_alpha*sigma1s
2019!            xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
2020!            xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
2021!            exp_xenv1 = exp(-1.*xenv1**2)
2022!            exp_xenv2 = exp(-1.*xenv2**2)
2023!
2024!            !surface CF
2025!            cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
2026!
2027!            !volume CF and condensed water
2028!            IntJ=0.5*senvl*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
2029!            IntJ_CF=0.5*(1.-1.*erf(xenv2))
2030!
2031!            IF (deltasenv .LT. 1.e-10) THEN
2032!              qcenv(ind1,ind2)=IntJ
2033!              cenv_vol(ind1,ind2)=IntJ_CF
2034!            ELSE
2035!              IntI1=(((senvl+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
2036!              IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
2037!              IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
2038!              IntI1_CF=((senvl+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
2039!              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
2040!              qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
2041!              cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
2042!            ENDIF
2043!
2044!            qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
2045!            cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
2046!
2047!
2048!           
2049!            ! Thermals + environment
2050!            ! =======================
2051!            ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
2052!            qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
2053!            ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
2054!            IF (qcth(ind1,ind2) .GT. 0) THEN
2055!               icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
2056!                    /(fraca(ind1,ind2)*qcth(ind1,ind2) &
2057!                    +(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2))
2058!                icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
2059!            ELSE
2060!                icefrac(ind1,ind2)=0.
2061!            ENDIF
2062!
2063!            IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
2064!                ctot(ind1,ind2)=0.
2065!                ctot_vol(ind1,ind2)=0.
2066!                qincloud(ind1)=0.
2067!                qcloud(ind1)=zqsatenv(ind1)
2068!            ELSE               
2069!                qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
2070!                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
2071!                qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
2072!                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
2073!            ENDIF
2074!
2075!        ENDIF ! mpc_bl_points
2076
2077
2078
2079!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2080SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp,pres,qth,qsith,qlth,deltazlev,vith,fraca,qith)
2081
2082! parameterization of ice for boundary
2083! layer mixed-phase clouds assuming a stationary system
2084!
2085! Note that vapor deposition on ice crystals and riming of liquid droplets
2086! depend on the ice number concentration Ni
2087! One could assume that Ni depends on qi, e.g.,  Ni=beta*(qi*rho)**xi
2088! and use values from Hong et al. 2004, MWR for instance
2089! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962
2090! One could also think of a more complex expression of Ni;
2091! function of qi, T, the concentration in aerosols or INP ..
2092! Here we prefer fixing Ni to a tuning parameter
2093! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
2094! in Mioche et al. 2017
2095!
2096!
2097! References:
2098!------------
2099! This parameterization is thoroughly described in Vignon et al.
2100!
2101! More specifically
2102! for the Water vapor deposition process:
2103!
2104! Rotstayn, L. D., 1997: A physically based scheme for the treat-
2105! ment of stratiform cloudfs and precipitation in large-scale
2106! models. I: Description and evaluation of the microphysical
2107! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.
2108!
2109!  Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
2110!  stratiform cloud microphysics scheme in the NCAR Com-
2111!  munity Atmosphere Model (CAM3). Part I: Description and
2112!  numerical tests. J. Climate, 21, 3642–3659
2113!
2114! for the Riming process:
2115!
2116! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
2117! scale structure and organization of clouds and precipitation in
2118! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
2119! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206
2120!
2121! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
2122! forecasts of winter precipitation using an improved bulk
2123! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
2124! forecasts of winter precipitation using an improved bulk
2125! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542
2126!
2127! For the formation of clouds by thermals:
2128!
2129! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
2130! the Atmospheric Sciences, 65, 407–425.
2131!
2132! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
2133! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3
2134!
2135!
2136!
2137! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
2138!=============================================================================
2139
2140    use phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm
2141
2142    USE yomcst_mod_h
2143IMPLICIT none
2144
2145
2146
2147    INTEGER, INTENT(IN) :: ind1,ind2, klev           ! horizontal and vertical indices and dimensions
2148    INTEGER, INTENT(IN) :: iflag_topthermals         ! uppermost layer of thermals ?
2149    REAL, INTENT(IN)    :: Ni                        ! ice number concentration [m-3]
2150    REAL, INTENT(IN)    :: Ei                        ! ice-droplet collision efficiency
2151    REAL, INTENT(IN)    :: C_cap                     ! ice crystal capacitance
2152    REAL, INTENT(IN)    :: d_top                     ! cloud-top ice crystal mixing parameter
2153    REAL,  DIMENSION(klev), INTENT(IN) :: temp       ! temperature [K] within thermals
2154    REAL,  DIMENSION(klev), INTENT(IN) :: pres       ! pressure [Pa]
2155    REAL,  DIMENSION(klev), INTENT(IN) :: qth        ! mean specific water content in thermals [kg/kg]
2156    REAL,  DIMENSION(klev), INTENT(IN) :: qsith        ! saturation specific humidity wrt ice in thermals [kg/kg]
2157    REAL,  DIMENSION(klev), INTENT(IN) :: qlth       ! condensed liquid water in thermals, approximated value [kg/kg]
2158    REAL,  DIMENSION(klev), INTENT(IN) :: deltazlev  ! layer thickness [m]
2159    REAL,  DIMENSION(klev), INTENT(IN) :: vith       ! ice crystal fall velocity [m/s]
2160    REAL,  DIMENSION(klev+1), INTENT(IN) :: fraca      ! fraction of the mesh covered by thermals
2161    REAL,  DIMENSION(klev), INTENT(INOUT) :: qith       ! condensed ice water , thermals [kg/kg]
2162
2163
2164    INTEGER ind2p1,ind2p2
2165    REAL rho(klev)
2166    REAL unsurtaudet, unsurtaustardep, unsurtaurim
2167    REAL qi, AA, BB, Ka, Dv, rhoi
2168    REAL p0, t0, fp1, fp2
2169    REAL alpha, flux_term
2170    REAL det_term, precip_term, rim_term, dep_term
2171
2172
2173    ind2p1=ind2+1
2174    ind2p2=ind2+2
2175
2176    rho=pres/temp/RD  ! air density kg/m3
2177
2178    Ka=2.4e-2      ! thermal conductivity of the air, SI
2179    p0=101325.0    ! ref pressure
2180    T0=273.15      ! ref temp
2181    rhoi=500.0     ! cloud ice density following Reisner et al. 1998
2182    alpha=700.     ! fallvelocity param
2183
2184
2185    IF (iflag_topthermals .GT. 0) THEN ! uppermost thermals levels
2186
2187    Dv=0.0001*0.211*(p0/pres(ind2))*((temp(ind2)/T0)**1.94) ! water vapor diffusivity in air, SI
2188
2189    ! Detrainment term:
2190
2191    unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2)
2192 
2193
2194    ! Deposition term
2195    AA=RLSTT/Ka/temp(ind2)*(RLSTT/RV/temp(ind2)-1.)
2196    BB=1./(rho(ind2)*Dv*qsith(ind2))
2197    unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2)-qsith(ind2))/qsith(ind2) &
2198                    *4.*RPI/(AA+BB)*(6.*rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33)
2199
2200    ! Riming term  neglected at this level
2201    !unsurtaurim=rho(ind2)*alpha*3./rhoi/2.*Ei*qlth(ind2)*((p0/pres(ind2))**0.4)
2202
2203    qi=fraca(ind2)*unsurtaustardep/MAX((d_top*unsurtaudet),1E-12)
2204    qi=MAX(qi,0.)**(3./2.)
2205
2206    ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2
2207
2208    Dv=0.0001*0.211*(p0/pres(ind2p1))*((temp(ind2p1)/T0)**1.94) ! water vapor diffusivity in air, SI
2209
2210    ! Detrainment term:
2211
2212    unsurtaudet=detr_therm(ind1,ind2p1)/rho(ind2p1)/deltazlev(ind2p1)
2213    det_term=-unsurtaudet*qith(ind2p1)*rho(ind2p1)
2214   
2215   
2216    ! Deposition term
2217
2218    AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.)
2219    BB=1./(rho(ind2p1)*Dv*qsith(ind2p1))
2220    unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2p1)-qsith(ind2p1)) &
2221         /qsith(ind2p1)*4.*RPI/(AA+BB) &
2222         *(6.*rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33)
2223    dep_term=rho(ind2p1)*fraca(ind2p1)*(qith(ind2p1)**0.33)*unsurtaustardep
2224 
2225    ! Riming term
2226
2227    unsurtaurim=rho(ind2p1)*alpha*3./rhoi/2.*Ei*qlth(ind2p1)*((p0/pres(ind2p1))**0.4)
2228    rim_term=rho(ind2p1)*fraca(ind2p1)*qith(ind2p1)*unsurtaurim
2229
2230    ! Precip term
2231
2232    ! We assume that there is no solid precipitation outside thermals
2233    ! so the precipitation flux within thermals is equal to the precipitation flux
2234    ! at mesh-scale divided by  thermals fraction
2235   
2236    fp2=0.
2237    fp1=0.
2238    IF (fraca(ind2p1) .GT. 0.) THEN
2239    fp2=-qith(ind2p2)*rho(ind2p2)*vith(ind2p2)*fraca(ind2p2)! flux defined positive upward
2240    fp1=-qith(ind2p1)*rho(ind2p1)*vith(ind2p1)*fraca(ind2p1)
2241    ENDIF
2242
2243    precip_term=-1./deltazlev(ind2p1)*(fp2-fp1)
2244
2245    ! Calculation in a top-to-bottom loop
2246
2247    IF (fm_therm(ind1,ind2p1) .GT. 0.) THEN
2248          qi= 1./fm_therm(ind1,ind2p1)* &
2249              (deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + &
2250              fm_therm(ind1,ind2p2)*(qith(ind2p1)))
2251    END IF
2252
2253    ENDIF ! top thermals
2254
2255    qith(ind2)=MAX(0.,qi)
2256
2257    RETURN
2258
2259END SUBROUTINE ICE_MPC_BL_CLOUDS
2260
2261END MODULE lmdz_cloudth
Note: See TracBrowser for help on using the repository browser.