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

Last change on this file since 2284 was 2267, checked in by fhourdin, 10 years ago

Modification du modèle du thermique pour diminuer entre autres la
sensibilité au pas de temps. Convergence numérique si iflag_thermals_ed=7.
Ne concerne que la version iflag_thermals_ed=8

Premier essaie de prise en compte d'une variabilité verticale sous maille
de l'eau pour accroître la converture nuageuse.
Sous la clé iflag_coudth_vert=1 (1 par défaut)

Modification of the thermal plume model in case iflag_thermals_ed=8
Modification of the bigaussian scheme for clouds associated with

thermals. Activated if iflag_cloudth_vert=1
  • 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: 17.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 "iniprint.h"
22#include "thermcell.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      REAL, SAVE :: iflag_cloudth_vert, iflag_cloudth_vert_omp=0
66
67
68      LOGICAL, SAVE :: first=.true.
69
70
71
72
73
74!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
75! Astuce pour gérer deux versions de cloudth en attendant
76! de converger sur une version nouvelle
77!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
78      IF (first) THEN
79     !$OMP MASTER
80     CALL getin('iflag_cloudth_vert',iflag_cloudth_vert_omp)
81     !$OMP END MASTER
82     !$OMP BARRIER
83     iflag_cloudth_vert=iflag_cloudth_vert_omp
84      first=.false.
85      ENDIF
86       IF (iflag_cloudth_vert==1) THEN
87       CALL cloudth_vert(ngrid,klev,ind2,  &
88     &           ztv,po,zqta,fraca, &
89     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
90     &           ratqs,zqs,t)
91       RETURN
92       ENDIF
93!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
94
95
96
97!-------------------------------------------------------------------------------
98! Initialisation des variables r?elles
99!-------------------------------------------------------------------------------
100      sigma1(:,:)=0.
101      sigma2(:,:)=0.
102      qlth(:,:)=0.
103      qlenv(:,:)=0. 
104      qltot(:,:)=0.
105      rneb(:,:)=0.
106      qcloud(:)=0.
107      cth(:,:)=0.
108      cenv(:,:)=0.
109      ctot(:,:)=0.
110      qsatmmussig1=0.
111      qsatmmussig2=0.
112      rdd=287.04
113      cppd=1005.7
114      pi=3.14159
115      Lv=2.5e6
116      sqrt2pi=sqrt(2.*pi)
117
118
119
120!-------------------------------------------------------------------------------
121! Calcul de la fraction du thermique et des ?cart-types des distributions
122!-------------------------------------------------------------------------------                 
123      do ind1=1,ngrid
124
125      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
126
127      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
128
129
130!      zqenv(ind1)=po(ind1)
131      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
132      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
133      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
134      qsatbef=MIN(0.5,qsatbef)
135      zcor=1./(1.-retv*qsatbef)
136      qsatbef=qsatbef*zcor
137      zqsatenv(ind1,ind2)=qsatbef
138
139
140
141
142      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
143      aenv=1./(1.+(alenv*Lv/cppd))
144      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
145
146
147
148
149      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
150      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
151      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
152      qsatbef=MIN(0.5,qsatbef)
153      zcor=1./(1.-retv*qsatbef)
154      qsatbef=qsatbef*zcor
155      zqsatth(ind1,ind2)=qsatbef
156           
157      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
158      ath=1./(1.+(alth*Lv/cppd))
159      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
160     
161     
162
163!------------------------------------------------------------------------------
164! Calcul des ?cart-types pour s
165!------------------------------------------------------------------------------
166
167!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
168!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
169!       if (paprs(ind1,ind2).gt.90000) then
170!       ratqs(ind1,ind2)=0.002
171!       else
172!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
173!       endif
174       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
175       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
176!       sigma1s=ratqs(ind1,ind2)*po(ind1)
177!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
178 
179!------------------------------------------------------------------------------
180! Calcul de l'eau condens?e et de la couverture nuageuse
181!------------------------------------------------------------------------------
182      sqrt2pi=sqrt(2.*pi)
183      xth=sth/(sqrt(2.)*sigma2s)
184      xenv=senv/(sqrt(2.)*sigma1s)
185      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
186      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
187      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
188!      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
189
190
191
192      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
193      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
194      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
195!      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
196     
197
198!      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
199
200!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
201      if (ctot(ind1,ind2).lt.1.e-10) then
202      ctot(ind1,ind2)=0.
203      qcloud(ind1)=zqsatenv(ind1,ind2)
204
205      else   
206               
207      ctot(ind1,ind2)=ctot(ind1,ind2)
208      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
209
210      endif                           
211     
212         
213!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
214
215
216      else  ! gaussienne environnement seule
217     
218      zqenv(ind1)=po(ind1)
219      Tbef=t(ind1,ind2)
220      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
221      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
222      qsatbef=MIN(0.5,qsatbef)
223      zcor=1./(1.-retv*qsatbef)
224      qsatbef=qsatbef*zcor
225      zqsatenv(ind1,ind2)=qsatbef
226     
227
228!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
229      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
230      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
231      aenv=1./(1.+(alenv*Lv/cppd))
232      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
233     
234
235      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
236
237      sqrt2pi=sqrt(2.*pi)
238      xenv=senv/(sqrt(2.)*sigma1s)
239      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
240      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
241     
242      if (ctot(ind1,ind2).lt.1.e-3) then
243      ctot(ind1,ind2)=0.
244      qcloud(ind1)=zqsatenv(ind1,ind2)
245
246      else   
247               
248      ctot(ind1,ind2)=ctot(ind1,ind2)
249      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
250
251      endif   
252 
253 
254 
255 
256 
257 
258      endif   
259      enddo
260     
261      return
262      end
263
264
265
266!===========================================================================
267     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
268     &           ztv,po,zqta,fraca, &
269     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
270     &           ratqs,zqs,t)
271
272
273      IMPLICIT NONE
274
275
276!===========================================================================
277! Auteur : Arnaud Octavio Jam (LMD/CNRS)
278! Date : 25 Mai 2010
279! Objet : calcule les valeurs de qc et rneb dans les thermiques
280!===========================================================================
281
282
283#include "YOMCST.h"
284#include "YOETHF.h"
285#include "FCTTRE.h"
286#include "iniprint.h"
287#include "thermcell.h"
288
289      INTEGER itap,ind1,ind2
290      INTEGER ngrid,klev,klon,l,ig
291     
292      REAL ztv(ngrid,klev)
293      REAL po(ngrid)
294      REAL zqenv(ngrid)   
295      REAL zqta(ngrid,klev)
296         
297      REAL fraca(ngrid,klev+1)
298      REAL zpspsk(ngrid,klev)
299      REAL paprs(ngrid,klev+1)
300      REAL ztla(ngrid,klev)
301      REAL zthl(ngrid,klev)
302
303      REAL zqsatth(ngrid,klev)
304      REAL zqsatenv(ngrid,klev)
305     
306     
307      REAL sigma1(ngrid,klev)                                                         
308      REAL sigma2(ngrid,klev)
309      REAL qlth(ngrid,klev)
310      REAL qlenv(ngrid,klev)
311      REAL qltot(ngrid,klev)
312      REAL cth(ngrid,klev) 
313      REAL cenv(ngrid,klev)   
314      REAL ctot(ngrid,klev)
315      REAL rneb(ngrid,klev)
316      REAL t(ngrid,klev)                                                                 
317      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
318      REAL rdd,cppd,Lv,sqrt2,sqrtpi
319      REAL alth,alenv,ath,aenv
320      REAL sth,senv,sigma1s,sigma2s,xth,xenv
321      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
322      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
323      REAL Tbef,zdelta,qsatbef,zcor
324      REAL alpha,qlbef 
325      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
326     
327      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
328      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
329      REAL zqs(ngrid), qcloud(ngrid)
330      REAL erf
331
332
333
334
335
336!------------------------------------------------------------------------------
337! Initialisation des variables r?elles
338!------------------------------------------------------------------------------
339      sigma1(:,:)=0.
340      sigma2(:,:)=0.
341      qlth(:,:)=0.
342      qlenv(:,:)=0. 
343      qltot(:,:)=0.
344      rneb(:,:)=0.
345      qcloud(:)=0.
346      cth(:,:)=0.
347      cenv(:,:)=0.
348      ctot(:,:)=0.
349      qsatmmussig1=0.
350      qsatmmussig2=0.
351      rdd=287.04
352      cppd=1005.7
353      pi=3.14159
354      Lv=2.5e6
355      sqrt2pi=sqrt(2.*pi)
356      sqrt2=sqrt(2.)
357      sqrtpi=sqrt(pi)
358
359
360
361!-------------------------------------------------------------------------------
362! Calcul de la fraction du thermique et des ?cart-types des distributions
363!-------------------------------------------------------------------------------                 
364      do ind1=1,ngrid
365
366      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
367
368      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
369
370
371!      zqenv(ind1)=po(ind1)
372      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
373      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
374      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
375      qsatbef=MIN(0.5,qsatbef)
376      zcor=1./(1.-retv*qsatbef)
377      qsatbef=qsatbef*zcor
378      zqsatenv(ind1,ind2)=qsatbef
379
380
381
382
383      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
384      aenv=1./(1.+(alenv*Lv/cppd))
385      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
386
387
388
389
390      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
391      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
392      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
393      qsatbef=MIN(0.5,qsatbef)
394      zcor=1./(1.-retv*qsatbef)
395      qsatbef=qsatbef*zcor
396      zqsatth(ind1,ind2)=qsatbef
397           
398      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
399      ath=1./(1.+(alth*Lv/cppd))
400      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
401     
402     
403
404!------------------------------------------------------------------------------
405! Calcul des ?cart-types pour s
406!------------------------------------------------------------------------------
407
408      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
409      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
410!       if (paprs(ind1,ind2).gt.90000) then
411!       ratqs(ind1,ind2)=0.002
412!       else
413!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
414!       endif
415!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
416!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
417!       sigma1s=ratqs(ind1,ind2)*po(ind1)
418!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
419 
420!------------------------------------------------------------------------------
421! Calcul de l'eau condens?e et de la couverture nuageuse
422!------------------------------------------------------------------------------
423      sqrt2pi=sqrt(2.*pi)
424      xth=sth/(sqrt(2.)*sigma2s)
425      xenv=senv/(sqrt(2.)*sigma1s)
426      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
427      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
428      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
429!      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
430
431
432
433      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
434      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
435      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
436!      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
437     
438
439!      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
440
441
442!-------------------------------------------------------------------------------
443!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
444!-------------------------------------------------------------------------------
445!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
446!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
447      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
448      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
449!      deltasenv=aenv*0.01*po(ind1)
450!     deltasth=ath*0.01*zqta(ind1,ind2)   
451      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
452      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
453      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
454      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
455      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
456      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
457     
458      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
459      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
460      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
461
462      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
463      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
464      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
465      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
466
467      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
468!      qlenv(ind1,ind2)=IntJ
469!      print*, qlenv(ind1,ind2),'VERIF EAU'
470
471
472      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
473!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
474!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
475      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
476      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
477      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
478      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
479!      qlth(ind1,ind2)=IntJ
480!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
481      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
482
483!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
484      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
485      ctot(ind1,ind2)=0.
486      qcloud(ind1)=zqsatenv(ind1,ind2)
487
488      else
489               
490      ctot(ind1,ind2)=ctot(ind1,ind2)
491      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
492!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
493!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
494
495      endif 
496                       
497     
498         
499!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
500
501
502      else  ! gaussienne environnement seule
503     
504      zqenv(ind1)=po(ind1)
505      Tbef=t(ind1,ind2)
506      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
507      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
508      qsatbef=MIN(0.5,qsatbef)
509      zcor=1./(1.-retv*qsatbef)
510      qsatbef=qsatbef*zcor
511      zqsatenv(ind1,ind2)=qsatbef
512     
513
514!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
515      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
516      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
517      aenv=1./(1.+(alenv*Lv/cppd))
518      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
519     
520
521      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
522
523      sqrt2pi=sqrt(2.*pi)
524      xenv=senv/(sqrt(2.)*sigma1s)
525      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
526      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
527     
528      if (ctot(ind1,ind2).lt.1.e-3) then
529      ctot(ind1,ind2)=0.
530      qcloud(ind1)=zqsatenv(ind1,ind2)
531
532      else   
533               
534      ctot(ind1,ind2)=ctot(ind1,ind2)
535      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
536
537      endif   
538 
539 
540 
541 
542 
543 
544      endif   
545      enddo
546     
547      return
548      end
Note: See TracBrowser for help on using the repository browser.