source: LMDZ5/trunk/libf/phylmd/cloudth.F90 @ 2548

Last change on this file since 2548 was 2547, checked in by jbmadeleine, 8 years ago

Ajout de l'option iflag_cloudth_vert=2 developpe par Jean Jouhaud pour l heterogeneite verticale sous maille des
nuages. Pour le moment, l ecart type de la pdf verticale vaut la moitie de l ecart type de la pdf horizontale.
Le iflag_cloudth_vert=1 conserve l heterogeneite verticale telle qu introduite par Arnaud Jam en Avril 2015.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 19.6 KB
Line 
1       SUBROUTINE cloudth(ngrid,klev,ind2,  &
2     &           ztv,po,zqta,fraca, &
3     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
4     &           ratqs,zqs,t)
5
6
7      USE IOIPSL, ONLY : getin
8      IMPLICIT NONE
9
10
11!===========================================================================
12! Auteur : Arnaud Octavio Jam (LMD/CNRS)
13! Date : 25 Mai 2010
14! Objet : calcule les valeurs de qc et rneb dans les thermiques
15!===========================================================================
16
17
18#include "YOMCST.h"
19#include "YOETHF.h"
20#include "FCTTRE.h"
21#include "thermcell.h"
22#include "nuage.h"
23
24      INTEGER itap,ind1,ind2
25      INTEGER ngrid,klev,klon,l,ig
26     
27      REAL ztv(ngrid,klev)
28      REAL po(ngrid)
29      REAL zqenv(ngrid)   
30      REAL zqta(ngrid,klev)
31         
32      REAL fraca(ngrid,klev+1)
33      REAL zpspsk(ngrid,klev)
34      REAL paprs(ngrid,klev+1)
35      REAL ztla(ngrid,klev)
36      REAL zthl(ngrid,klev)
37
38      REAL zqsatth(ngrid,klev)
39      REAL zqsatenv(ngrid,klev)
40     
41     
42      REAL sigma1(ngrid,klev)
43      REAL sigma2(ngrid,klev)
44      REAL qlth(ngrid,klev)
45      REAL qlenv(ngrid,klev)
46      REAL qltot(ngrid,klev)
47      REAL cth(ngrid,klev) 
48      REAL cenv(ngrid,klev)   
49      REAL ctot(ngrid,klev)
50      REAL rneb(ngrid,klev)
51      REAL t(ngrid,klev)
52      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
53      REAL rdd,cppd,Lv
54      REAL alth,alenv,ath,aenv
55      REAL sth,senv,sigma1s,sigma2s,xth,xenv
56      REAL Tbef,zdelta,qsatbef,zcor
57      REAL alpha,qlbef 
58      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
59     
60      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
61      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
62      REAL zqs(ngrid), qcloud(ngrid)
63      REAL erf
64
65
66
67
68!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
69! Gestion de deux versions de cloudth
70!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
71
72      IF (iflag_cloudth_vert.GE.1) THEN
73      CALL cloudth_vert(ngrid,klev,ind2,  &
74     &           ztv,po,zqta,fraca, &
75     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
76     &           ratqs,zqs,t)
77      RETURN
78      ENDIF
79!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
80
81
82!-------------------------------------------------------------------------------
83! Initialisation des variables r?elles
84!-------------------------------------------------------------------------------
85      sigma1(:,:)=0.
86      sigma2(:,:)=0.
87      qlth(:,:)=0.
88      qlenv(:,:)=0. 
89      qltot(:,:)=0.
90      rneb(:,:)=0.
91      qcloud(:)=0.
92      cth(:,:)=0.
93      cenv(:,:)=0.
94      ctot(:,:)=0.
95      qsatmmussig1=0.
96      qsatmmussig2=0.
97      rdd=287.04
98      cppd=1005.7
99      pi=3.14159
100      Lv=2.5e6
101      sqrt2pi=sqrt(2.*pi)
102
103
104
105!-------------------------------------------------------------------------------
106! Calcul de la fraction du thermique et des ?cart-types des distributions
107!-------------------------------------------------------------------------------                 
108      do ind1=1,ngrid
109
110      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
111
112      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
113
114
115!      zqenv(ind1)=po(ind1)
116      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
117      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
118      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
119      qsatbef=MIN(0.5,qsatbef)
120      zcor=1./(1.-retv*qsatbef)
121      qsatbef=qsatbef*zcor
122      zqsatenv(ind1,ind2)=qsatbef
123
124
125
126
127      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
128      aenv=1./(1.+(alenv*Lv/cppd))
129      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
130
131
132
133
134      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
135      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
136      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
137      qsatbef=MIN(0.5,qsatbef)
138      zcor=1./(1.-retv*qsatbef)
139      qsatbef=qsatbef*zcor
140      zqsatth(ind1,ind2)=qsatbef
141           
142      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
143      ath=1./(1.+(alth*Lv/cppd))
144      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
145     
146     
147
148!------------------------------------------------------------------------------
149! Calcul des ?cart-types pour s
150!------------------------------------------------------------------------------
151
152!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
153!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
154!       if (paprs(ind1,ind2).gt.90000) then
155!       ratqs(ind1,ind2)=0.002
156!       else
157!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
158!       endif
159       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
160       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
161!       sigma1s=ratqs(ind1,ind2)*po(ind1)
162!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
163 
164!------------------------------------------------------------------------------
165! Calcul de l'eau condens?e et de la couverture nuageuse
166!------------------------------------------------------------------------------
167      sqrt2pi=sqrt(2.*pi)
168      xth=sth/(sqrt(2.)*sigma2s)
169      xenv=senv/(sqrt(2.)*sigma1s)
170      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
171      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
172      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
173!      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
174
175
176
177      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
178      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
179      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
180!      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
181     
182
183!      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
184
185!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
186      if (ctot(ind1,ind2).lt.1.e-10) then
187      ctot(ind1,ind2)=0.
188      qcloud(ind1)=zqsatenv(ind1,ind2)
189
190      else   
191               
192      ctot(ind1,ind2)=ctot(ind1,ind2)
193      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
194
195      endif                           
196     
197         
198!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
199
200
201      else  ! gaussienne environnement seule
202     
203      zqenv(ind1)=po(ind1)
204      Tbef=t(ind1,ind2)
205      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
206      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
207      qsatbef=MIN(0.5,qsatbef)
208      zcor=1./(1.-retv*qsatbef)
209      qsatbef=qsatbef*zcor
210      zqsatenv(ind1,ind2)=qsatbef
211     
212
213!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
214      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
215      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
216      aenv=1./(1.+(alenv*Lv/cppd))
217      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
218     
219
220      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
221
222      sqrt2pi=sqrt(2.*pi)
223      xenv=senv/(sqrt(2.)*sigma1s)
224      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
225      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
226     
227      if (ctot(ind1,ind2).lt.1.e-3) then
228      ctot(ind1,ind2)=0.
229      qcloud(ind1)=zqsatenv(ind1,ind2)
230
231      else   
232               
233      ctot(ind1,ind2)=ctot(ind1,ind2)
234      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
235
236      endif   
237 
238 
239 
240 
241 
242 
243      endif   
244      enddo
245     
246      return
247      end
248
249
250
251!===========================================================================
252     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
253     &           ztv,po,zqta,fraca, &
254     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
255     &           ratqs,zqs,t)
256
257
258      IMPLICIT NONE
259
260
261!===========================================================================
262! Auteur : Arnaud Octavio Jam (LMD/CNRS)
263! Date : 25 Mai 2010
264! Objet : calcule les valeurs de qc et rneb dans les thermiques
265!===========================================================================
266
267
268#include "YOMCST.h"
269#include "YOETHF.h"
270#include "FCTTRE.h"
271#include "thermcell.h"
272#include "nuage.h"
273     
274      INTEGER itap,ind1,ind2
275      INTEGER ngrid,klev,klon,l,ig
276     
277      REAL ztv(ngrid,klev)
278      REAL po(ngrid)
279      REAL zqenv(ngrid)   
280      REAL zqta(ngrid,klev)
281         
282      REAL fraca(ngrid,klev+1)
283      REAL zpspsk(ngrid,klev)
284      REAL paprs(ngrid,klev+1)
285      REAL ztla(ngrid,klev)
286      REAL zthl(ngrid,klev)
287
288      REAL zqsatth(ngrid,klev)
289      REAL zqsatenv(ngrid,klev)
290     
291     
292      REAL sigma1(ngrid,klev)                                                         
293      REAL sigma2(ngrid,klev)
294      REAL qlth(ngrid,klev)
295      REAL qlenv(ngrid,klev)
296      REAL qltot(ngrid,klev)
297      REAL cth(ngrid,klev) 
298      REAL cenv(ngrid,klev)   
299      REAL ctot(ngrid,klev)
300      REAL rneb(ngrid,klev)
301      REAL t(ngrid,klev)                                                                 
302      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
303      REAL rdd,cppd,Lv,sqrt2,sqrtpi
304      REAL alth,alenv,ath,aenv
305      REAL sth,senv,sigma1s,sigma2s,xth,xenv
306      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
307      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
308      REAL Tbef,zdelta,qsatbef,zcor
309      REAL alpha,qlbef 
310      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
311     
312      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
313      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
314      REAL zqs(ngrid), qcloud(ngrid)
315      REAL erf
316
317
318
319
320
321!------------------------------------------------------------------------------
322! Initialisation des variables r?elles
323!------------------------------------------------------------------------------
324      sigma1(:,:)=0.
325      sigma2(:,:)=0.
326      qlth(:,:)=0.
327      qlenv(:,:)=0. 
328      qltot(:,:)=0.
329      rneb(:,:)=0.
330      qcloud(:)=0.
331      cth(:,:)=0.
332      cenv(:,:)=0.
333      ctot(:,:)=0.
334      qsatmmussig1=0.
335      qsatmmussig2=0.
336      rdd=287.04
337      cppd=1005.7
338      pi=3.14159
339      Lv=2.5e6
340      sqrt2pi=sqrt(2.*pi)
341      sqrt2=sqrt(2.)
342      sqrtpi=sqrt(pi)
343
344
345
346!-------------------------------------------------------------------------------
347! Calcul de la fraction du thermique et des ?cart-types des distributions
348!-------------------------------------------------------------------------------                 
349      do ind1=1,ngrid
350
351      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
352
353      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
354
355
356!      zqenv(ind1)=po(ind1)
357      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
358      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
359      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
360      qsatbef=MIN(0.5,qsatbef)
361      zcor=1./(1.-retv*qsatbef)
362      qsatbef=qsatbef*zcor
363      zqsatenv(ind1,ind2)=qsatbef
364
365
366
367
368      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
369      aenv=1./(1.+(alenv*Lv/cppd))
370      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
371
372
373
374
375      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
376      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
377      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
378      qsatbef=MIN(0.5,qsatbef)
379      zcor=1./(1.-retv*qsatbef)
380      qsatbef=qsatbef*zcor
381      zqsatth(ind1,ind2)=qsatbef
382           
383      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
384      ath=1./(1.+(alth*Lv/cppd))
385      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
386     
387     
388
389!------------------------------------------------------------------------------
390! Calcul des ?cart-types pour s
391!------------------------------------------------------------------------------
392
393      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
394      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
395!       if (paprs(ind1,ind2).gt.90000) then
396!       ratqs(ind1,ind2)=0.002
397!       else
398!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
399!       endif
400!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
401!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
402!       sigma1s=ratqs(ind1,ind2)*po(ind1)
403!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
404 
405!------------------------------------------------------------------------------
406! Calcul de l'eau condens?e et de la couverture nuageuse
407!------------------------------------------------------------------------------
408      sqrt2pi=sqrt(2.*pi)
409      xth=sth/(sqrt(2.)*sigma2s)
410      xenv=senv/(sqrt(2.)*sigma1s)
411      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
412      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
413      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
414!      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
415
416
417
418      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
419      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
420      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
421!      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
422     
423
424!      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
425
426
427       IF (iflag_cloudth_vert == 1) THEN
428!-------------------------------------------------------------------------------
429!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
430!-------------------------------------------------------------------------------
431!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
432!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
433      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
434      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
435!      deltasenv=aenv*0.01*po(ind1)
436!     deltasth=ath*0.01*zqta(ind1,ind2)   
437      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
438      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
439      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
440      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
441      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
442      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
443     
444      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
445      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
446      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
447
448      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
449      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
450      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
451      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
452
453      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
454!      qlenv(ind1,ind2)=IntJ
455!      print*, qlenv(ind1,ind2),'VERIF EAU'
456
457
458      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
459!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
460!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
461      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
462      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
463      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
464      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
465!      qlth(ind1,ind2)=IntJ
466!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
467      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
468
469      ELSE IF (iflag_cloudth_vert == 2) THEN
470
471!-------------------------------------------------------------------------------
472!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
473!-------------------------------------------------------------------------------
474!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
475!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
476!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
477!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
478      deltasenv=aenv*0.5*sigma1s
479      deltasth=ath*0.5*sigma2s
480     
481      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
482      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
483      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
484      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
485!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
486!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
487     
488      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
489      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
490      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
491
492      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
493      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
494      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
495      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
496
497!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
498!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
499!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
500
501      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
502!      qlenv(ind1,ind2)=IntJ
503!      print*, qlenv(ind1,ind2),'VERIF EAU'
504
505      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
506      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
507      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
508      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
509     
510      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
511!      qlth(ind1,ind2)=IntJ
512!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
513      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
514     
515
516
517
518      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
519
520!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
521
522      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
523      ctot(ind1,ind2)=0.
524      qcloud(ind1)=zqsatenv(ind1,ind2)
525
526      else
527               
528      ctot(ind1,ind2)=ctot(ind1,ind2)
529      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
530!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
531!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
532
533      endif 
534                       
535     
536         
537!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
538
539
540      else  ! gaussienne environnement seule
541     
542      zqenv(ind1)=po(ind1)
543      Tbef=t(ind1,ind2)
544      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
545      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
546      qsatbef=MIN(0.5,qsatbef)
547      zcor=1./(1.-retv*qsatbef)
548      qsatbef=qsatbef*zcor
549      zqsatenv(ind1,ind2)=qsatbef
550     
551
552!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
553      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
554      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
555      aenv=1./(1.+(alenv*Lv/cppd))
556      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
557     
558
559      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
560
561      sqrt2pi=sqrt(2.*pi)
562      xenv=senv/(sqrt(2.)*sigma1s)
563      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
564      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
565     
566      if (ctot(ind1,ind2).lt.1.e-3) then
567      ctot(ind1,ind2)=0.
568      qcloud(ind1)=zqsatenv(ind1,ind2)
569
570      else   
571               
572      ctot(ind1,ind2)=ctot(ind1,ind2)
573      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
574
575      endif   
576 
577 
578 
579 
580 
581 
582      endif   
583      enddo
584     
585      return
586      end
Note: See TracBrowser for help on using the repository browser.