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

Last change on this file since 3709 was 3570, checked in by fhourdin, 5 years ago

Sortie de nouveaux parametres pour controler la distribution des nuages.
File size: 57.1 KB
RevLine 
[2686]1MODULE cloudth_mod
2
3  IMPLICIT NONE
4
5CONTAINS
6
7       SUBROUTINE cloudth(ngrid,klev,ind2,  &
8     &           ztv,po,zqta,fraca, &
[3493]9     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
[2686]10     &           ratqs,zqs,t)
11
12
13      IMPLICIT NONE
14
15
16!===========================================================================
17! Auteur : Arnaud Octavio Jam (LMD/CNRS)
18! Date : 25 Mai 2010
19! Objet : calcule les valeurs de qc et rneb dans les thermiques
20!===========================================================================
21
22
23#include "YOMCST.h"
24#include "YOETHF.h"
25#include "FCTTRE.h"
26#include "thermcell.h"
27#include "nuage.h"
28
29      INTEGER itap,ind1,ind2
30      INTEGER ngrid,klev,klon,l,ig
31     
32      REAL ztv(ngrid,klev)
33      REAL po(ngrid)
34      REAL zqenv(ngrid)   
35      REAL zqta(ngrid,klev)
36         
37      REAL fraca(ngrid,klev+1)
38      REAL zpspsk(ngrid,klev)
39      REAL paprs(ngrid,klev+1)
[3493]40      REAL pplay(ngrid,klev)
[2686]41      REAL ztla(ngrid,klev)
42      REAL zthl(ngrid,klev)
43
44      REAL zqsatth(ngrid,klev)
45      REAL zqsatenv(ngrid,klev)
46     
47     
48      REAL sigma1(ngrid,klev)
49      REAL sigma2(ngrid,klev)
50      REAL qlth(ngrid,klev)
51      REAL qlenv(ngrid,klev)
52      REAL qltot(ngrid,klev)
53      REAL cth(ngrid,klev) 
54      REAL cenv(ngrid,klev)   
55      REAL ctot(ngrid,klev)
56      REAL rneb(ngrid,klev)
57      REAL t(ngrid,klev)
58      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
59      REAL rdd,cppd,Lv
60      REAL alth,alenv,ath,aenv
61      REAL sth,senv,sigma1s,sigma2s,xth,xenv
62      REAL Tbef,zdelta,qsatbef,zcor
63      REAL qlbef 
64      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
65     
66      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
67      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
68      REAL zqs(ngrid), qcloud(ngrid)
69      REAL erf
70
71
72
73
74!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
75! Gestion de deux versions de cloudth
76!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
77
78      IF (iflag_cloudth_vert.GE.1) THEN
79      CALL cloudth_vert(ngrid,klev,ind2,  &
80     &           ztv,po,zqta,fraca, &
[3493]81     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
[2686]82     &           ratqs,zqs,t)
83      RETURN
84      ENDIF
85!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
86
87
88!-------------------------------------------------------------------------------
89! Initialisation des variables r?elles
90!-------------------------------------------------------------------------------
91      sigma1(:,:)=0.
92      sigma2(:,:)=0.
93      qlth(:,:)=0.
94      qlenv(:,:)=0. 
95      qltot(:,:)=0.
96      rneb(:,:)=0.
97      qcloud(:)=0.
98      cth(:,:)=0.
99      cenv(:,:)=0.
100      ctot(:,:)=0.
101      qsatmmussig1=0.
102      qsatmmussig2=0.
103      rdd=287.04
104      cppd=1005.7
105      pi=3.14159
106      Lv=2.5e6
107      sqrt2pi=sqrt(2.*pi)
108
109
110
111!-------------------------------------------------------------------------------
112! Calcul de la fraction du thermique et des ?cart-types des distributions
113!-------------------------------------------------------------------------------                 
114      do ind1=1,ngrid
115
116      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
117
118      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
119
120
121!      zqenv(ind1)=po(ind1)
122      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
123      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
124      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
125      qsatbef=MIN(0.5,qsatbef)
126      zcor=1./(1.-retv*qsatbef)
127      qsatbef=qsatbef*zcor
128      zqsatenv(ind1,ind2)=qsatbef
129
130
131
132
133      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
134      aenv=1./(1.+(alenv*Lv/cppd))
135      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
136
137
138
139
140      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
141      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
142      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
143      qsatbef=MIN(0.5,qsatbef)
144      zcor=1./(1.-retv*qsatbef)
145      qsatbef=qsatbef*zcor
146      zqsatth(ind1,ind2)=qsatbef
147           
148      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
149      ath=1./(1.+(alth*Lv/cppd))
150      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
151     
152     
153
154!------------------------------------------------------------------------------
155! Calcul des ?cart-types pour s
156!------------------------------------------------------------------------------
157
158!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
159!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
160!       if (paprs(ind1,ind2).gt.90000) then
161!       ratqs(ind1,ind2)=0.002
162!       else
163!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
164!       endif
165       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
166       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
167!       sigma1s=ratqs(ind1,ind2)*po(ind1)
168!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
169 
170!------------------------------------------------------------------------------
171! Calcul de l'eau condens?e et de la couverture nuageuse
172!------------------------------------------------------------------------------
173      sqrt2pi=sqrt(2.*pi)
174      xth=sth/(sqrt(2.)*sigma2s)
175      xenv=senv/(sqrt(2.)*sigma1s)
176      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
177      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
178      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
179
180      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
181      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
182      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
183
184!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
185      if (ctot(ind1,ind2).lt.1.e-10) then
186      ctot(ind1,ind2)=0.
187      qcloud(ind1)=zqsatenv(ind1,ind2)
188
189      else   
190               
191      ctot(ind1,ind2)=ctot(ind1,ind2)
192      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
193
194      endif                           
195     
196         
197!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
198
199
200      else  ! gaussienne environnement seule
201     
202      zqenv(ind1)=po(ind1)
203      Tbef=t(ind1,ind2)
204      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
205      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
206      qsatbef=MIN(0.5,qsatbef)
207      zcor=1./(1.-retv*qsatbef)
208      qsatbef=qsatbef*zcor
209      zqsatenv(ind1,ind2)=qsatbef
210     
211
212!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
213      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
214      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
215      aenv=1./(1.+(alenv*Lv/cppd))
216      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
217     
218
219      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
220
221      sqrt2pi=sqrt(2.*pi)
222      xenv=senv/(sqrt(2.)*sigma1s)
223      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
224      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
225     
226      if (ctot(ind1,ind2).lt.1.e-3) then
227      ctot(ind1,ind2)=0.
228      qcloud(ind1)=zqsatenv(ind1,ind2)
229
230      else   
231               
232      ctot(ind1,ind2)=ctot(ind1,ind2)
233      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
234
235      endif   
236 
237 
238 
239 
240 
241 
242      endif   
243      enddo
244     
245      return
246!     end
247END SUBROUTINE cloudth
248
249
250
251!===========================================================================
252     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
253     &           ztv,po,zqta,fraca, &
[3493]254     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
[2686]255     &           ratqs,zqs,t)
256
257!===========================================================================
258! Auteur : Arnaud Octavio Jam (LMD/CNRS)
259! Date : 25 Mai 2010
260! Objet : calcule les valeurs de qc et rneb dans les thermiques
261!===========================================================================
262
263
264      USE ioipsl_getin_p_mod, ONLY : getin_p
265
266      IMPLICIT NONE
267
268#include "YOMCST.h"
269#include "YOETHF.h"
270#include "FCTTRE.h"
271#include "thermcell.h"
272#include "nuage.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)
[3493]285      REAL pplay(ngrid,klev)
[2686]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
[2909]315      !$OMP THREADPRIVATE(vert_alpha)
[2686]316      LOGICAL, SAVE :: firstcall = .TRUE.
[2909]317      !$OMP THREADPRIVATE(firstcall)
[2686]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
[3493]589
590
591
[2686]592       SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
593     &           ztv,po,zqta,fraca, &
[3493]594     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
[2686]595     &           ratqs,zqs,t)
596
597
598      IMPLICIT NONE
599
600
601!===========================================================================
602! Author : Arnaud Octavio Jam (LMD/CNRS)
603! Date : 25 Mai 2010
604! Objet : calcule les valeurs de qc et rneb dans les thermiques
605!===========================================================================
606
607
608#include "YOMCST.h"
609#include "YOETHF.h"
610#include "FCTTRE.h"
611#include "thermcell.h"
612#include "nuage.h"
613
614      INTEGER itap,ind1,ind2
615      INTEGER ngrid,klev,klon,l,ig
616     
617      REAL ztv(ngrid,klev)
618      REAL po(ngrid)
619      REAL zqenv(ngrid)   
620      REAL zqta(ngrid,klev)
621         
622      REAL fraca(ngrid,klev+1)
623      REAL zpspsk(ngrid,klev)
624      REAL paprs(ngrid,klev+1)
[3493]625      REAL pplay(ngrid,klev)
[2686]626      REAL ztla(ngrid,klev)
627      REAL zthl(ngrid,klev)
628
629      REAL zqsatth(ngrid,klev)
630      REAL zqsatenv(ngrid,klev)
631     
[2945]632      REAL sigma1(ngrid,klev)                                                         
[2686]633      REAL sigma2(ngrid,klev)
634      REAL qlth(ngrid,klev)
635      REAL qlenv(ngrid,klev)
636      REAL qltot(ngrid,klev)
[2945]637      REAL cth(ngrid,klev)
[2686]638      REAL cenv(ngrid,klev)   
639      REAL ctot(ngrid,klev)
[2945]640      REAL cth_vol(ngrid,klev)
641      REAL cenv_vol(ngrid,klev)
642      REAL ctot_vol(ngrid,klev)
643      REAL rneb(ngrid,klev)     
[2686]644      REAL t(ngrid,klev)
645      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
646      REAL rdd,cppd,Lv
647      REAL alth,alenv,ath,aenv
648      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
[3493]649      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
[2686]650      REAL Tbef,zdelta,qsatbef,zcor
651      REAL qlbef 
652      REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
653      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
654      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
655      REAL zqs(ngrid), qcloud(ngrid)
656      REAL erf
657
658
659
660      IF (iflag_cloudth_vert.GE.1) THEN
661      CALL cloudth_vert_v3(ngrid,klev,ind2,  &
662     &           ztv,po,zqta,fraca, &
[3493]663     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
[2686]664     &           ratqs,zqs,t)
665      RETURN
666      ENDIF
667!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
668
669
670!-------------------------------------------------------------------------------
671! Initialisation des variables r?elles
672!-------------------------------------------------------------------------------
673      sigma1(:,:)=0.
674      sigma2(:,:)=0.
675      qlth(:,:)=0.
676      qlenv(:,:)=0. 
677      qltot(:,:)=0.
678      rneb(:,:)=0.
679      qcloud(:)=0.
680      cth(:,:)=0.
681      cenv(:,:)=0.
682      ctot(:,:)=0.
[2945]683      cth_vol(:,:)=0.
684      cenv_vol(:,:)=0.
685      ctot_vol(:,:)=0.
[2686]686      qsatmmussig1=0.
687      qsatmmussig2=0.
688      rdd=287.04
689      cppd=1005.7
690      pi=3.14159
691      Lv=2.5e6
692      sqrt2pi=sqrt(2.*pi)
693      sqrt2=sqrt(2.)
694      sqrtpi=sqrt(pi)
695
696
697!-------------------------------------------------------------------------------
698! Cloud fraction in the thermals and standard deviation of the PDFs
699!-------------------------------------------------------------------------------                 
700      do ind1=1,ngrid
701
702      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
703
704      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
705
706      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
707      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
708      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
709      qsatbef=MIN(0.5,qsatbef)
710      zcor=1./(1.-retv*qsatbef)
711      qsatbef=qsatbef*zcor
712      zqsatenv(ind1,ind2)=qsatbef
713
714
715      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
716      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
717      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
718
719!po = qt de l'environnement ET des thermique
720!zqenv = qt environnement
721!zqsatenv = qsat environnement
722!zthl = Tl environnement
723
724
725      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
726      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
727      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
728      qsatbef=MIN(0.5,qsatbef)
729      zcor=1./(1.-retv*qsatbef)
730      qsatbef=qsatbef*zcor
731      zqsatth(ind1,ind2)=qsatbef
732           
733      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
734      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
735      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
736     
737!zqta = qt thermals
738!zqsatth = qsat thermals
739!ztla = Tl thermals
740
741!------------------------------------------------------------------------------
742! s standard deviations
743!------------------------------------------------------------------------------
744
745!     tests
746!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
747!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
748!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
749!     final option
750      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
751      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
752 
753!------------------------------------------------------------------------------
754! Condensed water and cloud cover
755!------------------------------------------------------------------------------
756      xth=sth/(sqrt2*sigma2s)
757      xenv=senv/(sqrt2*sigma1s)
758      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
759      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
760      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
[2945]761      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[2686]762
763      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
764      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
765      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
766
767      if (ctot(ind1,ind2).lt.1.e-10) then
768      ctot(ind1,ind2)=0.
769      qcloud(ind1)=zqsatenv(ind1,ind2)
770      else
771      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
772      endif
773
774      else  ! Environnement only, follow the if l.110
775     
776      zqenv(ind1)=po(ind1)
777      Tbef=t(ind1,ind2)
778      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
779      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
780      qsatbef=MIN(0.5,qsatbef)
781      zcor=1./(1.-retv*qsatbef)
782      qsatbef=qsatbef*zcor
783      zqsatenv(ind1,ind2)=qsatbef
784
785!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
786      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
787      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
788      aenv=1./(1.+(alenv*Lv/cppd))
789      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
790     
791      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
792
793      xenv=senv/(sqrt2*sigma1s)
794      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
[2945]795      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[2686]796      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
797
798      if (ctot(ind1,ind2).lt.1.e-3) then
799      ctot(ind1,ind2)=0.
800      qcloud(ind1)=zqsatenv(ind1,ind2)
801      else   
802      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
803      endif
804
805
806      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
807      enddo       ! from the loop on ngrid l.108
808      return
809!     end
810END SUBROUTINE cloudth_v3
811
812
813
814!===========================================================================
815     SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
816     &           ztv,po,zqta,fraca, &
[3493]817     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
[2686]818     &           ratqs,zqs,t)
819
820!===========================================================================
821! Auteur : Arnaud Octavio Jam (LMD/CNRS)
822! Date : 25 Mai 2010
823! Objet : calcule les valeurs de qc et rneb dans les thermiques
824!===========================================================================
825
826
827      USE ioipsl_getin_p_mod, ONLY : getin_p
[2958]828      USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
829     &                                cloudth_sigmath,cloudth_sigmaenv
[2686]830
831      IMPLICIT NONE
832
833#include "YOMCST.h"
834#include "YOETHF.h"
835#include "FCTTRE.h"
836#include "thermcell.h"
837#include "nuage.h"
838     
839      INTEGER itap,ind1,ind2
840      INTEGER ngrid,klev,klon,l,ig
841     
842      REAL ztv(ngrid,klev)
843      REAL po(ngrid)
844      REAL zqenv(ngrid)   
845      REAL zqta(ngrid,klev)
846         
847      REAL fraca(ngrid,klev+1)
848      REAL zpspsk(ngrid,klev)
849      REAL paprs(ngrid,klev+1)
[3493]850      REAL pplay(ngrid,klev)
[2686]851      REAL ztla(ngrid,klev)
852      REAL zthl(ngrid,klev)
853
854      REAL zqsatth(ngrid,klev)
855      REAL zqsatenv(ngrid,klev)
856     
857      REAL sigma1(ngrid,klev)                                                         
858      REAL sigma2(ngrid,klev)
859      REAL qlth(ngrid,klev)
860      REAL qlenv(ngrid,klev)
861      REAL qltot(ngrid,klev)
[2945]862      REAL cth(ngrid,klev)
[2686]863      REAL cenv(ngrid,klev)   
864      REAL ctot(ngrid,klev)
[2945]865      REAL cth_vol(ngrid,klev)
866      REAL cenv_vol(ngrid,klev)
867      REAL ctot_vol(ngrid,klev)
[2686]868      REAL rneb(ngrid,klev)
869      REAL t(ngrid,klev)                                                                 
870      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
871      REAL rdd,cppd,Lv
872      REAL alth,alenv,ath,aenv
[2957]873      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
[3493]874      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
[2957]875      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
[2686]876      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
[2945]877      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
[2686]878      REAL Tbef,zdelta,qsatbef,zcor
879      REAL qlbef 
880      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
881      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
882      ! (J Jouhaud, JL Dufresne, JB Madeleine)
[2911]883      REAL,SAVE :: vert_alpha, vert_alpha_th
884      !$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
[2957]885      REAL,SAVE :: sigma1s_factor=1.1
886      REAL,SAVE :: sigma1s_power=0.6
[3570]887      REAL,SAVE :: sigma2s_factor=0.09
888      REAL,SAVE :: sigma2s_power=0.5
[2960]889      REAL,SAVE :: cloudth_ratqsmin=-1.
[3570]890      !$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin)
[2957]891      INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
892      !$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
893
[2686]894      LOGICAL, SAVE :: firstcall = .TRUE.
[2910]895      !$OMP THREADPRIVATE(firstcall)
[2686]896
897      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
898      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
899      REAL zqs(ngrid), qcloud(ngrid)
900      REAL erf
901
[3493]902      REAL rhodz(ngrid,klev)
903      REAL zrho(ngrid,klev)
904      REAL dz(ngrid,klev)
905
906      DO ind1 = 1, ngrid
907        !Layer calculation
908        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
909        zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
910        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
911      END DO
912
913
[2686]914!------------------------------------------------------------------------------
915! Initialize
916!------------------------------------------------------------------------------
917      sigma1(:,:)=0.
918      sigma2(:,:)=0.
919      qlth(:,:)=0.
920      qlenv(:,:)=0. 
921      qltot(:,:)=0.
922      rneb(:,:)=0.
923      qcloud(:)=0.
924      cth(:,:)=0.
925      cenv(:,:)=0.
926      ctot(:,:)=0.
[2945]927      cth_vol(:,:)=0.
928      cenv_vol(:,:)=0.
929      ctot_vol(:,:)=0.
[2686]930      qsatmmussig1=0.
931      qsatmmussig2=0.
932      rdd=287.04
933      cppd=1005.7
934      pi=3.14159
935      Lv=2.5e6
936      sqrt2pi=sqrt(2.*pi)
937      sqrt2=sqrt(2.)
938      sqrtpi=sqrt(pi)
939
940      IF (firstcall) THEN
941        vert_alpha=0.5
942        CALL getin_p('cloudth_vert_alpha',vert_alpha)
943        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
[2911]944        ! The factor used for the thermal is equal to that of the environment
945        !   if nothing is explicitly specified in the def file
946        vert_alpha_th=vert_alpha
947        CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
948        WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th
[2957]949        ! Factor used in the calculation of sigma1s
950        CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor)
951        WRITE(*,*) 'cloudth_sigma1s_factor = ', sigma1s_factor
952        ! Power used in the calculation of sigma1s
953        CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
954        WRITE(*,*) 'cloudth_sigma1s_power = ', sigma1s_power
[3570]955        ! Factor used in the calculation of sigma2s
956        CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor)
957        WRITE(*,*) 'cloudth_sigma2s_factor = ', sigma2s_factor
958        ! Power used in the calculation of sigma2s
959        CALL getin_p('cloudth_sigma2s_power',sigma2s_power)
960        WRITE(*,*) 'cloudth_sigma2s_power = ', sigma2s_power
[2960]961        ! Minimum value for the environmental air subgrid water distrib
962        CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
963        WRITE(*,*) 'cloudth_ratqsmin = ', cloudth_ratqsmin
[2957]964        ! Remove the dependency to ratqs from the variance of the vertical PDF
965        CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs)
966        WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
967
[2686]968        firstcall=.FALSE.
969      ENDIF
970
971!-------------------------------------------------------------------------------
972! Calcul de la fraction du thermique et des ecart-types des distributions
973!-------------------------------------------------------------------------------                 
974      do ind1=1,ngrid
975
976      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
977
978      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
979
980
981      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
982      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
983      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
984      qsatbef=MIN(0.5,qsatbef)
985      zcor=1./(1.-retv*qsatbef)
986      qsatbef=qsatbef*zcor
987      zqsatenv(ind1,ind2)=qsatbef
988
989
990      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
991      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
992      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
993
994!zqenv = qt environnement
995!zqsatenv = qsat environnement
996!zthl = Tl environnement
997
998
999      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1000      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1001      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1002      qsatbef=MIN(0.5,qsatbef)
1003      zcor=1./(1.-retv*qsatbef)
1004      qsatbef=qsatbef*zcor
1005      zqsatth(ind1,ind2)=qsatbef
1006           
1007      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
1008      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
1009      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
1010     
1011     
1012!zqta = qt thermals
1013!zqsatth = qsat thermals
1014!ztla = Tl thermals
1015
1016!------------------------------------------------------------------------------
1017! s standard deviation
1018!------------------------------------------------------------------------------
1019
[2957]1020      sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
1021     &                (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
1022!     sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
[2960]1023      IF (cloudth_ratqsmin>0.) THEN
1024         sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
1025      ELSE
1026         sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
1027      ENDIF
[2957]1028      sigma1s = sigma1s_fraca + sigma1s_ratqs
[3570]1029      sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
[2686]1030!      tests
1031!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
1032!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
1033!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
1034!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
1035!       if (paprs(ind1,ind2).gt.90000) then
1036!       ratqs(ind1,ind2)=0.002
1037!       else
1038!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
1039!       endif
1040!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
1041!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
1042!       sigma1s=ratqs(ind1,ind2)*po(ind1)
1043!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
1044 
1045       IF (iflag_cloudth_vert == 1) THEN
1046!-------------------------------------------------------------------------------
1047!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
1048!-------------------------------------------------------------------------------
1049
1050      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1051      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1052
1053      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
1054      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
1055      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
1056      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
1057      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
1058      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
1059     
1060      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
1061      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
1062      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
1063
1064      ! Environment
1065      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
1066      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
1067      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
1068      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
1069
1070      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1071
1072      ! Thermal
1073      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
1074      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
1075      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
1076      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
1077      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1078      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1079
[2911]1080      ELSE IF (iflag_cloudth_vert >= 3) THEN
[3493]1081      IF (iflag_cloudth_vert < 5) THEN
[2686]1082!-------------------------------------------------------------------------------
1083!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1084!-------------------------------------------------------------------------------
1085!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
1086!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
1087!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1088!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
[2911]1089      IF (iflag_cloudth_vert == 3) THEN
1090        deltasenv=aenv*vert_alpha*sigma1s
1091        deltasth=ath*vert_alpha_th*sigma2s
1092      ELSE IF (iflag_cloudth_vert == 4) THEN
[2957]1093        IF (iflag_cloudth_vert_noratqs == 1) THEN
[2959]1094          deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
[2957]1095          deltasth=vert_alpha_th*sigma2s
1096        ELSE
1097          deltasenv=vert_alpha*sigma1s
1098          deltasth=vert_alpha_th*sigma2s
1099        ENDIF
[2911]1100      ENDIF
[2686]1101     
1102      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1103      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1104      exp_xenv1 = exp(-1.*xenv1**2)
1105      exp_xenv2 = exp(-1.*xenv2**2)
1106      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1107      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1108      exp_xth1 = exp(-1.*xth1**2)
1109      exp_xth2 = exp(-1.*xth2**2)
1110     
[2945]1111      !CF_surfacique
[2686]1112      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
1113      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
[2945]1114      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
[2686]1115
[2945]1116
1117      !CF_volumique & eau condense
[2686]1118      !environnement
1119      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
[2945]1120      IntJ_CF=0.5*(1.-1.*erf(xenv2))
[2686]1121      if (deltasenv .lt. 1.e-10) then
1122      qlenv(ind1,ind2)=IntJ
[2945]1123      cenv_vol(ind1,ind2)=IntJ_CF
[2686]1124      else
1125      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
1126      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
1127      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
[2945]1128      IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
1129      IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
[2686]1130      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
[2945]1131      cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
[2686]1132      endif
1133
1134      !thermique
1135      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
[2945]1136      IntJ_CF=0.5*(1.-1.*erf(xth2))
1137      if (deltasth .lt. 1.e-10) then
[2686]1138      qlth(ind1,ind2)=IntJ
[2945]1139      cth_vol(ind1,ind2)=IntJ_CF
[2686]1140      else
1141      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1142      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1143      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
[2945]1144      IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
1145      IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
[2686]1146      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
[2945]1147      cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
[2686]1148      endif
1149
1150      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
[2945]1151      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
[2686]1152
[3493]1153      ELSE IF (iflag_cloudth_vert == 5) THEN
1154      sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5)+ratqs(ind1,ind2)*po(ind1) !Environment
1155      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
1156      !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
1157      !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
1158      xth=sth/(sqrt(2.)*sigma2s)
1159      xenv=senv/(sqrt(2.)*sigma1s)
1160
1161      !Volumique
1162      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
1163      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1164      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1165      !print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
1166
1167      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2))
1168      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2)) 
1169      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1170
1171      !Surfacique
1172      !Neggers
1173      !beta=0.0044
1174      !inverse_rho=1.+beta*dz(ind1,ind2)
1175      !print *,'jeanjean : beta=',beta
1176      !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
1177      !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
1178      !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1179
1180      !Brooks
1181      a_Brooks=0.6694
1182      b_Brooks=0.1882
1183      A_Maj_Brooks=0.1635 !-- sans shear
1184      !A_Maj_Brooks=0.17   !-- ARM LES
1185      !A_Maj_Brooks=0.18   !-- RICO LES
1186      !A_Maj_Brooks=0.19   !-- BOMEX LES
1187      Dx_Brooks=200000.
1188      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
1189      !print *,'jeanjean_f=',f_Brooks
1190
1191      cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
1192      cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
1193      ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1194      !print *,'JJ_ctot_1',ctot(ind1,ind2)
1195
1196
1197
1198
1199
1200      ENDIF ! of if (iflag_cloudth_vert<5)
[2911]1201      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
[2686]1202
[3495]1203!      if (ctot(ind1,ind2).lt.1.e-10) then
1204      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
[2686]1205      ctot(ind1,ind2)=0.
[2945]1206      ctot_vol(ind1,ind2)=0.
[3493]1207      qcloud(ind1)=zqsatenv(ind1,ind2)
[2686]1208
[3493]1209      else
[2686]1210               
1211      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1212!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1213!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1214
1215      endif 
1216
[3493]1217      else  ! gaussienne environnement seule
[2686]1218     
1219      zqenv(ind1)=po(ind1)
1220      Tbef=t(ind1,ind2)
1221      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1222      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1223      qsatbef=MIN(0.5,qsatbef)
1224      zcor=1./(1.-retv*qsatbef)
1225      qsatbef=qsatbef*zcor
1226      zqsatenv(ind1,ind2)=qsatbef
1227     
1228
[3493]1229!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
[2686]1230      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
[3495]1231      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
[2686]1232      aenv=1./(1.+(alenv*Lv/cppd))
[3495]1233      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1234      sth=0.
[3493]1235     
[2686]1236
1237      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
[3495]1238      sigma2s=0.
[2686]1239
[3493]1240      sqrt2pi=sqrt(2.*pi)
1241      xenv=senv/(sqrt(2.)*sigma1s)
[2686]1242      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
[3495]1243      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[3493]1244      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
[2686]1245     
1246      if (ctot(ind1,ind2).lt.1.e-3) then
1247      ctot(ind1,ind2)=0.
[3495]1248      qcloud(ind1)=zqsatenv(ind1,ind2)
[2686]1249
1250      else   
1251               
[3495]1252!      ctot(ind1,ind2)=ctot(ind1,ind2)
[2686]1253      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1254
[3493]1255      endif 
[2686]1256 
[3493]1257
1258
1259
[2686]1260      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
[2958]1261      ! Outputs used to check the PDFs
1262      cloudth_senv(ind1,ind2) = senv
1263      cloudth_sth(ind1,ind2) = sth
1264      cloudth_sigmaenv(ind1,ind2) = sigma1s
1265      cloudth_sigmath(ind1,ind2) = sigma2s
1266
[2686]1267      enddo       ! from the loop on ngrid l.333
1268      return
1269!     end
1270END SUBROUTINE cloudth_vert_v3
1271!
[3493]1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283       SUBROUTINE cloudth_v6(ngrid,klev,ind2,  &
1284     &           ztv,po,zqta,fraca, &
1285     &           qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1286     &           ratqs,zqs,T)
1287
1288
1289      USE ioipsl_getin_p_mod, ONLY : getin_p
1290      USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
1291     &                                cloudth_sigmath,cloudth_sigmaenv
1292
1293      IMPLICIT NONE
1294
1295#include "YOMCST.h"
1296#include "YOETHF.h"
1297#include "FCTTRE.h"
1298#include "thermcell.h"
1299#include "nuage.h"
1300
1301
1302        !Domain variables
1303      INTEGER ngrid !indice Max lat-lon
1304      INTEGER klev  !indice Max alt
1305      INTEGER ind1  !indice in [1:ngrid]
1306      INTEGER ind2  !indice in [1:klev]
1307        !thermal plume fraction
1308      REAL fraca(ngrid,klev+1)   !thermal plumes fraction in the gridbox
1309        !temperatures
1310      REAL T(ngrid,klev)       !temperature
1311      REAL zpspsk(ngrid,klev)  !factor (p/p0)**kappa (used for potential variables)
1312      REAL ztv(ngrid,klev)     !potential temperature (voir thermcell_env.F90)
1313      REAL ztla(ngrid,klev)    !liquid temperature in the thermals (Tl_th)
1314      REAL zthl(ngrid,klev)    !liquid temperature in the environment (Tl_env)
1315        !pressure
1316      REAL paprs(ngrid,klev+1)   !pressure at the interface of levels
1317      REAL pplay(ngrid,klev)     !pressure at the middle of the level
1318        !humidity
1319      REAL ratqs(ngrid,klev)   !width of the total water subgrid-scale distribution
1320      REAL po(ngrid)           !total water (qt)
1321      REAL zqenv(ngrid)        !total water in the environment (qt_env)
1322      REAL zqta(ngrid,klev)    !total water in the thermals (qt_th)
1323      REAL zqsatth(ngrid,klev)   !water saturation level in the thermals (q_sat_th)
1324      REAL zqsatenv(ngrid,klev)  !water saturation level in the environment (q_sat_env)
1325      REAL qlth(ngrid,klev)    !condensed water in the thermals
1326      REAL qlenv(ngrid,klev)   !condensed water in the environment
1327      REAL qltot(ngrid,klev)   !condensed water in the gridbox
1328        !cloud fractions
1329      REAL cth_vol(ngrid,klev)   !cloud fraction by volume in the thermals
1330      REAL cenv_vol(ngrid,klev)  !cloud fraction by volume in the environment
1331      REAL ctot_vol(ngrid,klev)  !cloud fraction by volume in the gridbox
1332      REAL cth_surf(ngrid,klev)  !cloud fraction by surface in the thermals
1333      REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment 
1334      REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
1335        !PDF of saturation deficit variables
1336      REAL rdd,cppd,Lv
1337      REAL Tbef,zdelta,qsatbef,zcor
1338      REAL alth,alenv,ath,aenv
1339      REAL sth,senv              !saturation deficits in the thermals and environment
1340      REAL sigma_env,sigma_th    !standard deviations of the biGaussian PDF
1341        !cloud fraction variables
1342      REAL xth,xenv
1343      REAL inverse_rho,beta                                  !Neggers et al. (2011) method
1344      REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
1345        !Incloud total water variables
1346      REAL zqs(ngrid)    !q_sat
1347      REAL qcloud(ngrid) !eau totale dans le nuage
1348        !Some arithmetic variables
1349      REAL erf,pi,sqrt2,sqrt2pi
1350        !Depth of the layer
1351      REAL dz(ngrid,klev)    !epaisseur de la couche en metre
1352      REAL rhodz(ngrid,klev)
1353      REAL zrho(ngrid,klev)
1354      DO ind1 = 1, ngrid
1355        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
1356        zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd          ![kg/m3]
1357        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2)            ![m]
1358      END DO
1359
1360!------------------------------------------------------------------------------
1361! Initialization
1362!------------------------------------------------------------------------------
1363      qlth(:,:)=0.
1364      qlenv(:,:)=0. 
1365      qltot(:,:)=0.
1366      cth_vol(:,:)=0.
1367      cenv_vol(:,:)=0.
1368      ctot_vol(:,:)=0.
1369      cth_surf(:,:)=0.
1370      cenv_surf(:,:)=0.
1371      ctot_surf(:,:)=0.
1372      qcloud(:)=0.
1373      rdd=287.04
1374      cppd=1005.7
1375      pi=3.14159
1376      Lv=2.5e6
1377      sqrt2=sqrt(2.)
1378      sqrt2pi=sqrt(2.*pi)
1379
1380
1381      DO ind1=1,ngrid
1382!-------------------------------------------------------------------------------
1383!Both thermal and environment in the gridbox
1384!-------------------------------------------------------------------------------
1385      IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
1386        !--------------------------------------------
1387        !calcul de qsat_env
1388        !--------------------------------------------
1389      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1390      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1391      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1392      qsatbef=MIN(0.5,qsatbef)
1393      zcor=1./(1.-retv*qsatbef)
1394      qsatbef=qsatbef*zcor
1395      zqsatenv(ind1,ind2)=qsatbef
1396        !--------------------------------------------
1397        !calcul de s_env
1398        !--------------------------------------------
1399      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
1400      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
1401      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
1402        !--------------------------------------------
1403        !calcul de qsat_th
1404        !--------------------------------------------
1405      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1406      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1407      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1408      qsatbef=MIN(0.5,qsatbef)
1409      zcor=1./(1.-retv*qsatbef)
1410      qsatbef=qsatbef*zcor
1411      zqsatth(ind1,ind2)=qsatbef
1412        !--------------------------------------------
1413        !calcul de s_th 
1414        !--------------------------------------------
1415      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84 these Arnaud Jam
1416      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84 these Arnaud Jam
1417      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84 these Arnaud Jam
1418        !--------------------------------------------
1419        !calcul standard deviations bi-Gaussian PDF
1420        !--------------------------------------------
1421      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)
1422      sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5)+ratqs(ind1,ind2)*po(ind1)
1423      xth=sth/(sqrt2*sigma_th)
1424      xenv=senv/(sqrt2*sigma_env)
1425        !--------------------------------------------
1426        !Cloud fraction by volume CF_vol
1427        !--------------------------------------------
1428      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
1429      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1430      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1431        !--------------------------------------------
1432        !Condensed water qc
1433        !--------------------------------------------
1434      qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2))
1435      qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2)) 
1436      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1437        !--------------------------------------------
1438        !Cloud fraction by surface CF_surf
1439        !--------------------------------------------
1440        !Method Neggers et al. (2011) : ok for cumulus clouds only
1441      !beta=0.0044 (Jouhaud et al.2018)
1442      !inverse_rho=1.+beta*dz(ind1,ind2)
1443      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1444        !Method Brooks et al. (2005) : ok for all types of clouds
1445      a_Brooks=0.6694
1446      b_Brooks=0.1882
1447      A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
1448      Dx_Brooks=200000.   !-- si l'on considere des mailles de 200km de cote
1449      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
1450      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1451        !--------------------------------------------
1452        !Incloud Condensed water qcloud
1453        !--------------------------------------------
1454      if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
1455      ctot_vol(ind1,ind2)=0.
1456      ctot_surf(ind1,ind2)=0.
1457      qcloud(ind1)=zqsatenv(ind1,ind2)
1458      else
1459      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
1460      endif
1461
1462
1463
1464!-------------------------------------------------------------------------------
1465!Environment only in the gridbox
1466!-------------------------------------------------------------------------------
1467      ELSE
1468        !--------------------------------------------
1469        !calcul de qsat_env
1470        !--------------------------------------------
1471      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1472      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1473      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1474      qsatbef=MIN(0.5,qsatbef)
1475      zcor=1./(1.-retv*qsatbef)
1476      qsatbef=qsatbef*zcor
1477      zqsatenv(ind1,ind2)=qsatbef
1478        !--------------------------------------------
1479        !calcul de s_env
1480        !--------------------------------------------
1481      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
1482      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
1483      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
1484        !--------------------------------------------
1485        !calcul standard deviations Gaussian PDF
1486        !--------------------------------------------
1487      zqenv(ind1)=po(ind1)
1488      sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
1489      xenv=senv/(sqrt2*sigma_env)
1490        !--------------------------------------------
1491        !Cloud fraction by volume CF_vol
1492        !--------------------------------------------
1493      ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1494        !--------------------------------------------
1495        !Condensed water qc
1496        !--------------------------------------------
1497      qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2))
1498        !--------------------------------------------
1499        !Cloud fraction by surface CF_surf
1500        !--------------------------------------------
1501        !Method Neggers et al. (2011) : ok for cumulus clouds only
1502      !beta=0.0044 (Jouhaud et al.2018)
1503      !inverse_rho=1.+beta*dz(ind1,ind2)
1504      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1505        !Method Brooks et al. (2005) : ok for all types of clouds
1506      a_Brooks=0.6694
1507      b_Brooks=0.1882
1508      A_Maj_Brooks=0.1635 !-- sans dependence au shear
1509      Dx_Brooks=200000.
1510      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
1511      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1512        !--------------------------------------------
1513        !Incloud Condensed water qcloud
1514        !--------------------------------------------
1515      if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
1516      ctot_vol(ind1,ind2)=0.
1517      ctot_surf(ind1,ind2)=0.
1518      qcloud(ind1)=zqsatenv(ind1,ind2)
1519      else
1520      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
1521      endif
1522
1523
1524      END IF  ! From the separation (thermal/envrionnement) et (environnement only)
1525
1526      ! Outputs used to check the PDFs
1527      cloudth_senv(ind1,ind2) = senv
1528      cloudth_sth(ind1,ind2) = sth
1529      cloudth_sigmaenv(ind1,ind2) = sigma_env
1530      cloudth_sigmath(ind1,ind2) = sigma_th
1531
1532      END DO  ! From the loop on ngrid
1533      return
1534
1535END SUBROUTINE cloudth_v6
[2686]1536END MODULE cloudth_mod
[3493]1537
1538
1539
1540
Note: See TracBrowser for help on using the repository browser.