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

Last change on this file since 2639 was 2586, checked in by jbmadeleine, 8 years ago

Added a parameter called cloudth_vert_alpha defined in physiq.def;
it can be used to change the width of the PDF used for the cloud
vertical subgrid scale heterogeneity.

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