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

Last change on this file since 2671 was 2662, checked in by jbmadeleine, 8 years ago

Updated cloudth.F90 following changes by J. Jouhaud

  • improved calculation of the PDF for vertical subgrid scale heterogeneity

(used when iflag_cloudth > 0)

  • use of ratqs in the computation of sigma1s (even when iflag_cloudth.EQ.0)

Sigma1s should have depended on ratqs before; this was probably a remnant of previous tests.
  • 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: 18.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      IMPLICIT NONE
8
9
10!===========================================================================
11! Author : 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,sqrt2,sqrtpi,pi
52      REAL rdd,cppd,Lv
53      REAL alth,alenv,ath,aenv
54      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
55      REAL Tbef,zdelta,qsatbef,zcor
56      REAL qlbef 
57      REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
58      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
59      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
60      REAL zqs(ngrid), qcloud(ngrid)
61      REAL erf
62
63
64
65      IF (iflag_cloudth_vert.GE.1) THEN
66      CALL cloudth_vert(ngrid,klev,ind2,  &
67     &           ztv,po,zqta,fraca, &
68     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
69     &           ratqs,zqs,t)
70      RETURN
71      ENDIF
72!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
73
74
75!-------------------------------------------------------------------------------
76! Initialisation des variables r?elles
77!-------------------------------------------------------------------------------
78      sigma1(:,:)=0.
79      sigma2(:,:)=0.
80      qlth(:,:)=0.
81      qlenv(:,:)=0. 
82      qltot(:,:)=0.
83      rneb(:,:)=0.
84      qcloud(:)=0.
85      cth(:,:)=0.
86      cenv(:,:)=0.
87      ctot(:,:)=0.
88      qsatmmussig1=0.
89      qsatmmussig2=0.
90      rdd=287.04
91      cppd=1005.7
92      pi=3.14159
93      Lv=2.5e6
94      sqrt2pi=sqrt(2.*pi)
95      sqrt2=sqrt(2.)
96      sqrtpi=sqrt(pi)
97
98
99!-------------------------------------------------------------------------------
100! Cloud fraction in the thermals and standard deviation of the PDFs
101!-------------------------------------------------------------------------------                 
102      do ind1=1,ngrid
103
104      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
105
106      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
107
108      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
109      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
110      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
111      qsatbef=MIN(0.5,qsatbef)
112      zcor=1./(1.-retv*qsatbef)
113      qsatbef=qsatbef*zcor
114      zqsatenv(ind1,ind2)=qsatbef
115
116
117      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
118      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
119      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
120
121!po = qt de l'environnement ET des thermique
122!zqenv = qt environnement
123!zqsatenv = qsat environnement
124!zthl = Tl environnement
125
126
127      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
128      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
129      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
130      qsatbef=MIN(0.5,qsatbef)
131      zcor=1./(1.-retv*qsatbef)
132      qsatbef=qsatbef*zcor
133      zqsatth(ind1,ind2)=qsatbef
134           
135      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
136      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
137      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
138     
139!zqta = qt thermals
140!zqsatth = qsat thermals
141!ztla = Tl thermals
142
143!------------------------------------------------------------------------------
144! s standard deviations
145!------------------------------------------------------------------------------
146
147!     tests
148!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
149!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
150!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
151!     final option
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.01)**0.4+0.002*zqta(ind1,ind2)
154 
155!------------------------------------------------------------------------------
156! Condensed water and cloud cover
157!------------------------------------------------------------------------------
158      xth=sth/(sqrt2*sigma2s)
159      xenv=senv/(sqrt2*sigma1s)
160      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
161      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
162      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
163
164      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
165      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
166      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
167
168      if (ctot(ind1,ind2).lt.1.e-10) then
169      ctot(ind1,ind2)=0.
170      qcloud(ind1)=zqsatenv(ind1,ind2)
171      else
172      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
173      endif
174
175      else  ! Environnement only, follow the if l.110
176     
177      zqenv(ind1)=po(ind1)
178      Tbef=t(ind1,ind2)
179      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
180      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
181      qsatbef=MIN(0.5,qsatbef)
182      zcor=1./(1.-retv*qsatbef)
183      qsatbef=qsatbef*zcor
184      zqsatenv(ind1,ind2)=qsatbef
185
186!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
187      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
188      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
189      aenv=1./(1.+(alenv*Lv/cppd))
190      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
191     
192      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
193
194      xenv=senv/(sqrt2*sigma1s)
195      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
196      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
197
198      if (ctot(ind1,ind2).lt.1.e-3) then
199      ctot(ind1,ind2)=0.
200      qcloud(ind1)=zqsatenv(ind1,ind2)
201      else   
202      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
203      endif
204
205
206      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
207      enddo       ! from the loop on ngrid l.108
208      return
209      end
210
211
212
213!===========================================================================
214     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
215     &           ztv,po,zqta,fraca, &
216     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
217     &           ratqs,zqs,t)
218
219!===========================================================================
220! Auteur : Arnaud Octavio Jam (LMD/CNRS)
221! Date : 25 Mai 2010
222! Objet : calcule les valeurs de qc et rneb dans les thermiques
223!===========================================================================
224
225
226      USE ioipsl_getin_p_mod, ONLY : getin_p
227
228      IMPLICIT NONE
229
230#include "YOMCST.h"
231#include "YOETHF.h"
232#include "FCTTRE.h"
233#include "thermcell.h"
234#include "nuage.h"
235     
236      INTEGER itap,ind1,ind2
237      INTEGER ngrid,klev,klon,l,ig
238     
239      REAL ztv(ngrid,klev)
240      REAL po(ngrid)
241      REAL zqenv(ngrid)   
242      REAL zqta(ngrid,klev)
243         
244      REAL fraca(ngrid,klev+1)
245      REAL zpspsk(ngrid,klev)
246      REAL paprs(ngrid,klev+1)
247      REAL ztla(ngrid,klev)
248      REAL zthl(ngrid,klev)
249
250      REAL zqsatth(ngrid,klev)
251      REAL zqsatenv(ngrid,klev)
252     
253     
254      REAL sigma1(ngrid,klev)                                                         
255      REAL sigma2(ngrid,klev)
256      REAL qlth(ngrid,klev)
257      REAL qlenv(ngrid,klev)
258      REAL qltot(ngrid,klev)
259      REAL cth(ngrid,klev) 
260      REAL cenv(ngrid,klev)   
261      REAL ctot(ngrid,klev)
262      REAL rneb(ngrid,klev)
263      REAL t(ngrid,klev)                                                                 
264      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
265      REAL rdd,cppd,Lv
266      REAL alth,alenv,ath,aenv
267      REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
268      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
269      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
270      REAL Tbef,zdelta,qsatbef,zcor
271      REAL qlbef 
272      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
273      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
274      ! (J Jouhaud, JL Dufresne, JB Madeleine)
275      REAL,SAVE :: vert_alpha
276      LOGICAL, SAVE :: firstcall = .TRUE.
277
278      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
279      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
280      REAL zqs(ngrid), qcloud(ngrid)
281      REAL erf
282
283!------------------------------------------------------------------------------
284! Initialize
285!------------------------------------------------------------------------------
286      sigma1(:,:)=0.
287      sigma2(:,:)=0.
288      qlth(:,:)=0.
289      qlenv(:,:)=0. 
290      qltot(:,:)=0.
291      rneb(:,:)=0.
292      qcloud(:)=0.
293      cth(:,:)=0.
294      cenv(:,:)=0.
295      ctot(:,:)=0.
296      qsatmmussig1=0.
297      qsatmmussig2=0.
298      rdd=287.04
299      cppd=1005.7
300      pi=3.14159
301      Lv=2.5e6
302      sqrt2pi=sqrt(2.*pi)
303      sqrt2=sqrt(2.)
304      sqrtpi=sqrt(pi)
305
306      IF (firstcall) THEN
307        vert_alpha=0.5
308        CALL getin_p('cloudth_vert_alpha',vert_alpha)
309        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
310        firstcall=.FALSE.
311      ENDIF
312
313!-------------------------------------------------------------------------------
314! Calcul de la fraction du thermique et des ecart-types des distributions
315!-------------------------------------------------------------------------------                 
316      do ind1=1,ngrid
317
318      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
319
320      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
321
322
323      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
324      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
325      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
326      qsatbef=MIN(0.5,qsatbef)
327      zcor=1./(1.-retv*qsatbef)
328      qsatbef=qsatbef*zcor
329      zqsatenv(ind1,ind2)=qsatbef
330
331
332      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
333      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
334      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
335
336!zqenv = qt environnement
337!zqsatenv = qsat environnement
338!zthl = Tl environnement
339
340
341      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
342      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
343      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
344      qsatbef=MIN(0.5,qsatbef)
345      zcor=1./(1.-retv*qsatbef)
346      qsatbef=qsatbef*zcor
347      zqsatth(ind1,ind2)=qsatbef
348           
349      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
350      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
351      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
352     
353     
354!zqta = qt thermals
355!zqsatth = qsat thermals
356!ztla = Tl thermals
357
358!------------------------------------------------------------------------------
359! s standard deviation
360!------------------------------------------------------------------------------
361
362      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
363      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
364!      tests
365!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
366!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
367!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
368!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
369!       if (paprs(ind1,ind2).gt.90000) then
370!       ratqs(ind1,ind2)=0.002
371!       else
372!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
373!       endif
374!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
375!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
376!       sigma1s=ratqs(ind1,ind2)*po(ind1)
377!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
378 
379       IF (iflag_cloudth_vert == 1) THEN
380!-------------------------------------------------------------------------------
381!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
382!-------------------------------------------------------------------------------
383
384      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
385      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
386
387      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
388      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
389      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
390      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
391      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
392      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
393     
394      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
395      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
396      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
397
398      ! Environment
399      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
400      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
401      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
402      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
403
404      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
405
406      ! Thermal
407      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
408      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
409      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
410      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
411      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
412      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
413
414      ELSE IF (iflag_cloudth_vert == 2) THEN
415
416!-------------------------------------------------------------------------------
417!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
418!-------------------------------------------------------------------------------
419!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
420!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
421!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
422!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
423      deltasenv=aenv*vert_alpha*sigma1s
424      deltasth=ath*vert_alpha*sigma2s
425     
426      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
427      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
428      exp_xenv1 = exp(-1.*xenv1**2)
429      exp_xenv2 = exp(-1.*xenv2**2)
430      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
431      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
432      exp_xth1 = exp(-1.*xth1**2)
433      exp_xth2 = exp(-1.*xth2**2)
434     
435      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
436      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
437      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
438
439     
440      !environnement
441      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
442      if (deltasenv .lt. 1.e-10) then
443      qlenv(ind1,ind2)=IntJ
444      else
445      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
446      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
447      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
448      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
449      endif
450
451
452      !thermique
453      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
454      if (deltasth .lt. 1.e-10) then ! Seuil a choisir !!!
455      qlth(ind1,ind2)=IntJ
456      else
457      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
458      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
459      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
460      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
461      endif
462
463      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
464
465
466      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
467
468      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
469      ctot(ind1,ind2)=0.
470      qcloud(ind1)=zqsatenv(ind1,ind2)
471
472      else
473               
474      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
475!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
476!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
477
478      endif 
479
480      else  ! Environment only
481     
482      zqenv(ind1)=po(ind1)
483      Tbef=t(ind1,ind2)
484      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
485      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
486      qsatbef=MIN(0.5,qsatbef)
487      zcor=1./(1.-retv*qsatbef)
488      qsatbef=qsatbef*zcor
489      zqsatenv(ind1,ind2)=qsatbef
490     
491
492!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
493      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
494      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
495      aenv=1./(1.+(alenv*Lv/cppd))
496      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
497     
498
499      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
500
501      xenv=senv/(sqrt2*sigma1s)
502      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
503      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
504     
505      if (ctot(ind1,ind2).lt.1.e-3) then
506      ctot(ind1,ind2)=0.
507      qcloud(ind1)=zqsatenv(ind1,ind2)
508
509      else   
510               
511      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
512
513      endif   
514 
515      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
516      enddo       ! from the loop on ngrid l.333
517     
518      return
519      end
Note: See TracBrowser for help on using the repository browser.