source: LMDZ6/trunk/libf/phylmd/cloudth_mod.F90 @ 4564

Last change on this file since 4564 was 4535, checked in by evignon, 18 months ago

poursuite de la replay-isation de lscp en vue de la session
de reecriture de lscp_mod en juin

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