source: LMDZ6/branches/DYNAMICO-conv-GC/libf/phylmd/cloudth_mod.F90 @ 3603

Last change on this file since 3603 was 2960, checked in by fhourdin, 7 years ago

Un lievre leve, sans doute.
Un tout petit seuil sur les largeurs de distributions sous maille
du schéma de nuage, qui produit un gros effet sur le contraste
strato-cumulus / cumulus ...
File size: 42.1 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      USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
822     &                                cloudth_sigmath,cloudth_sigmaenv
823
824      IMPLICIT NONE
825
826#include "YOMCST.h"
827#include "YOETHF.h"
828#include "FCTTRE.h"
829#include "thermcell.h"
830#include "nuage.h"
831     
832      INTEGER itap,ind1,ind2
833      INTEGER ngrid,klev,klon,l,ig
834     
835      REAL ztv(ngrid,klev)
836      REAL po(ngrid)
837      REAL zqenv(ngrid)   
838      REAL zqta(ngrid,klev)
839         
840      REAL fraca(ngrid,klev+1)
841      REAL zpspsk(ngrid,klev)
842      REAL paprs(ngrid,klev+1)
843      REAL ztla(ngrid,klev)
844      REAL zthl(ngrid,klev)
845
846      REAL zqsatth(ngrid,klev)
847      REAL zqsatenv(ngrid,klev)
848     
849      REAL sigma1(ngrid,klev)                                                         
850      REAL sigma2(ngrid,klev)
851      REAL qlth(ngrid,klev)
852      REAL qlenv(ngrid,klev)
853      REAL qltot(ngrid,klev)
854      REAL cth(ngrid,klev)
855      REAL cenv(ngrid,klev)   
856      REAL ctot(ngrid,klev)
857      REAL cth_vol(ngrid,klev)
858      REAL cenv_vol(ngrid,klev)
859      REAL ctot_vol(ngrid,klev)
860      REAL rneb(ngrid,klev)
861      REAL t(ngrid,klev)                                                                 
862      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
863      REAL rdd,cppd,Lv
864      REAL alth,alenv,ath,aenv
865      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
866      REAL 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      REAL,SAVE :: sigma1s_factor=1.1
877      REAL,SAVE :: sigma1s_power=0.6
878      REAL,SAVE :: cloudth_ratqsmin=-1.
879      !$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,cloudth_ratqsmin)
880      INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
881      !$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
882
883      LOGICAL, SAVE :: firstcall = .TRUE.
884      !$OMP THREADPRIVATE(firstcall)
885
886      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
887      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
888      REAL zqs(ngrid), qcloud(ngrid)
889      REAL erf
890
891!------------------------------------------------------------------------------
892! Initialize
893!------------------------------------------------------------------------------
894      sigma1(:,:)=0.
895      sigma2(:,:)=0.
896      qlth(:,:)=0.
897      qlenv(:,:)=0. 
898      qltot(:,:)=0.
899      rneb(:,:)=0.
900      qcloud(:)=0.
901      cth(:,:)=0.
902      cenv(:,:)=0.
903      ctot(:,:)=0.
904      cth_vol(:,:)=0.
905      cenv_vol(:,:)=0.
906      ctot_vol(:,:)=0.
907      qsatmmussig1=0.
908      qsatmmussig2=0.
909      rdd=287.04
910      cppd=1005.7
911      pi=3.14159
912      Lv=2.5e6
913      sqrt2pi=sqrt(2.*pi)
914      sqrt2=sqrt(2.)
915      sqrtpi=sqrt(pi)
916
917      IF (firstcall) THEN
918        vert_alpha=0.5
919        CALL getin_p('cloudth_vert_alpha',vert_alpha)
920        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
921        ! The factor used for the thermal is equal to that of the environment
922        !   if nothing is explicitly specified in the def file
923        vert_alpha_th=vert_alpha
924        CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
925        WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th
926        ! Factor used in the calculation of sigma1s
927        CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor)
928        WRITE(*,*) 'cloudth_sigma1s_factor = ', sigma1s_factor
929        ! Power used in the calculation of sigma1s
930        CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
931        WRITE(*,*) 'cloudth_sigma1s_power = ', sigma1s_power
932        ! Minimum value for the environmental air subgrid water distrib
933        CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
934        WRITE(*,*) 'cloudth_ratqsmin = ', cloudth_ratqsmin
935        ! Remove the dependency to ratqs from the variance of the vertical PDF
936        CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs)
937        WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
938
939        firstcall=.FALSE.
940      ENDIF
941
942!-------------------------------------------------------------------------------
943! Calcul de la fraction du thermique et des ecart-types des distributions
944!-------------------------------------------------------------------------------                 
945      do ind1=1,ngrid
946
947      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
948
949      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
950
951
952      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
953      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
954      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
955      qsatbef=MIN(0.5,qsatbef)
956      zcor=1./(1.-retv*qsatbef)
957      qsatbef=qsatbef*zcor
958      zqsatenv(ind1,ind2)=qsatbef
959
960
961      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
962      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
963      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
964
965!zqenv = qt environnement
966!zqsatenv = qsat environnement
967!zthl = Tl environnement
968
969
970      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
971      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
972      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
973      qsatbef=MIN(0.5,qsatbef)
974      zcor=1./(1.-retv*qsatbef)
975      qsatbef=qsatbef*zcor
976      zqsatth(ind1,ind2)=qsatbef
977           
978      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
979      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
980      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
981     
982     
983!zqta = qt thermals
984!zqsatth = qsat thermals
985!ztla = Tl thermals
986
987!------------------------------------------------------------------------------
988! s standard deviation
989!------------------------------------------------------------------------------
990
991      sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
992     &                (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
993!     sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
994      IF (cloudth_ratqsmin>0.) THEN
995         sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
996      ELSE
997         sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
998      ENDIF
999      sigma1s = sigma1s_fraca + sigma1s_ratqs
1000      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
1001!      tests
1002!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
1003!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
1004!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
1005!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
1006!       if (paprs(ind1,ind2).gt.90000) then
1007!       ratqs(ind1,ind2)=0.002
1008!       else
1009!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
1010!       endif
1011!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
1012!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
1013!       sigma1s=ratqs(ind1,ind2)*po(ind1)
1014!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
1015 
1016       IF (iflag_cloudth_vert == 1) THEN
1017!-------------------------------------------------------------------------------
1018!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
1019!-------------------------------------------------------------------------------
1020
1021      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1022      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1023
1024      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
1025      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
1026      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
1027      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
1028      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
1029      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
1030     
1031      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
1032      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
1033      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
1034
1035      ! Environment
1036      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
1037      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
1038      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
1039      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
1040
1041      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1042
1043      ! Thermal
1044      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
1045      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
1046      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
1047      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
1048      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1049      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1050
1051      ELSE IF (iflag_cloudth_vert >= 3) THEN
1052
1053!-------------------------------------------------------------------------------
1054!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1055!-------------------------------------------------------------------------------
1056!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
1057!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
1058!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1059!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1060      IF (iflag_cloudth_vert == 3) THEN
1061        deltasenv=aenv*vert_alpha*sigma1s
1062        deltasth=ath*vert_alpha_th*sigma2s
1063      ELSE IF (iflag_cloudth_vert == 4) THEN
1064        IF (iflag_cloudth_vert_noratqs == 1) THEN
1065          deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
1066          deltasth=vert_alpha_th*sigma2s
1067        ELSE
1068          deltasenv=vert_alpha*sigma1s
1069          deltasth=vert_alpha_th*sigma2s
1070        ENDIF
1071      ENDIF
1072     
1073      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1074      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1075      exp_xenv1 = exp(-1.*xenv1**2)
1076      exp_xenv2 = exp(-1.*xenv2**2)
1077      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1078      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1079      exp_xth1 = exp(-1.*xth1**2)
1080      exp_xth2 = exp(-1.*xth2**2)
1081     
1082      !CF_surfacique
1083      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
1084      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
1085      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1086
1087
1088      !CF_volumique & eau condense
1089      !environnement
1090      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1091      IntJ_CF=0.5*(1.-1.*erf(xenv2))
1092      if (deltasenv .lt. 1.e-10) then
1093      qlenv(ind1,ind2)=IntJ
1094      cenv_vol(ind1,ind2)=IntJ_CF
1095      else
1096      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
1097      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
1098      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
1099      IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
1100      IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
1101      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1102      cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1103      endif
1104
1105      !thermique
1106      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1107      IntJ_CF=0.5*(1.-1.*erf(xth2))
1108      if (deltasth .lt. 1.e-10) then
1109      qlth(ind1,ind2)=IntJ
1110      cth_vol(ind1,ind2)=IntJ_CF
1111      else
1112      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1113      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1114      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
1115      IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
1116      IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
1117      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1118      cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1119      endif
1120
1121
1122      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1123      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1124
1125      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1126
1127
1128      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1129      ctot(ind1,ind2)=0.
1130      ctot_vol(ind1,ind2)=0.
1131      qcloud(ind1)=zqsatenv(ind1,ind2)
1132
1133      else
1134               
1135      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1136!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1137!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1138
1139      endif 
1140
1141      else  ! Environment only
1142     
1143      zqenv(ind1)=po(ind1)
1144      Tbef=t(ind1,ind2)
1145      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1146      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1147      qsatbef=MIN(0.5,qsatbef)
1148      zcor=1./(1.-retv*qsatbef)
1149      qsatbef=qsatbef*zcor
1150      zqsatenv(ind1,ind2)=qsatbef
1151     
1152
1153!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1154      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
1155      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
1156      aenv=1./(1.+(alenv*Lv/cppd))
1157      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1158      sth=0.
1159
1160      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1161      sigma2s=0.
1162
1163      xenv=senv/(sqrt2*sigma1s)
1164      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1165      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1166      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
1167     
1168      if (ctot(ind1,ind2).lt.1.e-3) then
1169      ctot(ind1,ind2)=0.
1170      qcloud(ind1)=zqsatenv(ind1,ind2)
1171
1172      else   
1173               
1174      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1175
1176      endif   
1177 
1178      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1179      ! Outputs used to check the PDFs
1180      cloudth_senv(ind1,ind2) = senv
1181      cloudth_sth(ind1,ind2) = sth
1182      cloudth_sigmaenv(ind1,ind2) = sigma1s
1183      cloudth_sigmath(ind1,ind2) = sigma2s
1184
1185      enddo       ! from the loop on ngrid l.333
1186     
1187      return
1188!     end
1189END SUBROUTINE cloudth_vert_v3
1190!
1191END MODULE cloudth_mod
Note: See TracBrowser for help on using the repository browser.