source: LMDZ5/branches/AI-cosp/libf/phylmd/cloudth.F90 @ 5444

Last change on this file since 5444 was 2311, checked in by Ehouarn Millour, 10 years ago

Further modifications to enforce physics/dynamics separation:

  • moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
  • created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
  • created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
  • reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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