source: LMDZ5/trunk/libf/phylmd/cloudth_mod.F90 @ 2945

Last change on this file since 2945 was 2945, checked in by jbmadeleine, 7 years ago
  • Added a new output called rneblsvol which is the cloud fraction by volume

computed in the thermals (see cloudth_vert in cloudth_mod.F90)

  • Added an option called iflag_rain_incloud_vol that computes the conversion

of cloud water to rain using the cloud fraction by volume instead of the cloud
fraction by area, which is larger and otherwise erroneously reduces the in-cloud
water content; iflag_rain_incloud_vol can only be used for iflag_cloudth_vert>=3
File size: 40.2 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,ctot_vol,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      REAL sigma1(ngrid,klev)                                                         
627      REAL sigma2(ngrid,klev)
628      REAL qlth(ngrid,klev)
629      REAL qlenv(ngrid,klev)
630      REAL qltot(ngrid,klev)
631      REAL cth(ngrid,klev)
632      REAL cenv(ngrid,klev)   
633      REAL ctot(ngrid,klev)
634      REAL cth_vol(ngrid,klev)
635      REAL cenv_vol(ngrid,klev)
636      REAL ctot_vol(ngrid,klev)
637      REAL rneb(ngrid,klev)     
638      REAL t(ngrid,klev)
639      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
640      REAL rdd,cppd,Lv
641      REAL alth,alenv,ath,aenv
642      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
643      REAL Tbef,zdelta,qsatbef,zcor
644      REAL qlbef 
645      REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
646      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
647      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
648      REAL zqs(ngrid), qcloud(ngrid)
649      REAL erf
650
651
652
653      IF (iflag_cloudth_vert.GE.1) THEN
654      CALL cloudth_vert_v3(ngrid,klev,ind2,  &
655     &           ztv,po,zqta,fraca, &
656     &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
657     &           ratqs,zqs,t)
658      RETURN
659      ENDIF
660!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
661
662
663!-------------------------------------------------------------------------------
664! Initialisation des variables r?elles
665!-------------------------------------------------------------------------------
666      sigma1(:,:)=0.
667      sigma2(:,:)=0.
668      qlth(:,:)=0.
669      qlenv(:,:)=0. 
670      qltot(:,:)=0.
671      rneb(:,:)=0.
672      qcloud(:)=0.
673      cth(:,:)=0.
674      cenv(:,:)=0.
675      ctot(:,:)=0.
676      cth_vol(:,:)=0.
677      cenv_vol(:,:)=0.
678      ctot_vol(:,:)=0.
679      qsatmmussig1=0.
680      qsatmmussig2=0.
681      rdd=287.04
682      cppd=1005.7
683      pi=3.14159
684      Lv=2.5e6
685      sqrt2pi=sqrt(2.*pi)
686      sqrt2=sqrt(2.)
687      sqrtpi=sqrt(pi)
688
689
690!-------------------------------------------------------------------------------
691! Cloud fraction in the thermals and standard deviation of the PDFs
692!-------------------------------------------------------------------------------                 
693      do ind1=1,ngrid
694
695      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
696
697      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
698
699      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
700      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
701      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
702      qsatbef=MIN(0.5,qsatbef)
703      zcor=1./(1.-retv*qsatbef)
704      qsatbef=qsatbef*zcor
705      zqsatenv(ind1,ind2)=qsatbef
706
707
708      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
709      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
710      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
711
712!po = qt de l'environnement ET des thermique
713!zqenv = qt environnement
714!zqsatenv = qsat environnement
715!zthl = Tl environnement
716
717
718      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
719      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
720      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
721      qsatbef=MIN(0.5,qsatbef)
722      zcor=1./(1.-retv*qsatbef)
723      qsatbef=qsatbef*zcor
724      zqsatth(ind1,ind2)=qsatbef
725           
726      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
727      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
728      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
729     
730!zqta = qt thermals
731!zqsatth = qsat thermals
732!ztla = Tl thermals
733
734!------------------------------------------------------------------------------
735! s standard deviations
736!------------------------------------------------------------------------------
737
738!     tests
739!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
740!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
741!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
742!     final option
743      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
744      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
745 
746!------------------------------------------------------------------------------
747! Condensed water and cloud cover
748!------------------------------------------------------------------------------
749      xth=sth/(sqrt2*sigma2s)
750      xenv=senv/(sqrt2*sigma1s)
751      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
752      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
753      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
754      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
755
756      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
757      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
758      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
759
760      if (ctot(ind1,ind2).lt.1.e-10) then
761      ctot(ind1,ind2)=0.
762      qcloud(ind1)=zqsatenv(ind1,ind2)
763      else
764      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
765      endif
766
767      else  ! Environnement only, follow the if l.110
768     
769      zqenv(ind1)=po(ind1)
770      Tbef=t(ind1,ind2)
771      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
772      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
773      qsatbef=MIN(0.5,qsatbef)
774      zcor=1./(1.-retv*qsatbef)
775      qsatbef=qsatbef*zcor
776      zqsatenv(ind1,ind2)=qsatbef
777
778!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
779      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
780      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
781      aenv=1./(1.+(alenv*Lv/cppd))
782      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
783     
784      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
785
786      xenv=senv/(sqrt2*sigma1s)
787      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
788      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
789      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
790
791      if (ctot(ind1,ind2).lt.1.e-3) then
792      ctot(ind1,ind2)=0.
793      qcloud(ind1)=zqsatenv(ind1,ind2)
794      else   
795      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
796      endif
797
798
799      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
800      enddo       ! from the loop on ngrid l.108
801      return
802!     end
803END SUBROUTINE cloudth_v3
804
805
806
807!===========================================================================
808     SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
809     &           ztv,po,zqta,fraca, &
810     &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
811     &           ratqs,zqs,t)
812
813!===========================================================================
814! Auteur : Arnaud Octavio Jam (LMD/CNRS)
815! Date : 25 Mai 2010
816! Objet : calcule les valeurs de qc et rneb dans les thermiques
817!===========================================================================
818
819
820      USE ioipsl_getin_p_mod, ONLY : getin_p
821
822      IMPLICIT NONE
823
824#include "YOMCST.h"
825#include "YOETHF.h"
826#include "FCTTRE.h"
827#include "thermcell.h"
828#include "nuage.h"
829     
830      INTEGER itap,ind1,ind2
831      INTEGER ngrid,klev,klon,l,ig
832     
833      REAL ztv(ngrid,klev)
834      REAL po(ngrid)
835      REAL zqenv(ngrid)   
836      REAL zqta(ngrid,klev)
837         
838      REAL fraca(ngrid,klev+1)
839      REAL zpspsk(ngrid,klev)
840      REAL paprs(ngrid,klev+1)
841      REAL ztla(ngrid,klev)
842      REAL zthl(ngrid,klev)
843
844      REAL zqsatth(ngrid,klev)
845      REAL zqsatenv(ngrid,klev)
846     
847   
848
849
850      REAL sigma1(ngrid,klev)                                                         
851      REAL sigma2(ngrid,klev)
852      REAL qlth(ngrid,klev)
853      REAL qlenv(ngrid,klev)
854      REAL qltot(ngrid,klev)
855      REAL cth(ngrid,klev)
856      REAL cenv(ngrid,klev)   
857      REAL ctot(ngrid,klev)
858      REAL cth_vol(ngrid,klev)
859      REAL cenv_vol(ngrid,klev)
860      REAL ctot_vol(ngrid,klev)
861      REAL rneb(ngrid,klev)
862      REAL t(ngrid,klev)                                                                 
863      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
864      REAL rdd,cppd,Lv
865      REAL alth,alenv,ath,aenv
866      REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
867      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
868      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
869      REAL Tbef,zdelta,qsatbef,zcor
870      REAL qlbef 
871      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
872      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
873      ! (J Jouhaud, JL Dufresne, JB Madeleine)
874      REAL,SAVE :: vert_alpha, vert_alpha_th
875      !$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
876      LOGICAL, SAVE :: firstcall = .TRUE.
877      !$OMP THREADPRIVATE(firstcall)
878
879      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
880      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
881      REAL zqs(ngrid), qcloud(ngrid)
882      REAL erf
883
884!------------------------------------------------------------------------------
885! Initialize
886!------------------------------------------------------------------------------
887      sigma1(:,:)=0.
888      sigma2(:,:)=0.
889      qlth(:,:)=0.
890      qlenv(:,:)=0. 
891      qltot(:,:)=0.
892      rneb(:,:)=0.
893      qcloud(:)=0.
894      cth(:,:)=0.
895      cenv(:,:)=0.
896      ctot(:,:)=0.
897      cth_vol(:,:)=0.
898      cenv_vol(:,:)=0.
899      ctot_vol(:,:)=0.
900      qsatmmussig1=0.
901      qsatmmussig2=0.
902      rdd=287.04
903      cppd=1005.7
904      pi=3.14159
905      Lv=2.5e6
906      sqrt2pi=sqrt(2.*pi)
907      sqrt2=sqrt(2.)
908      sqrtpi=sqrt(pi)
909
910      IF (firstcall) THEN
911        vert_alpha=0.5
912        CALL getin_p('cloudth_vert_alpha',vert_alpha)
913        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
914        ! The factor used for the thermal is equal to that of the environment
915        !   if nothing is explicitly specified in the def file
916        vert_alpha_th=vert_alpha
917        CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
918        WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th
919        firstcall=.FALSE.
920      ENDIF
921
922!-------------------------------------------------------------------------------
923! Calcul de la fraction du thermique et des ecart-types des distributions
924!-------------------------------------------------------------------------------                 
925      do ind1=1,ngrid
926
927      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
928
929      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
930
931
932      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
933      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
934      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
935      qsatbef=MIN(0.5,qsatbef)
936      zcor=1./(1.-retv*qsatbef)
937      qsatbef=qsatbef*zcor
938      zqsatenv(ind1,ind2)=qsatbef
939
940
941      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
942      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
943      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
944
945!zqenv = qt environnement
946!zqsatenv = qsat environnement
947!zthl = Tl environnement
948
949
950      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
951      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
952      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
953      qsatbef=MIN(0.5,qsatbef)
954      zcor=1./(1.-retv*qsatbef)
955      qsatbef=qsatbef*zcor
956      zqsatth(ind1,ind2)=qsatbef
957           
958      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
959      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
960      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
961     
962     
963!zqta = qt thermals
964!zqsatth = qsat thermals
965!ztla = Tl thermals
966
967!------------------------------------------------------------------------------
968! s standard deviation
969!------------------------------------------------------------------------------
970
971      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
972      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
973!      tests
974!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
975!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
976!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
977!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
978!       if (paprs(ind1,ind2).gt.90000) then
979!       ratqs(ind1,ind2)=0.002
980!       else
981!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
982!       endif
983!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
984!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
985!       sigma1s=ratqs(ind1,ind2)*po(ind1)
986!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
987 
988       IF (iflag_cloudth_vert == 1) THEN
989!-------------------------------------------------------------------------------
990!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
991!-------------------------------------------------------------------------------
992
993      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
994      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
995
996      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
997      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
998      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
999      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
1000      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
1001      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
1002     
1003      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
1004      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
1005      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
1006
1007      ! Environment
1008      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
1009      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
1010      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
1011      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
1012
1013      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1014
1015      ! Thermal
1016      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
1017      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
1018      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
1019      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
1020      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1021      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1022
1023      ELSE IF (iflag_cloudth_vert >= 3) THEN
1024
1025!-------------------------------------------------------------------------------
1026!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1027!-------------------------------------------------------------------------------
1028!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
1029!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
1030!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1031!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1032      IF (iflag_cloudth_vert == 3) THEN
1033        deltasenv=aenv*vert_alpha*sigma1s
1034        deltasth=ath*vert_alpha_th*sigma2s
1035      ELSE IF (iflag_cloudth_vert == 4) THEN
1036        deltasenv=vert_alpha*sigma1s
1037        deltasth=vert_alpha_th*sigma2s
1038      ENDIF
1039     
1040      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1041      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1042      exp_xenv1 = exp(-1.*xenv1**2)
1043      exp_xenv2 = exp(-1.*xenv2**2)
1044      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1045      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1046      exp_xth1 = exp(-1.*xth1**2)
1047      exp_xth2 = exp(-1.*xth2**2)
1048     
1049      !CF_surfacique
1050      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
1051      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
1052      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1053
1054
1055      !CF_volumique & eau condense
1056      !environnement
1057      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1058      IntJ_CF=0.5*(1.-1.*erf(xenv2))
1059      if (deltasenv .lt. 1.e-10) then
1060      qlenv(ind1,ind2)=IntJ
1061      cenv_vol(ind1,ind2)=IntJ_CF
1062      else
1063      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
1064      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
1065      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
1066      IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
1067      IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
1068      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1069      cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1070      endif
1071
1072      !thermique
1073      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1074      IntJ_CF=0.5*(1.-1.*erf(xth2))
1075      if (deltasth .lt. 1.e-10) then
1076      qlth(ind1,ind2)=IntJ
1077      cth_vol(ind1,ind2)=IntJ_CF
1078      else
1079      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1080      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1081      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
1082      IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
1083      IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
1084      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1085      cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1086      endif
1087
1088
1089      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1090      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1091
1092      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1093
1094
1095      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1096      ctot(ind1,ind2)=0.
1097      ctot_vol(ind1,ind2)=0.
1098      qcloud(ind1)=zqsatenv(ind1,ind2)
1099
1100      else
1101               
1102      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1103!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1104!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1105
1106      endif 
1107
1108      else  ! Environment only
1109     
1110      zqenv(ind1)=po(ind1)
1111      Tbef=t(ind1,ind2)
1112      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1113      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1114      qsatbef=MIN(0.5,qsatbef)
1115      zcor=1./(1.-retv*qsatbef)
1116      qsatbef=qsatbef*zcor
1117      zqsatenv(ind1,ind2)=qsatbef
1118     
1119
1120!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1121      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
1122      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
1123      aenv=1./(1.+(alenv*Lv/cppd))
1124      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1125     
1126
1127      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1128
1129      xenv=senv/(sqrt2*sigma1s)
1130      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1131      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1132      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
1133     
1134      if (ctot(ind1,ind2).lt.1.e-3) then
1135      ctot(ind1,ind2)=0.
1136      qcloud(ind1)=zqsatenv(ind1,ind2)
1137
1138      else   
1139               
1140      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1141
1142      endif   
1143 
1144      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1145      enddo       ! from the loop on ngrid l.333
1146     
1147      return
1148!     end
1149END SUBROUTINE cloudth_vert_v3
1150!
1151END MODULE cloudth_mod
Note: See TracBrowser for help on using the repository browser.