Context Navigation

cloudth.F90 @ 2267

Last change on this file since 2267 was 2267, checked in by fhourdin, 9 years ago

Modification du modèle du thermique pour diminuer entre autres la
sensibilité au pas de temps. Convergence numérique si iflag_thermals_ed=7.
Ne concerne que la version iflag_thermals_ed=8

Premier essaie de prise en compte d'une variabilité verticale sous maille
de l'eau pour accroître la converture nuageuse.
Sous la clé iflag_coudth_vert=1 (1 par défaut)

Modification of the thermal plume model in case iflag_thermals_ed=8
Modification of the bigaussian scheme for clouds associated with

thermals. Activated if iflag_cloudth_vert=1

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format