source: LMDZ5/branches/testing/libf/phylmd/cloudth_mod.F90 @ 3998

Last change on this file since 3998 was 2720, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2664:2719 into testing branch
File size: 38.3 KB
RevLine 
[2686]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      LOGICAL, SAVE :: firstcall = .TRUE.
314     
315      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
316      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
317      REAL zqs(ngrid), qcloud(ngrid)
318      REAL erf
319
320!------------------------------------------------------------------------------
321! Initialisation des variables r?elles
322!------------------------------------------------------------------------------
323      sigma1(:,:)=0.
324      sigma2(:,:)=0.
325      qlth(:,:)=0.
326      qlenv(:,:)=0. 
327      qltot(:,:)=0.
328      rneb(:,:)=0.
329      qcloud(:)=0.
330      cth(:,:)=0.
331      cenv(:,:)=0.
332      ctot(:,:)=0.
333      qsatmmussig1=0.
334      qsatmmussig2=0.
335      rdd=287.04
336      cppd=1005.7
337      pi=3.14159
338      Lv=2.5e6
339      sqrt2pi=sqrt(2.*pi)
340      sqrt2=sqrt(2.)
341      sqrtpi=sqrt(pi)
342
343      IF (firstcall) THEN
344        vert_alpha=0.5
345        CALL getin_p('cloudth_vert_alpha',vert_alpha)
346        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
347        firstcall=.FALSE.
348      ENDIF
349
350!-------------------------------------------------------------------------------
351! Calcul de la fraction du thermique et des ?cart-types des distributions
352!-------------------------------------------------------------------------------                 
353      do ind1=1,ngrid
354
355      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
356
357      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
358
359
360!      zqenv(ind1)=po(ind1)
361      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
362      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
363      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
364      qsatbef=MIN(0.5,qsatbef)
365      zcor=1./(1.-retv*qsatbef)
366      qsatbef=qsatbef*zcor
367      zqsatenv(ind1,ind2)=qsatbef
368
369
370
371
372      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
373      aenv=1./(1.+(alenv*Lv/cppd))
374      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
375
376
377
378
379      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
380      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
381      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
382      qsatbef=MIN(0.5,qsatbef)
383      zcor=1./(1.-retv*qsatbef)
384      qsatbef=qsatbef*zcor
385      zqsatth(ind1,ind2)=qsatbef
386           
387      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
388      ath=1./(1.+(alth*Lv/cppd))
389      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
390     
391     
392
393!------------------------------------------------------------------------------
394! Calcul des ?cart-types pour s
395!------------------------------------------------------------------------------
396
397      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
398      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
399!       if (paprs(ind1,ind2).gt.90000) then
400!       ratqs(ind1,ind2)=0.002
401!       else
402!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
403!       endif
404!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
405!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
406!       sigma1s=ratqs(ind1,ind2)*po(ind1)
407!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
408 
409!------------------------------------------------------------------------------
410! Calcul de l'eau condens?e et de la couverture nuageuse
411!------------------------------------------------------------------------------
412      sqrt2pi=sqrt(2.*pi)
413      xth=sth/(sqrt(2.)*sigma2s)
414      xenv=senv/(sqrt(2.)*sigma1s)
415      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
416      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
417      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
418
419      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
420      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
421      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
422     
423       IF (iflag_cloudth_vert == 1) THEN
424!-------------------------------------------------------------------------------
425!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
426!-------------------------------------------------------------------------------
427!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
428!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
429      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
430      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
431!      deltasenv=aenv*0.01*po(ind1)
432!     deltasth=ath*0.01*zqta(ind1,ind2)   
433      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
434      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
435      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
436      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
437      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
438      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
439     
440      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
441      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
442      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
443
444      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
445      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
446      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
447      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
448
449      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
450!      qlenv(ind1,ind2)=IntJ
451!      print*, qlenv(ind1,ind2),'VERIF EAU'
452
453
454      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
455!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
456!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
457      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
458      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
459      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
460      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
461!      qlth(ind1,ind2)=IntJ
462!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
463      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
464
465      ELSE IF (iflag_cloudth_vert == 2) THEN
466
467!-------------------------------------------------------------------------------
468!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
469!-------------------------------------------------------------------------------
470!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
471!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
472!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
473!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
474      deltasenv=aenv*vert_alpha*sigma1s
475      deltasth=ath*vert_alpha*sigma2s
476     
477      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
478      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
479      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
480      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
481!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
482!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
483     
484      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
485      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
486      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
487
488      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
489      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
490      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
491      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
492
493!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
494!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
495!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
496
497      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
498!      qlenv(ind1,ind2)=IntJ
499!      print*, qlenv(ind1,ind2),'VERIF EAU'
500
501      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
502      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
503      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
504      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
505     
506      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
507!      qlth(ind1,ind2)=IntJ
508!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
509      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
510     
511
512
513
514      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
515
516!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
517
518      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
519      ctot(ind1,ind2)=0.
520      qcloud(ind1)=zqsatenv(ind1,ind2)
521
522      else
523               
524      ctot(ind1,ind2)=ctot(ind1,ind2)
525      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
526!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
527!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
528
529      endif 
530                       
531     
532         
533!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
534
535
536      else  ! gaussienne environnement seule
537     
538      zqenv(ind1)=po(ind1)
539      Tbef=t(ind1,ind2)
540      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
541      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
542      qsatbef=MIN(0.5,qsatbef)
543      zcor=1./(1.-retv*qsatbef)
544      qsatbef=qsatbef*zcor
545      zqsatenv(ind1,ind2)=qsatbef
546     
547
548!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
549      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
550      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
551      aenv=1./(1.+(alenv*Lv/cppd))
552      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
553     
554
555      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
556
557      sqrt2pi=sqrt(2.*pi)
558      xenv=senv/(sqrt(2.)*sigma1s)
559      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
560      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
561     
562      if (ctot(ind1,ind2).lt.1.e-3) then
563      ctot(ind1,ind2)=0.
564      qcloud(ind1)=zqsatenv(ind1,ind2)
565
566      else   
567               
568      ctot(ind1,ind2)=ctot(ind1,ind2)
569      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
570
571      endif   
572 
573 
574 
575 
576 
577 
578      endif   
579      enddo
580     
581      return
582!     end
583END SUBROUTINE cloudth_vert
584
585       SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
586     &           ztv,po,zqta,fraca, &
587     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
588     &           ratqs,zqs,t)
589
590
591      IMPLICIT NONE
592
593
594!===========================================================================
595! Author : Arnaud Octavio Jam (LMD/CNRS)
596! Date : 25 Mai 2010
597! Objet : calcule les valeurs de qc et rneb dans les thermiques
598!===========================================================================
599
600
601#include "YOMCST.h"
602#include "YOETHF.h"
603#include "FCTTRE.h"
604#include "thermcell.h"
605#include "nuage.h"
606
607      INTEGER itap,ind1,ind2
608      INTEGER ngrid,klev,klon,l,ig
609     
610      REAL ztv(ngrid,klev)
611      REAL po(ngrid)
612      REAL zqenv(ngrid)   
613      REAL zqta(ngrid,klev)
614         
615      REAL fraca(ngrid,klev+1)
616      REAL zpspsk(ngrid,klev)
617      REAL paprs(ngrid,klev+1)
618      REAL ztla(ngrid,klev)
619      REAL zthl(ngrid,klev)
620
621      REAL zqsatth(ngrid,klev)
622      REAL zqsatenv(ngrid,klev)
623     
624     
625      REAL sigma1(ngrid,klev)
626      REAL sigma2(ngrid,klev)
627      REAL qlth(ngrid,klev)
628      REAL qlenv(ngrid,klev)
629      REAL qltot(ngrid,klev)
630      REAL cth(ngrid,klev) 
631      REAL cenv(ngrid,klev)   
632      REAL ctot(ngrid,klev)
633      REAL rneb(ngrid,klev)
634      REAL t(ngrid,klev)
635      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
636      REAL rdd,cppd,Lv
637      REAL alth,alenv,ath,aenv
638      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
639      REAL Tbef,zdelta,qsatbef,zcor
640      REAL qlbef 
641      REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
642      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
643      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
644      REAL zqs(ngrid), qcloud(ngrid)
645      REAL erf
646
647
648
649      IF (iflag_cloudth_vert.GE.1) THEN
650      CALL cloudth_vert_v3(ngrid,klev,ind2,  &
651     &           ztv,po,zqta,fraca, &
652     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
653     &           ratqs,zqs,t)
654      RETURN
655      ENDIF
656!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
657
658
659!-------------------------------------------------------------------------------
660! Initialisation des variables r?elles
661!-------------------------------------------------------------------------------
662      sigma1(:,:)=0.
663      sigma2(:,:)=0.
664      qlth(:,:)=0.
665      qlenv(:,:)=0. 
666      qltot(:,:)=0.
667      rneb(:,:)=0.
668      qcloud(:)=0.
669      cth(:,:)=0.
670      cenv(:,:)=0.
671      ctot(:,:)=0.
672      qsatmmussig1=0.
673      qsatmmussig2=0.
674      rdd=287.04
675      cppd=1005.7
676      pi=3.14159
677      Lv=2.5e6
678      sqrt2pi=sqrt(2.*pi)
679      sqrt2=sqrt(2.)
680      sqrtpi=sqrt(pi)
681
682
683!-------------------------------------------------------------------------------
684! Cloud fraction in the thermals and standard deviation of the PDFs
685!-------------------------------------------------------------------------------                 
686      do ind1=1,ngrid
687
688      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
689
690      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
691
692      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
693      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
694      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
695      qsatbef=MIN(0.5,qsatbef)
696      zcor=1./(1.-retv*qsatbef)
697      qsatbef=qsatbef*zcor
698      zqsatenv(ind1,ind2)=qsatbef
699
700
701      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
702      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
703      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
704
705!po = qt de l'environnement ET des thermique
706!zqenv = qt environnement
707!zqsatenv = qsat environnement
708!zthl = Tl environnement
709
710
711      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
712      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
713      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
714      qsatbef=MIN(0.5,qsatbef)
715      zcor=1./(1.-retv*qsatbef)
716      qsatbef=qsatbef*zcor
717      zqsatth(ind1,ind2)=qsatbef
718           
719      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
720      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
721      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
722     
723!zqta = qt thermals
724!zqsatth = qsat thermals
725!ztla = Tl thermals
726
727!------------------------------------------------------------------------------
728! s standard deviations
729!------------------------------------------------------------------------------
730
731!     tests
732!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
733!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
734!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
735!     final option
736      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
737      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
738 
739!------------------------------------------------------------------------------
740! Condensed water and cloud cover
741!------------------------------------------------------------------------------
742      xth=sth/(sqrt2*sigma2s)
743      xenv=senv/(sqrt2*sigma1s)
744      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
745      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
746      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
747
748      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
749      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
750      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
751
752      if (ctot(ind1,ind2).lt.1.e-10) then
753      ctot(ind1,ind2)=0.
754      qcloud(ind1)=zqsatenv(ind1,ind2)
755      else
756      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
757      endif
758
759      else  ! Environnement only, follow the if l.110
760     
761      zqenv(ind1)=po(ind1)
762      Tbef=t(ind1,ind2)
763      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
764      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
765      qsatbef=MIN(0.5,qsatbef)
766      zcor=1./(1.-retv*qsatbef)
767      qsatbef=qsatbef*zcor
768      zqsatenv(ind1,ind2)=qsatbef
769
770!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
771      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
772      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
773      aenv=1./(1.+(alenv*Lv/cppd))
774      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
775     
776      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
777
778      xenv=senv/(sqrt2*sigma1s)
779      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
780      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
781
782      if (ctot(ind1,ind2).lt.1.e-3) then
783      ctot(ind1,ind2)=0.
784      qcloud(ind1)=zqsatenv(ind1,ind2)
785      else   
786      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
787      endif
788
789
790      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
791      enddo       ! from the loop on ngrid l.108
792      return
793!     end
794END SUBROUTINE cloudth_v3
795
796
797
798!===========================================================================
799     SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
800     &           ztv,po,zqta,fraca, &
801     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
802     &           ratqs,zqs,t)
803
804!===========================================================================
805! Auteur : Arnaud Octavio Jam (LMD/CNRS)
806! Date : 25 Mai 2010
807! Objet : calcule les valeurs de qc et rneb dans les thermiques
808!===========================================================================
809
810
811      USE ioipsl_getin_p_mod, ONLY : getin_p
812
813      IMPLICIT NONE
814
815#include "YOMCST.h"
816#include "YOETHF.h"
817#include "FCTTRE.h"
818#include "thermcell.h"
819#include "nuage.h"
820     
821      INTEGER itap,ind1,ind2
822      INTEGER ngrid,klev,klon,l,ig
823     
824      REAL ztv(ngrid,klev)
825      REAL po(ngrid)
826      REAL zqenv(ngrid)   
827      REAL zqta(ngrid,klev)
828         
829      REAL fraca(ngrid,klev+1)
830      REAL zpspsk(ngrid,klev)
831      REAL paprs(ngrid,klev+1)
832      REAL ztla(ngrid,klev)
833      REAL zthl(ngrid,klev)
834
835      REAL zqsatth(ngrid,klev)
836      REAL zqsatenv(ngrid,klev)
837     
838     
839      REAL sigma1(ngrid,klev)                                                         
840      REAL sigma2(ngrid,klev)
841      REAL qlth(ngrid,klev)
842      REAL qlenv(ngrid,klev)
843      REAL qltot(ngrid,klev)
844      REAL cth(ngrid,klev) 
845      REAL cenv(ngrid,klev)   
846      REAL ctot(ngrid,klev)
847      REAL rneb(ngrid,klev)
848      REAL t(ngrid,klev)                                                                 
849      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
850      REAL rdd,cppd,Lv
851      REAL alth,alenv,ath,aenv
852      REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
853      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
854      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
855      REAL Tbef,zdelta,qsatbef,zcor
856      REAL qlbef 
857      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
858      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
859      ! (J Jouhaud, JL Dufresne, JB Madeleine)
860      REAL,SAVE :: vert_alpha
861      LOGICAL, SAVE :: firstcall = .TRUE.
862
863      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
864      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
865      REAL zqs(ngrid), qcloud(ngrid)
866      REAL erf
867
868!------------------------------------------------------------------------------
869! Initialize
870!------------------------------------------------------------------------------
871      sigma1(:,:)=0.
872      sigma2(:,:)=0.
873      qlth(:,:)=0.
874      qlenv(:,:)=0. 
875      qltot(:,:)=0.
876      rneb(:,:)=0.
877      qcloud(:)=0.
878      cth(:,:)=0.
879      cenv(:,:)=0.
880      ctot(:,:)=0.
881      qsatmmussig1=0.
882      qsatmmussig2=0.
883      rdd=287.04
884      cppd=1005.7
885      pi=3.14159
886      Lv=2.5e6
887      sqrt2pi=sqrt(2.*pi)
888      sqrt2=sqrt(2.)
889      sqrtpi=sqrt(pi)
890
891      IF (firstcall) THEN
892        vert_alpha=0.5
893        CALL getin_p('cloudth_vert_alpha',vert_alpha)
894        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
895        firstcall=.FALSE.
896      ENDIF
897
898!-------------------------------------------------------------------------------
899! Calcul de la fraction du thermique et des ecart-types des distributions
900!-------------------------------------------------------------------------------                 
901      do ind1=1,ngrid
902
903      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
904
905      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
906
907
908      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
909      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
910      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
911      qsatbef=MIN(0.5,qsatbef)
912      zcor=1./(1.-retv*qsatbef)
913      qsatbef=qsatbef*zcor
914      zqsatenv(ind1,ind2)=qsatbef
915
916
917      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
918      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
919      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
920
921!zqenv = qt environnement
922!zqsatenv = qsat environnement
923!zthl = Tl environnement
924
925
926      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
927      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
928      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
929      qsatbef=MIN(0.5,qsatbef)
930      zcor=1./(1.-retv*qsatbef)
931      qsatbef=qsatbef*zcor
932      zqsatth(ind1,ind2)=qsatbef
933           
934      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
935      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
936      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
937     
938     
939!zqta = qt thermals
940!zqsatth = qsat thermals
941!ztla = Tl thermals
942
943!------------------------------------------------------------------------------
944! s standard deviation
945!------------------------------------------------------------------------------
946
947      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
948      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
949!      tests
950!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
951!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
952!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
953!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
954!       if (paprs(ind1,ind2).gt.90000) then
955!       ratqs(ind1,ind2)=0.002
956!       else
957!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
958!       endif
959!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
960!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
961!       sigma1s=ratqs(ind1,ind2)*po(ind1)
962!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
963 
964       IF (iflag_cloudth_vert == 1) THEN
965!-------------------------------------------------------------------------------
966!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
967!-------------------------------------------------------------------------------
968
969      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
970      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
971
972      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
973      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
974      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
975      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
976      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
977      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
978     
979      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
980      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
981      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
982
983      ! Environment
984      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
985      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
986      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
987      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
988
989      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
990
991      ! Thermal
992      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
993      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
994      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
995      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
996      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
997      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
998
999      ELSE IF (iflag_cloudth_vert == 3) THEN
1000
1001!-------------------------------------------------------------------------------
1002!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1003!-------------------------------------------------------------------------------
1004!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
1005!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
1006!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1007!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1008      deltasenv=aenv*vert_alpha*sigma1s
1009      deltasth=ath*vert_alpha*sigma2s
1010     
1011      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1012      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1013      exp_xenv1 = exp(-1.*xenv1**2)
1014      exp_xenv2 = exp(-1.*xenv2**2)
1015      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1016      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1017      exp_xth1 = exp(-1.*xth1**2)
1018      exp_xth2 = exp(-1.*xth2**2)
1019     
1020      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
1021      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
1022      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1023
1024     
1025      !environnement
1026      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1027      if (deltasenv .lt. 1.e-10) then
1028      qlenv(ind1,ind2)=IntJ
1029      else
1030      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
1031      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
1032      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
1033      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1034      endif
1035
1036
1037      !thermique
1038      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1039      if (deltasth .lt. 1.e-10) then ! Seuil a choisir !!!
1040      qlth(ind1,ind2)=IntJ
1041      else
1042      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1043      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1044      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
1045      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1046      endif
1047
1048      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1049
1050
1051      ENDIF ! of if (iflag_cloudth_vert==1 or 3)
1052
1053      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1054      ctot(ind1,ind2)=0.
1055      qcloud(ind1)=zqsatenv(ind1,ind2)
1056
1057      else
1058               
1059      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1060!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1061!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1062
1063      endif 
1064
1065      else  ! Environment only
1066     
1067      zqenv(ind1)=po(ind1)
1068      Tbef=t(ind1,ind2)
1069      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1070      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1071      qsatbef=MIN(0.5,qsatbef)
1072      zcor=1./(1.-retv*qsatbef)
1073      qsatbef=qsatbef*zcor
1074      zqsatenv(ind1,ind2)=qsatbef
1075     
1076
1077!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1078      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
1079      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
1080      aenv=1./(1.+(alenv*Lv/cppd))
1081      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1082     
1083
1084      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1085
1086      xenv=senv/(sqrt2*sigma1s)
1087      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1088      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
1089     
1090      if (ctot(ind1,ind2).lt.1.e-3) then
1091      ctot(ind1,ind2)=0.
1092      qcloud(ind1)=zqsatenv(ind1,ind2)
1093
1094      else   
1095               
1096      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1097
1098      endif   
1099 
1100      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1101      enddo       ! from the loop on ngrid l.333
1102     
1103      return
1104!     end
1105END SUBROUTINE cloudth_vert_v3
1106!
1107END MODULE cloudth_mod
Note: See TracBrowser for help on using the repository browser.