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

Last change on this file since 2694 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
RevLine 
[1399]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!===========================================================================
[2662]11! Author : Arnaud Octavio Jam (LMD/CNRS)
[1399]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"
[2547]21#include "nuage.h"
[1399]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     
[2267]41      REAL sigma1(ngrid,klev)
[1399]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)
[2267]50      REAL t(ngrid,klev)
[2662]51      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
[1399]52      REAL rdd,cppd,Lv
53      REAL alth,alenv,ath,aenv
[2662]54      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
[1399]55      REAL Tbef,zdelta,qsatbef,zcor
[2586]56      REAL qlbef 
[2662]57      REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
[1399]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
[2547]65      IF (iflag_cloudth_vert.GE.1) THEN
66      CALL cloudth_vert(ngrid,klev,ind2,  &
[2267]67     &           ztv,po,zqta,fraca, &
68     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
69     &           ratqs,zqs,t)
[2547]70      RETURN
71      ENDIF
[2267]72!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
73
74
75!-------------------------------------------------------------------------------
[2140]76! Initialisation des variables r?elles
[2267]77!-------------------------------------------------------------------------------
[1399]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)
[2662]95      sqrt2=sqrt(2.)
96      sqrtpi=sqrt(pi)
[1399]97
98
[2267]99!-------------------------------------------------------------------------------
[2662]100! Cloud fraction in the thermals and standard deviation of the PDFs
[2267]101!-------------------------------------------------------------------------------                 
[1399]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))
[2662]110      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
[1399]111      qsatbef=MIN(0.5,qsatbef)
112      zcor=1./(1.-retv*qsatbef)
113      qsatbef=qsatbef*zcor
114      zqsatenv(ind1,ind2)=qsatbef
115
116
[2662]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
[1399]120
[2662]121!po = qt de l'environnement ET des thermique
122!zqenv = qt environnement
123!zqsatenv = qsat environnement
124!zthl = Tl environnement
[1399]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           
[2662]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
[1399]138     
[2662]139!zqta = qt thermals
140!zqsatth = qsat thermals
141!ztla = Tl thermals
[1399]142
[2267]143!------------------------------------------------------------------------------
[2662]144! s standard deviations
[2267]145!------------------------------------------------------------------------------
[1399]146
[2662]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)
[1399]154 
[2267]155!------------------------------------------------------------------------------
[2662]156! Condensed water and cloud cover
[2267]157!------------------------------------------------------------------------------
[2662]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)
[1399]163
[2662]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))
[1399]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)
[2662]171      else
[1399]172      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
[2662]173      endif
[1399]174
[2662]175      else  ! Environnement only, follow the if l.110
[1399]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
[2662]186!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
[1399]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))
[2662]190      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
[1399]191     
[1411]192      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
[1399]193
[2662]194      xenv=senv/(sqrt2*sigma1s)
[1399]195      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
[2662]196      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
197
[1399]198      if (ctot(ind1,ind2).lt.1.e-3) then
199      ctot(ind1,ind2)=0.
[2662]200      qcloud(ind1)=zqsatenv(ind1,ind2)
[1399]201      else   
202      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
[2662]203      endif
[1399]204
[2662]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
[1399]208      return
209      end
210
211
212
[2267]213!===========================================================================
214     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
215     &           ztv,po,zqta,fraca, &
216     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
217     &           ratqs,zqs,t)
[1399]218
[2267]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!===========================================================================
[1399]224
225
[2586]226      USE ioipsl_getin_p_mod, ONLY : getin_p
227
228      IMPLICIT NONE
229
[2267]230#include "YOMCST.h"
231#include "YOETHF.h"
232#include "FCTTRE.h"
233#include "thermcell.h"
[2547]234#include "nuage.h"
235     
[2267]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)                                                                 
[2662]264      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
265      REAL rdd,cppd,Lv
[2267]266      REAL alth,alenv,ath,aenv
[2662]267      REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
[2267]268      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
269      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
270      REAL Tbef,zdelta,qsatbef,zcor
[2586]271      REAL qlbef 
[2267]272      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
[2586]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.
[2662]277
[2267]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!------------------------------------------------------------------------------
[2662]284! Initialize
[2267]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
[2586]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
[2267]312
313!-------------------------------------------------------------------------------
[2662]314! Calcul de la fraction du thermique et des ecart-types des distributions
[2267]315!-------------------------------------------------------------------------------                 
316      do ind1=1,ngrid
317
[2662]318      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
[2267]319
[2662]320      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
[2267]321
322
323      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
324      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
[2662]325      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
[2267]326      qsatbef=MIN(0.5,qsatbef)
327      zcor=1./(1.-retv*qsatbef)
328      qsatbef=qsatbef*zcor
329      zqsatenv(ind1,ind2)=qsatbef
330
331
[2662]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
[2267]335
[2662]336!zqenv = qt environnement
337!zqsatenv = qsat environnement
338!zthl = Tl environnement
[2267]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           
[2662]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
[2267]352     
353     
[2662]354!zqta = qt thermals
355!zqsatth = qsat thermals
356!ztla = Tl thermals
[2267]357
358!------------------------------------------------------------------------------
[2662]359! s standard deviation
[2267]360!------------------------------------------------------------------------------
361
[2662]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)
[2267]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 
[2547]379       IF (iflag_cloudth_vert == 1) THEN
[2267]380!-------------------------------------------------------------------------------
[2662]381!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
[2267]382!-------------------------------------------------------------------------------
[2662]383
[2267]384      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
385      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
[2662]386
[2267]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
[2662]398      ! Environment
[2267]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
[2662]406      ! Thermal
[2267]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
[2547]414      ELSE IF (iflag_cloudth_vert == 2) THEN
415
416!-------------------------------------------------------------------------------
[2662]417!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
[2547]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)
[2586]423      deltasenv=aenv*vert_alpha*sigma1s
424      deltasth=ath*vert_alpha*sigma2s
[2547]425     
426      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
427      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
[2662]428      exp_xenv1 = exp(-1.*xenv1**2)
429      exp_xenv2 = exp(-1.*xenv2**2)
[2547]430      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
431      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
[2662]432      exp_xth1 = exp(-1.*xth1**2)
433      exp_xth2 = exp(-1.*xth2**2)
[2547]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
[2662]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
[2547]445      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
[2662]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
[2547]450
451
[2662]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
[2547]457      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
[2662]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)
[2547]460      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
[2662]461      endif
462
[2547]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
[2267]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
[2662]480      else  ! Environment only
[2267]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
[2662]492!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
[2267]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
[2662]501      xenv=senv/(sqrt2*sigma1s)
[2267]502      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
[2662]503      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
[2267]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 
[2662]515      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
516      enddo       ! from the loop on ngrid l.333
[2267]517     
518      return
519      end
Note: See TracBrowser for help on using the repository browser.