Context Navigation

thermcell_plume.F90 @ 1100

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

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

Download in other formats:

Original Format