LMDZ
cloudth.F90
Go to the documentation of this file.
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
cloudth_vert
subroutine cloudth_vert(ngrid, klev, ind2, ztv, po, zqta, fraca, qcloud, ctot, zpspsk, paprs, ztla, zthl, ratqs, zqs, t)
Definition: cloudth.F90:270
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
paprs
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL &zphi geo500!IM on interpole a chaque pas de temps le paprs
Definition: calcul_STDlev.h:144
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
cloudth
subroutine cloudth(ngrid, klev, ind2, ztv, po, zqta, fraca, qcloud, ctot, zpspsk, paprs, ztla, zthl, ratqs, zqs, t)
Definition: cloudth.F90:5