source: LMDZ5/branches/IPSLCM6.0.11/libf/phylmd/cloudth_mod.F90 @ 5416

Last change on this file since 5416 was 2911, checked in by jbmadeleine, 8 years ago

Changes in the vertical PDF used for shallow convective clouds.
Added vert_alpha_th which is the variance of the vertical PDF used for thermals.
If nothing is specified in the def files, it is set equal to vert_alpha.
Otherwise, it is now possible to have a different vertical PDF in the thermals
and in the environment.
Also added a fourth option iflag_cloudth_vert = 4 where
deltas=vert_alpha*sigmas instead of deltas=a*vert_alpha*sigmas.
File size: 38.9 KB
Line 
1MODULE cloudth_mod
2
3  IMPLICIT NONE
4
5CONTAINS
6
7       SUBROUTINE cloudth(ngrid,klev,ind2,  &
8     &           ztv,po,zqta,fraca, &
9     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
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)
40      REAL ztla(ngrid,klev)
41      REAL zthl(ngrid,klev)
42
43      REAL zqsatth(ngrid,klev)
44      REAL zqsatenv(ngrid,klev)
45     
46     
47      REAL sigma1(ngrid,klev)
48      REAL sigma2(ngrid,klev)
49      REAL qlth(ngrid,klev)
50      REAL qlenv(ngrid,klev)
51      REAL qltot(ngrid,klev)
52      REAL cth(ngrid,klev) 
53      REAL cenv(ngrid,klev)   
54      REAL ctot(ngrid,klev)
55      REAL rneb(ngrid,klev)
56      REAL t(ngrid,klev)
57      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
58      REAL rdd,cppd,Lv
59      REAL alth,alenv,ath,aenv
60      REAL sth,senv,sigma1s,sigma2s,xth,xenv
61      REAL Tbef,zdelta,qsatbef,zcor
62      REAL qlbef 
63      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
64     
65      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
66      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
67      REAL zqs(ngrid), qcloud(ngrid)
68      REAL erf
69
70
71
72
73!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
74! Gestion de deux versions de cloudth
75!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
76
77      IF (iflag_cloudth_vert.GE.1) THEN
78      CALL cloudth_vert(ngrid,klev,ind2,  &
79     &           ztv,po,zqta,fraca, &
80     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
81     &           ratqs,zqs,t)
82      RETURN
83      ENDIF
84!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
85
86
87!-------------------------------------------------------------------------------
88! Initialisation des variables r?elles
89!-------------------------------------------------------------------------------
90      sigma1(:,:)=0.
91      sigma2(:,:)=0.
92      qlth(:,:)=0.
93      qlenv(:,:)=0. 
94      qltot(:,:)=0.
95      rneb(:,:)=0.
96      qcloud(:)=0.
97      cth(:,:)=0.
98      cenv(:,:)=0.
99      ctot(:,:)=0.
100      qsatmmussig1=0.
101      qsatmmussig2=0.
102      rdd=287.04
103      cppd=1005.7
104      pi=3.14159
105      Lv=2.5e6
106      sqrt2pi=sqrt(2.*pi)
107
108
109
110!-------------------------------------------------------------------------------
111! Calcul de la fraction du thermique et des ?cart-types des distributions
112!-------------------------------------------------------------------------------                 
113      do ind1=1,ngrid
114
115      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
116
117      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
118
119
120!      zqenv(ind1)=po(ind1)
121      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
122      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
123      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
124      qsatbef=MIN(0.5,qsatbef)
125      zcor=1./(1.-retv*qsatbef)
126      qsatbef=qsatbef*zcor
127      zqsatenv(ind1,ind2)=qsatbef
128
129
130
131
132      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
133      aenv=1./(1.+(alenv*Lv/cppd))
134      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
135
136
137
138
139      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
140      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
141      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
142      qsatbef=MIN(0.5,qsatbef)
143      zcor=1./(1.-retv*qsatbef)
144      qsatbef=qsatbef*zcor
145      zqsatth(ind1,ind2)=qsatbef
146           
147      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
148      ath=1./(1.+(alth*Lv/cppd))
149      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
150     
151     
152
153!------------------------------------------------------------------------------
154! Calcul des ?cart-types pour s
155!------------------------------------------------------------------------------
156
157!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
158!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
159!       if (paprs(ind1,ind2).gt.90000) then
160!       ratqs(ind1,ind2)=0.002
161!       else
162!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
163!       endif
164       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
165       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
166!       sigma1s=ratqs(ind1,ind2)*po(ind1)
167!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
168 
169!------------------------------------------------------------------------------
170! Calcul de l'eau condens?e et de la couverture nuageuse
171!------------------------------------------------------------------------------
172      sqrt2pi=sqrt(2.*pi)
173      xth=sth/(sqrt(2.)*sigma2s)
174      xenv=senv/(sqrt(2.)*sigma1s)
175      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
176      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
177      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
178
179      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
180      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
181      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
182
183!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
184      if (ctot(ind1,ind2).lt.1.e-10) then
185      ctot(ind1,ind2)=0.
186      qcloud(ind1)=zqsatenv(ind1,ind2)
187
188      else   
189               
190      ctot(ind1,ind2)=ctot(ind1,ind2)
191      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
192
193      endif                           
194     
195         
196!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
197
198
199      else  ! gaussienne environnement seule
200     
201      zqenv(ind1)=po(ind1)
202      Tbef=t(ind1,ind2)
203      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
204      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
205      qsatbef=MIN(0.5,qsatbef)
206      zcor=1./(1.-retv*qsatbef)
207      qsatbef=qsatbef*zcor
208      zqsatenv(ind1,ind2)=qsatbef
209     
210
211!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
212      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
213      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
214      aenv=1./(1.+(alenv*Lv/cppd))
215      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
216     
217
218      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
219
220      sqrt2pi=sqrt(2.*pi)
221      xenv=senv/(sqrt(2.)*sigma1s)
222      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
223      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
224     
225      if (ctot(ind1,ind2).lt.1.e-3) then
226      ctot(ind1,ind2)=0.
227      qcloud(ind1)=zqsatenv(ind1,ind2)
228
229      else   
230               
231      ctot(ind1,ind2)=ctot(ind1,ind2)
232      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
233
234      endif   
235 
236 
237 
238 
239 
240 
241      endif   
242      enddo
243     
244      return
245!     end
246END SUBROUTINE cloudth
247
248
249
250!===========================================================================
251     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
252     &           ztv,po,zqta,fraca, &
253     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
254     &           ratqs,zqs,t)
255
256!===========================================================================
257! Auteur : Arnaud Octavio Jam (LMD/CNRS)
258! Date : 25 Mai 2010
259! Objet : calcule les valeurs de qc et rneb dans les thermiques
260!===========================================================================
261
262
263      USE ioipsl_getin_p_mod, ONLY : getin_p
264
265      IMPLICIT NONE
266
267#include "YOMCST.h"
268#include "YOETHF.h"
269#include "FCTTRE.h"
270#include "thermcell.h"
271#include "nuage.h"
272     
273      INTEGER itap,ind1,ind2
274      INTEGER ngrid,klev,klon,l,ig
275     
276      REAL ztv(ngrid,klev)
277      REAL po(ngrid)
278      REAL zqenv(ngrid)   
279      REAL zqta(ngrid,klev)
280         
281      REAL fraca(ngrid,klev+1)
282      REAL zpspsk(ngrid,klev)
283      REAL paprs(ngrid,klev+1)
284      REAL ztla(ngrid,klev)
285      REAL zthl(ngrid,klev)
286
287      REAL zqsatth(ngrid,klev)
288      REAL zqsatenv(ngrid,klev)
289     
290     
291      REAL sigma1(ngrid,klev)                                                         
292      REAL sigma2(ngrid,klev)
293      REAL qlth(ngrid,klev)
294      REAL qlenv(ngrid,klev)
295      REAL qltot(ngrid,klev)
296      REAL cth(ngrid,klev) 
297      REAL cenv(ngrid,klev)   
298      REAL ctot(ngrid,klev)
299      REAL rneb(ngrid,klev)
300      REAL t(ngrid,klev)                                                                 
301      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
302      REAL rdd,cppd,Lv,sqrt2,sqrtpi
303      REAL alth,alenv,ath,aenv
304      REAL sth,senv,sigma1s,sigma2s,xth,xenv
305      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
306      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
307      REAL Tbef,zdelta,qsatbef,zcor
308      REAL qlbef 
309      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
310      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
311      ! (J Jouhaud, JL Dufresne, JB Madeleine)
312      REAL,SAVE :: vert_alpha
313      !$OMP THREADPRIVATE(vert_alpha)
314      LOGICAL, SAVE :: firstcall = .TRUE.
315      !$OMP THREADPRIVATE(firstcall)
316     
317      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
318      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
319      REAL zqs(ngrid), qcloud(ngrid)
320      REAL erf
321
322!------------------------------------------------------------------------------
323! Initialisation des variables r?elles
324!------------------------------------------------------------------------------
325      sigma1(:,:)=0.
326      sigma2(:,:)=0.
327      qlth(:,:)=0.
328      qlenv(:,:)=0. 
329      qltot(:,:)=0.
330      rneb(:,:)=0.
331      qcloud(:)=0.
332      cth(:,:)=0.
333      cenv(:,:)=0.
334      ctot(:,:)=0.
335      qsatmmussig1=0.
336      qsatmmussig2=0.
337      rdd=287.04
338      cppd=1005.7
339      pi=3.14159
340      Lv=2.5e6
341      sqrt2pi=sqrt(2.*pi)
342      sqrt2=sqrt(2.)
343      sqrtpi=sqrt(pi)
344
345      IF (firstcall) THEN
346        vert_alpha=0.5
347        CALL getin_p('cloudth_vert_alpha',vert_alpha)
348        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
349        firstcall=.FALSE.
350      ENDIF
351
352!-------------------------------------------------------------------------------
353! Calcul de la fraction du thermique et des ?cart-types des distributions
354!-------------------------------------------------------------------------------                 
355      do ind1=1,ngrid
356
357      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
358
359      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
360
361
362!      zqenv(ind1)=po(ind1)
363      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
364      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
365      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
366      qsatbef=MIN(0.5,qsatbef)
367      zcor=1./(1.-retv*qsatbef)
368      qsatbef=qsatbef*zcor
369      zqsatenv(ind1,ind2)=qsatbef
370
371
372
373
374      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
375      aenv=1./(1.+(alenv*Lv/cppd))
376      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
377
378
379
380
381      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
382      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
383      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
384      qsatbef=MIN(0.5,qsatbef)
385      zcor=1./(1.-retv*qsatbef)
386      qsatbef=qsatbef*zcor
387      zqsatth(ind1,ind2)=qsatbef
388           
389      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
390      ath=1./(1.+(alth*Lv/cppd))
391      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
392     
393     
394
395!------------------------------------------------------------------------------
396! Calcul des ?cart-types pour s
397!------------------------------------------------------------------------------
398
399      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
400      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
401!       if (paprs(ind1,ind2).gt.90000) then
402!       ratqs(ind1,ind2)=0.002
403!       else
404!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
405!       endif
406!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
407!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
408!       sigma1s=ratqs(ind1,ind2)*po(ind1)
409!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
410 
411!------------------------------------------------------------------------------
412! Calcul de l'eau condens?e et de la couverture nuageuse
413!------------------------------------------------------------------------------
414      sqrt2pi=sqrt(2.*pi)
415      xth=sth/(sqrt(2.)*sigma2s)
416      xenv=senv/(sqrt(2.)*sigma1s)
417      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
418      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
419      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
420
421      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
422      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
423      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
424     
425       IF (iflag_cloudth_vert == 1) THEN
426!-------------------------------------------------------------------------------
427!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
428!-------------------------------------------------------------------------------
429!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
430!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
431      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
432      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
433!      deltasenv=aenv*0.01*po(ind1)
434!     deltasth=ath*0.01*zqta(ind1,ind2)   
435      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
436      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
437      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
438      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
439      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
440      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
441     
442      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
443      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
444      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
445
446      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
447      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
448      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
449      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
450
451      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
452!      qlenv(ind1,ind2)=IntJ
453!      print*, qlenv(ind1,ind2),'VERIF EAU'
454
455
456      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
457!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
458!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
459      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
460      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
461      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
462      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
463!      qlth(ind1,ind2)=IntJ
464!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
465      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
466
467      ELSE IF (iflag_cloudth_vert == 2) THEN
468
469!-------------------------------------------------------------------------------
470!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
471!-------------------------------------------------------------------------------
472!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
473!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
474!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
475!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
476      deltasenv=aenv*vert_alpha*sigma1s
477      deltasth=ath*vert_alpha*sigma2s
478     
479      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
480      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
481      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
482      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
483!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
484!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
485     
486      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
487      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
488      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
489
490      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
491      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
492      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
493      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
494
495!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
496!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
497!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
498
499      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
500!      qlenv(ind1,ind2)=IntJ
501!      print*, qlenv(ind1,ind2),'VERIF EAU'
502
503      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
504      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
505      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
506      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
507     
508      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
509!      qlth(ind1,ind2)=IntJ
510!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
511      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
512     
513
514
515
516      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
517
518!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
519
520      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
521      ctot(ind1,ind2)=0.
522      qcloud(ind1)=zqsatenv(ind1,ind2)
523
524      else
525               
526      ctot(ind1,ind2)=ctot(ind1,ind2)
527      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
528!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
529!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
530
531      endif 
532                       
533     
534         
535!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
536
537
538      else  ! gaussienne environnement seule
539     
540      zqenv(ind1)=po(ind1)
541      Tbef=t(ind1,ind2)
542      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
543      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
544      qsatbef=MIN(0.5,qsatbef)
545      zcor=1./(1.-retv*qsatbef)
546      qsatbef=qsatbef*zcor
547      zqsatenv(ind1,ind2)=qsatbef
548     
549
550!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
551      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
552      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
553      aenv=1./(1.+(alenv*Lv/cppd))
554      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
555     
556
557      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
558
559      sqrt2pi=sqrt(2.*pi)
560      xenv=senv/(sqrt(2.)*sigma1s)
561      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
562      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
563     
564      if (ctot(ind1,ind2).lt.1.e-3) then
565      ctot(ind1,ind2)=0.
566      qcloud(ind1)=zqsatenv(ind1,ind2)
567
568      else   
569               
570      ctot(ind1,ind2)=ctot(ind1,ind2)
571      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
572
573      endif   
574 
575 
576 
577 
578 
579 
580      endif   
581      enddo
582     
583      return
584!     end
585END SUBROUTINE cloudth_vert
586
587       SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
588     &           ztv,po,zqta,fraca, &
589     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
590     &           ratqs,zqs,t)
591
592
593      IMPLICIT NONE
594
595
596!===========================================================================
597! Author : Arnaud Octavio Jam (LMD/CNRS)
598! Date : 25 Mai 2010
599! Objet : calcule les valeurs de qc et rneb dans les thermiques
600!===========================================================================
601
602
603#include "YOMCST.h"
604#include "YOETHF.h"
605#include "FCTTRE.h"
606#include "thermcell.h"
607#include "nuage.h"
608
609      INTEGER itap,ind1,ind2
610      INTEGER ngrid,klev,klon,l,ig
611     
612      REAL ztv(ngrid,klev)
613      REAL po(ngrid)
614      REAL zqenv(ngrid)   
615      REAL zqta(ngrid,klev)
616         
617      REAL fraca(ngrid,klev+1)
618      REAL zpspsk(ngrid,klev)
619      REAL paprs(ngrid,klev+1)
620      REAL ztla(ngrid,klev)
621      REAL zthl(ngrid,klev)
622
623      REAL zqsatth(ngrid,klev)
624      REAL zqsatenv(ngrid,klev)
625     
626     
627      REAL sigma1(ngrid,klev)
628      REAL sigma2(ngrid,klev)
629      REAL qlth(ngrid,klev)
630      REAL qlenv(ngrid,klev)
631      REAL qltot(ngrid,klev)
632      REAL cth(ngrid,klev) 
633      REAL cenv(ngrid,klev)   
634      REAL ctot(ngrid,klev)
635      REAL rneb(ngrid,klev)
636      REAL t(ngrid,klev)
637      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
638      REAL rdd,cppd,Lv
639      REAL alth,alenv,ath,aenv
640      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
641      REAL Tbef,zdelta,qsatbef,zcor
642      REAL qlbef 
643      REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
644      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
645      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
646      REAL zqs(ngrid), qcloud(ngrid)
647      REAL erf
648
649
650
651      IF (iflag_cloudth_vert.GE.1) THEN
652      CALL cloudth_vert_v3(ngrid,klev,ind2,  &
653     &           ztv,po,zqta,fraca, &
654     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
655     &           ratqs,zqs,t)
656      RETURN
657      ENDIF
658!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
659
660
661!-------------------------------------------------------------------------------
662! Initialisation des variables r?elles
663!-------------------------------------------------------------------------------
664      sigma1(:,:)=0.
665      sigma2(:,:)=0.
666      qlth(:,:)=0.
667      qlenv(:,:)=0. 
668      qltot(:,:)=0.
669      rneb(:,:)=0.
670      qcloud(:)=0.
671      cth(:,:)=0.
672      cenv(:,:)=0.
673      ctot(:,:)=0.
674      qsatmmussig1=0.
675      qsatmmussig2=0.
676      rdd=287.04
677      cppd=1005.7
678      pi=3.14159
679      Lv=2.5e6
680      sqrt2pi=sqrt(2.*pi)
681      sqrt2=sqrt(2.)
682      sqrtpi=sqrt(pi)
683
684
685!-------------------------------------------------------------------------------
686! Cloud fraction in the thermals and standard deviation of the PDFs
687!-------------------------------------------------------------------------------                 
688      do ind1=1,ngrid
689
690      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
691
692      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
693
694      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
695      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
696      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
697      qsatbef=MIN(0.5,qsatbef)
698      zcor=1./(1.-retv*qsatbef)
699      qsatbef=qsatbef*zcor
700      zqsatenv(ind1,ind2)=qsatbef
701
702
703      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
704      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
705      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
706
707!po = qt de l'environnement ET des thermique
708!zqenv = qt environnement
709!zqsatenv = qsat environnement
710!zthl = Tl environnement
711
712
713      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
714      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
715      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
716      qsatbef=MIN(0.5,qsatbef)
717      zcor=1./(1.-retv*qsatbef)
718      qsatbef=qsatbef*zcor
719      zqsatth(ind1,ind2)=qsatbef
720           
721      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
722      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
723      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
724     
725!zqta = qt thermals
726!zqsatth = qsat thermals
727!ztla = Tl thermals
728
729!------------------------------------------------------------------------------
730! s standard deviations
731!------------------------------------------------------------------------------
732
733!     tests
734!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
735!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
736!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
737!     final option
738      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
739      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
740 
741!------------------------------------------------------------------------------
742! Condensed water and cloud cover
743!------------------------------------------------------------------------------
744      xth=sth/(sqrt2*sigma2s)
745      xenv=senv/(sqrt2*sigma1s)
746      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
747      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
748      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
749
750      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
751      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
752      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
753
754      if (ctot(ind1,ind2).lt.1.e-10) then
755      ctot(ind1,ind2)=0.
756      qcloud(ind1)=zqsatenv(ind1,ind2)
757      else
758      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
759      endif
760
761      else  ! Environnement only, follow the if l.110
762     
763      zqenv(ind1)=po(ind1)
764      Tbef=t(ind1,ind2)
765      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
766      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
767      qsatbef=MIN(0.5,qsatbef)
768      zcor=1./(1.-retv*qsatbef)
769      qsatbef=qsatbef*zcor
770      zqsatenv(ind1,ind2)=qsatbef
771
772!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
773      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
774      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
775      aenv=1./(1.+(alenv*Lv/cppd))
776      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
777     
778      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
779
780      xenv=senv/(sqrt2*sigma1s)
781      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
782      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
783
784      if (ctot(ind1,ind2).lt.1.e-3) then
785      ctot(ind1,ind2)=0.
786      qcloud(ind1)=zqsatenv(ind1,ind2)
787      else   
788      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
789      endif
790
791
792      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
793      enddo       ! from the loop on ngrid l.108
794      return
795!     end
796END SUBROUTINE cloudth_v3
797
798
799
800!===========================================================================
801     SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
802     &           ztv,po,zqta,fraca, &
803     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
804     &           ratqs,zqs,t)
805
806!===========================================================================
807! Auteur : Arnaud Octavio Jam (LMD/CNRS)
808! Date : 25 Mai 2010
809! Objet : calcule les valeurs de qc et rneb dans les thermiques
810!===========================================================================
811
812
813      USE ioipsl_getin_p_mod, ONLY : getin_p
814
815      IMPLICIT NONE
816
817#include "YOMCST.h"
818#include "YOETHF.h"
819#include "FCTTRE.h"
820#include "thermcell.h"
821#include "nuage.h"
822     
823      INTEGER itap,ind1,ind2
824      INTEGER ngrid,klev,klon,l,ig
825     
826      REAL ztv(ngrid,klev)
827      REAL po(ngrid)
828      REAL zqenv(ngrid)   
829      REAL zqta(ngrid,klev)
830         
831      REAL fraca(ngrid,klev+1)
832      REAL zpspsk(ngrid,klev)
833      REAL paprs(ngrid,klev+1)
834      REAL ztla(ngrid,klev)
835      REAL zthl(ngrid,klev)
836
837      REAL zqsatth(ngrid,klev)
838      REAL zqsatenv(ngrid,klev)
839     
840     
841      REAL sigma1(ngrid,klev)                                                         
842      REAL sigma2(ngrid,klev)
843      REAL qlth(ngrid,klev)
844      REAL qlenv(ngrid,klev)
845      REAL qltot(ngrid,klev)
846      REAL cth(ngrid,klev) 
847      REAL cenv(ngrid,klev)   
848      REAL ctot(ngrid,klev)
849      REAL rneb(ngrid,klev)
850      REAL t(ngrid,klev)                                                                 
851      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
852      REAL rdd,cppd,Lv
853      REAL alth,alenv,ath,aenv
854      REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
855      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
856      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
857      REAL Tbef,zdelta,qsatbef,zcor
858      REAL qlbef 
859      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
860      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
861      ! (J Jouhaud, JL Dufresne, JB Madeleine)
862      REAL,SAVE :: vert_alpha, vert_alpha_th
863      !$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
864      LOGICAL, SAVE :: firstcall = .TRUE.
865      !$OMP THREADPRIVATE(firstcall)
866
867      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
868      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
869      REAL zqs(ngrid), qcloud(ngrid)
870      REAL erf
871
872!------------------------------------------------------------------------------
873! Initialize
874!------------------------------------------------------------------------------
875      sigma1(:,:)=0.
876      sigma2(:,:)=0.
877      qlth(:,:)=0.
878      qlenv(:,:)=0. 
879      qltot(:,:)=0.
880      rneb(:,:)=0.
881      qcloud(:)=0.
882      cth(:,:)=0.
883      cenv(:,:)=0.
884      ctot(:,:)=0.
885      qsatmmussig1=0.
886      qsatmmussig2=0.
887      rdd=287.04
888      cppd=1005.7
889      pi=3.14159
890      Lv=2.5e6
891      sqrt2pi=sqrt(2.*pi)
892      sqrt2=sqrt(2.)
893      sqrtpi=sqrt(pi)
894
895      IF (firstcall) THEN
896        vert_alpha=0.5
897        CALL getin_p('cloudth_vert_alpha',vert_alpha)
898        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
899        ! The factor used for the thermal is equal to that of the environment
900        !   if nothing is explicitly specified in the def file
901        vert_alpha_th=vert_alpha
902        CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
903        WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th
904        firstcall=.FALSE.
905      ENDIF
906
907!-------------------------------------------------------------------------------
908! Calcul de la fraction du thermique et des ecart-types des distributions
909!-------------------------------------------------------------------------------                 
910      do ind1=1,ngrid
911
912      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
913
914      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
915
916
917      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
918      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
919      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
920      qsatbef=MIN(0.5,qsatbef)
921      zcor=1./(1.-retv*qsatbef)
922      qsatbef=qsatbef*zcor
923      zqsatenv(ind1,ind2)=qsatbef
924
925
926      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
927      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
928      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
929
930!zqenv = qt environnement
931!zqsatenv = qsat environnement
932!zthl = Tl environnement
933
934
935      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
936      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
937      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
938      qsatbef=MIN(0.5,qsatbef)
939      zcor=1./(1.-retv*qsatbef)
940      qsatbef=qsatbef*zcor
941      zqsatth(ind1,ind2)=qsatbef
942           
943      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
944      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
945      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
946     
947     
948!zqta = qt thermals
949!zqsatth = qsat thermals
950!ztla = Tl thermals
951
952!------------------------------------------------------------------------------
953! s standard deviation
954!------------------------------------------------------------------------------
955
956      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
957      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
958!      tests
959!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
960!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
961!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
962!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
963!       if (paprs(ind1,ind2).gt.90000) then
964!       ratqs(ind1,ind2)=0.002
965!       else
966!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
967!       endif
968!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
969!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
970!       sigma1s=ratqs(ind1,ind2)*po(ind1)
971!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
972 
973       IF (iflag_cloudth_vert == 1) THEN
974!-------------------------------------------------------------------------------
975!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
976!-------------------------------------------------------------------------------
977
978      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
979      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
980
981      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
982      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
983      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
984      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
985      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
986      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
987     
988      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
989      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
990      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
991
992      ! Environment
993      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
994      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
995      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
996      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
997
998      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
999
1000      ! Thermal
1001      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
1002      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
1003      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
1004      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
1005      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1006      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1007
1008      ELSE IF (iflag_cloudth_vert >= 3) THEN
1009
1010!-------------------------------------------------------------------------------
1011!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1012!-------------------------------------------------------------------------------
1013!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
1014!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
1015!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1016!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1017      IF (iflag_cloudth_vert == 3) THEN
1018        deltasenv=aenv*vert_alpha*sigma1s
1019        deltasth=ath*vert_alpha_th*sigma2s
1020      ELSE IF (iflag_cloudth_vert == 4) THEN
1021        deltasenv=vert_alpha*sigma1s
1022        deltasth=vert_alpha_th*sigma2s
1023      ENDIF
1024     
1025      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1026      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1027      exp_xenv1 = exp(-1.*xenv1**2)
1028      exp_xenv2 = exp(-1.*xenv2**2)
1029      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1030      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1031      exp_xth1 = exp(-1.*xth1**2)
1032      exp_xth2 = exp(-1.*xth2**2)
1033     
1034      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
1035      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
1036      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1037
1038     
1039      !environnement
1040      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1041      if (deltasenv .lt. 1.e-10) then
1042      qlenv(ind1,ind2)=IntJ
1043      else
1044      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
1045      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
1046      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
1047      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1048      endif
1049
1050
1051      !thermique
1052      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1053      if (deltasth .lt. 1.e-10) then ! Seuil a choisir !!!
1054      qlth(ind1,ind2)=IntJ
1055      else
1056      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1057      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1058      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
1059      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1060      endif
1061
1062      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1063
1064
1065      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1066
1067      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1068      ctot(ind1,ind2)=0.
1069      qcloud(ind1)=zqsatenv(ind1,ind2)
1070
1071      else
1072               
1073      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1074!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1075!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1076
1077      endif 
1078
1079      else  ! Environment only
1080     
1081      zqenv(ind1)=po(ind1)
1082      Tbef=t(ind1,ind2)
1083      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1084      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1085      qsatbef=MIN(0.5,qsatbef)
1086      zcor=1./(1.-retv*qsatbef)
1087      qsatbef=qsatbef*zcor
1088      zqsatenv(ind1,ind2)=qsatbef
1089     
1090
1091!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1092      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
1093      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
1094      aenv=1./(1.+(alenv*Lv/cppd))
1095      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1096     
1097
1098      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1099
1100      xenv=senv/(sqrt2*sigma1s)
1101      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1102      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
1103     
1104      if (ctot(ind1,ind2).lt.1.e-3) then
1105      ctot(ind1,ind2)=0.
1106      qcloud(ind1)=zqsatenv(ind1,ind2)
1107
1108      else   
1109               
1110      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1111
1112      endif   
1113 
1114      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1115      enddo       ! from the loop on ngrid l.333
1116     
1117      return
1118!     end
1119END SUBROUTINE cloudth_vert_v3
1120!
1121END MODULE cloudth_mod
Note: See TracBrowser for help on using the repository browser.