Context Navigation

cloudth.F90 @ 2543

Last change on this file since 2543 was 2543, checked in by jbmadeleine, 9 years ago

Ajout de l'option iflag_cloudth_vert=2 developpe par Jean Jouhaud pour l heterogeneite verticale sous maille des nuages.

Pour le moment, l ecart type de la pdf verticale vaut la moitie de l ecart type de la pdf horizontale.

Le iflag_cloudth_vert=1 conserve l heterogeneite verticale telle qu introduite par Arnaud Jam en Avril 2015.

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

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

Download in other formats:

Original Format