Context Navigation

thermcell_plume.F90 @ 1024

Last change on this file since 1024 was 972, checked in by lmdzadmin, 16 years ago
Version thermique FH/CRio Ajout tests cas physiques non pris en comptes et ajout/enleve prints Nouvelle routine thermcell_flux2.F90 IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 21.1 KB

Line
1	SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
2	& zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star, &
3	& lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, &
4	& ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
5
6	!--------------------------------------------------------------------------
7	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
8	!--------------------------------------------------------------------------
9
10	IMPLICIT NONE
11
12	#include "YOMCST.h"
13	#include "YOETHF.h"
14	#include "FCTTRE.h"
15	#include "iniprint.h"
16
17	INTEGER itap
18
19	INTEGER ngrid,klev
20	REAL ptimestep
21	REAL ztv(ngrid,klev)
22	REAL zthl(ngrid,klev)
23	REAL po(ngrid,klev)
24	REAL zl(ngrid,klev)
25	REAL rhobarz(ngrid,klev)
26	REAL zlev(ngrid,klev+1)
27	REAL pplev(ngrid,klev+1)
28	REAL pphi(ngrid,klev)
29	REAL zpspsk(ngrid,klev)
30	REAL alim_star(ngrid,klev)
31	REAL zmax_sec(ngrid)
32	REAL f0(ngrid)
33	REAL l_mix
34	REAL r_aspect
35	INTEGER lalim(ngrid)
36	integer lev_out ! niveau pour les print
37	real zcon2(ngrid)
38
39	REAL ztva(ngrid,klev)
40	REAL ztla(ngrid,klev)
41	REAL zqla(ngrid,klev)
42	REAL zqla0(ngrid,klev)
43	REAL zqta(ngrid,klev)
44	REAL zha(ngrid,klev)
45
46	REAL detr_star(ngrid,klev)
47	REAL coefc
48	REAL detr_stara(ngrid,klev)
49	REAL detr_starb(ngrid,klev)
50	REAL detr_starc(ngrid,klev)
51	REAL detr_star0(ngrid,klev)
52	REAL detr_star1(ngrid,klev)
53	REAL detr_star2(ngrid,klev)
54
55	REAL entr_star(ngrid,klev)
56	REAL entr_star1(ngrid,klev)
57	REAL entr_star2(ngrid,klev)
58	REAL detr(ngrid,klev)
59	REAL entr(ngrid,klev)
60
61	REAL zw2(ngrid,klev+1)
62	REAL w_est(ngrid,klev+1)
63	REAL f_star(ngrid,klev+1)
64	REAL wa_moy(ngrid,klev+1)
65
66	REAL ztva_est(ngrid,klev)
67	REAL zqla_est(ngrid,klev)
68	REAL zqsatth(ngrid,klev)
69
70	REAL linter(ngrid)
71	INTEGER lmix(ngrid)
72	INTEGER lmix_bis(ngrid)
73	REAL wmaxa(ngrid)
74
75	INTEGER ig,l,k
76
77	real zcor,zdelta,zcvm5,qlbef
78	real Tbef,qsatbef
79	real dqsat_dT,DT,num,denom
80	REAL REPS,RLvCp,DDT0
81	PARAMETER (DDT0=.01)
82	logical Zsat
83
84	Zsat=.false.
85	! Initialisation
86	RLvCp = RLVTT/RCPD
87
88	do l=1,klev
89	do ig=1,ngrid
90	zqla_est(ig,l)=0.
91	ztva_est(ig,l)=ztva(ig,l)
92	zqsatth(ig,l)=0.
93	enddo
94	enddo
95
96	!AM:initialisations du thermique
97	do k=1,klev
98	do ig=1,ngrid
99	ztva(ig,k)=ztv(ig,k)
100	ztla(ig,k)=zthl(ig,k)
101	zqla(ig,k)=0.
102	zqta(ig,k)=po(ig,k)
103	!
104	ztva(ig,k) = ztla(ig,k)zpspsk(ig,k)+RLvCpzqla(ig,k)
105	ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
106	zha(ig,k) = ztva(ig,k)
107	!
108	enddo
109	enddo
110	do k=1,klev
111	do ig=1,ngrid
112	detr_star(ig,k)=0.
113	entr_star(ig,k)=0.
114
115	detr_stara(ig,k)=0.
116	detr_starb(ig,k)=0.
117	detr_starc(ig,k)=0.
118	detr_star0(ig,k)=0.
119	zqla0(ig,k)=0.
120	detr_star1(ig,k)=0.
121	detr_star2(ig,k)=0.
122	entr_star1(ig,k)=0.
123	entr_star2(ig,k)=0.
124
125	detr(ig,k)=0.
126	entr(ig,k)=0.
127	enddo
128	enddo
129	if (prt_level.ge.1) print*,'7 OK convect8'
130	do k=1,klev+1
131	do ig=1,ngrid
132	zw2(ig,k)=0.
133	w_est(ig,k)=0.
134	f_star(ig,k)=0.
135	wa_moy(ig,k)=0.
136	enddo
137	enddo
138
139	if (prt_level.ge.1) print*,'8 OK convect8'
140	do ig=1,ngrid
141	linter(ig)=1.
142	lmix(ig)=1
143	wmaxa(ig)=0.
144	enddo
145
146	!-----------------------------------------------------------------------------------
147	!boucle de calcul de la vitesse verticale dans le thermique
148	!-----------------------------------------------------------------------------------
149	do l=1,klev-1
150	do ig=1,ngrid
151
152
153
154	! Calcul dans la premiere couche active du thermique (ce qu'on teste
155	! en disant que la couche est instable et que w2 en bas de la couche
156	! est nulle.
157
158	if (ztv(ig,l).gt.ztv(ig,l+1) &
159	& .and.alim_star(ig,l).gt.1.e-10 &
160	& .and.zw2(ig,l).lt.1e-10) then
161
162
163	! Le panache va prendre au debut les caracteristiques de l'air contenu
164	! dans cette couche.
165	ztla(ig,l)=zthl(ig,l)
166	zqta(ig,l)=po(ig,l)
167	zqla(ig,l)=zl(ig,l)
168	f_star(ig,l+1)=alim_star(ig,l)
169
170	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) &
171	& *(zlev(ig,l+1)-zlev(ig,l)) &
172	& 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
173	w_est(ig,l+1)=zw2(ig,l+1)
174	!
175
176
177	else if ((zw2(ig,l).ge.1e-10).and. &
178	& (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
179	!estimation du detrainement a partir de la geometrie du pas precedent
180	!tests sur la definition du detr
181	!calcul de detr_star et entr_star
182
183
184
185	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
186	! FH le test miraculeux de Catherine ? Le bout du tunel ?
187	! w_est(ig,3)=zw2(ig,2)* &
188	! & ((f_star(ig,2))**2) &
189	! & /(f_star(ig,2)+alim_star(ig,2))**2+ &
190	! & 2.RG(ztva(ig,1)-ztv(ig,2))/ztv(ig,2) &
191	! & *(zlev(ig,3)-zlev(ig,2))
192	! if (l.gt.2) then
193	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
194
195
196
197	! Premier calcul de la vitesse verticale a partir de la temperature
198	! potentielle virtuelle
199
200	! FH CESTQUOI CA ????
201	#define int1d2
202	!#undef int1d2
203	#ifdef int1d2
204	if (l.ge.2) then
205	#else
206	if (l.gt.2) then
207	#endif
208
209	if (1.eq.1) then
210	w_est(ig,3)=zw2(ig,2)* &
211	& ((f_star(ig,2))**2) &
212	& /(f_star(ig,2)+alim_star(ig,2))**2+ &
213	& 2.RG(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
214	& *(zlev(ig,3)-zlev(ig,2))
215	endif
216
217
218	!---------------------------------------------------------------------------
219	!calcul de l entrainement et du detrainement lateral
220	!---------------------------------------------------------------------------
221	!
222	!test:estimation de ztva_new_est sans entrainement
223
224	Tbef=ztla(ig,l-1)*zpspsk(ig,l)
225	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
226	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
227	qsatbef=MIN(0.5,qsatbef)
228	zcor=1./(1.-retv*qsatbef)
229	qsatbef=qsatbef*zcor
230	Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
231	if (Zsat) then
232	qlbef=max(0.,zqta(ig,l-1)-qsatbef)
233	DT = 0.5RLvCpqlbef
234	do while (abs(DT).gt.DDT0)
235	Tbef=Tbef+DT
236	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
237	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
238	qsatbef=MIN(0.5,qsatbef)
239	zcor=1./(1.-retv*qsatbef)
240	qsatbef=qsatbef*zcor
241	qlbef=zqta(ig,l-1)-qsatbef
242
243	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
244	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
245	zcor=1./(1.-retv*qsatbef)
246	dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
247	num=-Tbef+ztla(ig,l-1)zpspsk(ig,l)+RLvCpqlbef
248	denom=1.+RLvCp*dqsat_dT
249	DT=num/denom
250	enddo
251	zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef)
252	endif
253	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
254	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
255	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
256	& -zqla_est(ig,l))-zqla_est(ig,l))
257
258	w_est(ig,l+1)=zw2(ig,l)* &
259	& ((f_star(ig,l))**2) &
260	& /(f_star(ig,l)+alim_star(ig,l))**2+ &
261	& 2.RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
262	& *(zlev(ig,l+1)-zlev(ig,l))
263	if (w_est(ig,l+1).lt.0.) then
264	w_est(ig,l+1)=zw2(ig,l)
265	endif
266	!
267	!calcul du detrainement
268	!=======================
269
270
271	! Modifications du calcul du detrainement.
272	! Dans la version de la these de Catherine, on passe brusquement
273	! de la version seche a la version nuageuse pour le detrainement
274	! ce qui peut occasioner des oscillations.
275	! dans la nouvelle version, on commence par calculer un detrainement sec.
276	! Puis un autre en cas de nuages.
277	! Puis on combine les deux lineairement en fonction de la quantite d'eau.
278
279	#define int1d3
280	!#undef int1d3
281	#define RIO_TH
282	#ifdef RIO_TH
283	!1. Cas non nuageux
284	! 1.1 on est sous le zmax_sec et w croit
285	if ((w_est(ig,l+1).gt.w_est(ig,l)).and. &
286	& (zlev(ig,l+1).lt.zmax_sec(ig)).and. &
287	#ifdef int1d3
288	& (zqla_est(ig,l).lt.1.e-10)) then
289	#else
290	& (zqla(ig,l-1).lt.1.e-10)) then
291	#endif
292	detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1) &
293	& sqrt(w_est(ig,l+1))sqrt(l_mix*zlev(ig,l+1)) &
294	& -rhobarz(ig,l)sqrt(w_est(ig,l))sqrt(l_mix*zlev(ig,l))) &
295	& /(r_aspect*zmax_sec(ig)))
296	detr_stara(ig,l)=detr_star(ig,l)
297
298	if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
299
300	! 1.2 on est sous le zmax_sec et w decroit
301	else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and. &
302	#ifdef int1d3
303	& (zqla_est(ig,l).lt.1.e-10)) then
304	#else
305	& (zqla(ig,l-1).lt.1.e-10)) then
306	#endif
307	detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig)) &
308	& /(rhobarz(ig,lmix(ig))wmaxa(ig)) &
309	& (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1)) &
310	& *((zmax_sec(ig)-zlev(ig,l+1))/ &
311	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2. &
312	& -rhobarz(ig,l)*sqrt(w_est(ig,l)) &
313	& *((zmax_sec(ig)-zlev(ig,l))/ &
314	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
315	detr_starb(ig,l)=detr_star(ig,l)
316
317	if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
318
319	else
320
321	! 1.3 dans les autres cas
322	detr_star(ig,l)=0.002f0(ig)f_star(ig,l) &
323	& *(zlev(ig,l+1)-zlev(ig,l))
324	detr_starc(ig,l)=detr_star(ig,l)
325
326	if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
327
328	endif
329
330	#else
331
332	! 1.1 on est sous le zmax_sec et w croit
333	if ((w_est(ig,l+1).gt.w_est(ig,l)).and. &
334	& (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
335	detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1) &
336	& sqrt(w_est(ig,l+1))sqrt(l_mix*zlev(ig,l+1)) &
337	& -rhobarz(ig,l)sqrt(w_est(ig,l))sqrt(l_mix*zlev(ig,l))) &
338	& /(r_aspect*zmax_sec(ig)))
339
340	if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
341
342	! 1.2 on est sous le zmax_sec et w decroit
343	else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
344	detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig)) &
345	& /(rhobarz(ig,lmix(ig))wmaxa(ig)) &
346	& (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1)) &
347	& *((zmax_sec(ig)-zlev(ig,l+1))/ &
348	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2. &
349	& -rhobarz(ig,l)*sqrt(w_est(ig,l)) &
350	& *((zmax_sec(ig)-zlev(ig,l))/ &
351	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
352	if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
353
354	else
355	detr_star=0.
356	endif
357
358	! 1.3 dans les autres cas
359	detr_starc(ig,l)=0.002f0(ig)f_star(ig,l) &
360	& *(zlev(ig,l+1)-zlev(ig,l))
361
362	coefc=min(zqla(ig,l-1)/1.e-3,1.)
363	if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
364	coefc=1.
365	! il semble qu'il soit important de baser le calcul sur
366	! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
367	detr_star(ig,l)=detr_starc(ig,l)coefc+detr_star(ig,l)(1.-coefc)
368
369	if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
370
371	#endif
372
373
374	if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
375	!IM 730508 beg
376	! if(itap.GE.7200) THEN
377	! print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
378	! endif
379	!IM 730508 end
380
381	zqla0(ig,l)=zqla_est(ig,l)
382	detr_star0(ig,l)=detr_star(ig,l)
383	!IM 060508 beg
384	! if(detr_star(ig,l).GT.1.) THEN
385	! print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
386	! & ,detr_starc(ig,l),coefc
387	! endif
388	!IM 060508 end
389	!IM 160508 beg
390	!IM 160508 IF (f0(ig).NE.0.) THEN
391	detr_star(ig,l)=detr_star(ig,l)/f0(ig)
392	!IM 160508 ELSE IF(detr_star(ig,l).EQ.0.) THEN
393	!IM 160508 print*,'WARNING1 : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
394	!IM 160508 ELSE
395	!IM 160508 print*,'WARNING2 : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
396	!IM 160508 ENDIF
397	!IM 160508 end
398	!IM 060508 beg
399	! if(detr_star(ig,l).GT.1.) THEN
400	! print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
401	! & float(1)/f0(ig)
402	! endif
403	!IM 060508 end
404	if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
405	!
406	!calcul de entr_star
407
408	! #undef test2
409	! #ifdef test2
410	! La version test2 destabilise beaucoup le modele.
411	! Il semble donc que ca aide d'avoir un entrainement important sous
412	! le nuage.
413	! if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
414	! entr_star(ig,l)=0.4*detr_star(ig,l)
415	! else
416	! entr_star(ig,l)=0.
417	! endif
418	! #else
419	!
420	! Deplacement du calcul de entr_star pour eviter d'avoir aussi
421	! entr_star > fstar.
422	! Redeplacer suite a la transformation du cas detr>f
423	! FH
424
425	if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
426	#define int1d
427	!FH 070508 #define int1d4
428	!#undef int1d4
429	! L'option int1d4 correspond au choix dans le cas ou le detrainement
430	! devient trop grand.
431
432	#ifdef int1d
433
434	#ifdef int1d4
435	#else
436	detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
437	!FH 070508 plus
438	detr_star(ig,l)=min(detr_star(ig,l),1.)
439	#endif
440
441	entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
442
443	if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
444	#ifdef int1d4
445	! Si le detrainement excede le flux en bas + l'entrainement, le thermique
446	! doit disparaitre.
447	if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
448	detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
449	f_star(ig,l+1)=0.
450	linter(ig)=l+1
451	zw2(ig,l+1)=-1.e-10
452	endif
453	#endif
454
455
456	#else
457
458	if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
459	if(l.gt.lalim(ig)) then
460	entr_star(ig,l)=0.4*detr_star(ig,l)
461	else
462
463	! FH :
464	! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
465	! en haut de la couche d'alimentation.
466	! A remettre en questoin a la premiere occasion mais ca peut aider a
467	! ecrire un code robuste.
468	! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
469	! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
470	! d'alimentation, ce qui n'est pas forcement heureux.
471
472	if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
473	#undef pre_int1c
474	#ifdef pre_int1c
475	entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
476	detr_star(ig,l)=entr_star(ig,l)
477	#else
478	entr_star(ig,l)=0.
479	#endif
480
481	endif
482
483	#endif
484
485	if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
486	entr_star1(ig,l)=entr_star(ig,l)
487	detr_star1(ig,l)=detr_star(ig,l)
488	!
489
490	#ifdef int1d
491	#else
492	if (detr_star(ig,l).gt.f_star(ig,l)) then
493
494	! Ce test est là entre autres parce qu'on passe par des valeurs
495	! delirantes de detr_star.
496	! ca vaut sans doute le coup de verifier pourquoi.
497
498	detr_star(ig,l)=f_star(ig,l)
499	#ifdef pre_int1c
500	if (l.gt.lalim(ig)+1) then
501	entr_star(ig,l)=0.
502	alim_star(ig,l)=0.
503	! FH ajout pour forcer a stoper le thermique juste sous le sommet
504	! de la couche (voir calcul de finter)
505	zw2(ig,l+1)=-1.e-10
506	linter(ig)=l+1
507	else
508	entr_star(ig,l)=0.4*detr_star(ig,l)
509	endif
510	#else
511	entr_star(ig,l)=0.4*detr_star(ig,l)
512	#endif
513	endif
514	#endif
515
516	else !l > 2
517	detr_star(ig,l)=0.
518	entr_star(ig,l)=0.
519	endif
520
521	entr_star2(ig,l)=entr_star(ig,l)
522	detr_star2(ig,l)=detr_star(ig,l)
523	if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l
524
525	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
526	! FH inutile si on conserve comme on l'a fait plus haut entr=detr
527	! dans la couche d'alimentation
528	!pas d entrainement dans la couche alim
529	! if ((l.le.lalim(ig))) then
530	! entr_star(ig,l)=0.
531	! endif
532	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
533	!
534	!prise en compte du detrainement et de l entrainement dans le calcul du flux
535	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
536	& -detr_star(ig,l)
537
538	!test sur le signe de f_star
539	if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
540	if (f_star(ig,l+1).gt.1.e-10) then
541	!----------------------------------------------------------------------------
542	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
543	!---------------------------------------------------------------------------
544	!
545	Zsat=.false.
546	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
547	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
548	& /(f_star(ig,l+1)+detr_star(ig,l))
549	!
550	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
551	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
552	& /(f_star(ig,l+1)+detr_star(ig,l))
553	!
554	Tbef=ztla(ig,l)*zpspsk(ig,l)
555	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
556	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
557	qsatbef=MIN(0.5,qsatbef)
558	zcor=1./(1.-retv*qsatbef)
559	qsatbef=qsatbef*zcor
560	Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
561	if (Zsat) then
562	qlbef=max(0.,zqta(ig,l)-qsatbef)
563	DT = 0.5RLvCpqlbef
564	do while (abs(DT).gt.DDT0)
565	Tbef=Tbef+DT
566	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
567	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
568	qsatbef=MIN(0.5,qsatbef)
569	zcor=1./(1.-retv*qsatbef)
570	qsatbef=qsatbef*zcor
571	qlbef=zqta(ig,l)-qsatbef
572
573	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
574	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
575	zcor=1./(1.-retv*qsatbef)
576	dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
577	num=-Tbef+ztla(ig,l)zpspsk(ig,l)+RLvCpqlbef
578	denom=1.+RLvCp*dqsat_dT
579	DT=num/denom
580	enddo
581	zqla(ig,l) = max(0.,qlbef)
582	endif
583	!
584	if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
585	! on ecrit de maniere conservative (sat ou non)
586	! T = Tl +Lv/Cp ql
587	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
588	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
589	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
590	zha(ig,l) = ztva(ig,l)
591	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
592	& -zqla(ig,l))-zqla(ig,l))
593
594	!on ecrit zqsat
595	zqsatth(ig,l)=qsatbef
596	!calcul de vitesse
597	zw2(ig,l+1)=zw2(ig,l)* &
598	& ((f_star(ig,l))**2) &
599	& /(f_star(ig,l+1)+detr_star(ig,l))**2+ &
600	& 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
601	& *(zlev(ig,l+1)-zlev(ig,l))
602
603	endif
604	endif
605	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
606	!
607	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
608
609	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
610	! stop'On tombe sur le cas particulier de thermcell_dry'
611	print*,'On tombe sur le cas particulier de thermcell_plume'
612	zw2(ig,l+1)=0.
613	linter(ig)=l+1
614	endif
615
616
617	if (zw2(ig,l+1).lt.0.) then
618	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
619	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
620	zw2(ig,l+1)=0.
621	endif
622
623	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
624
625	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
626	! lmix est le niveau de la couche ou w (wa_moy) est maximum
627	lmix(ig)=l+1
628	wmaxa(ig)=wa_moy(ig,l+1)
629	endif
630	enddo
631	enddo
632
633	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
634	#ifdef troisD
635	call wrgradsfi(1,klev,zqla_est,'ql_es ','ql_es ')
636	call wrgradsfi(1,klev,entr_star1,'e_st1 ','e_st1 ')
637	call wrgradsfi(1,klev,entr_star2,'e_st2 ','e_st2 ')
638	call wrgradsfi(1,klev,detr_stara,'d_sta ','d_sta ')
639	call wrgradsfi(1,klev,detr_starb,'d_stb ','d_stb ')
640	call wrgradsfi(1,klev,detr_starc,'d_stc ','d_stc ')
641	call wrgradsfi(1,klev,zqla0 ,'zqla0 ','zqla0 ')
642	call wrgradsfi(1,klev,detr_star0,'d_st0 ','d_st0 ')
643	call wrgradsfi(1,klev,detr_star1,'d_st1 ','d_st1 ')
644	call wrgradsfi(1,klev,detr_star2,'d_st2 ','d_st2 ')
645	#endif
646
647	! print*,'thermcell_plume OK'
648
649	return
650	end

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

Download in other formats:

Original Format