Context Navigation

thermcell_plume.F90 @ 1088

Last change on this file since 1088 was 852, checked in by Laurent Fairhead, 17 years ago
Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 13.3 KB

Line
1	SUBROUTINE thermcell_plume(ngrid,klev,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
16	INTEGER ngrid,klev
17	REAL ztv(ngrid,klev)
18	REAL zthl(ngrid,klev)
19	REAL po(ngrid,klev)
20	REAL zl(ngrid,klev)
21	REAL rhobarz(ngrid,klev)
22	REAL zlev(ngrid,klev+1)
23	REAL pplev(ngrid,klev+1)
24	REAL pphi(ngrid,klev)
25	REAL zpspsk(ngrid,klev)
26	REAL alim_star(ngrid,klev)
27	REAL zmax_sec(ngrid)
28	REAL f0(ngrid)
29	REAL l_mix
30	REAL r_aspect
31	INTEGER lalim(ngrid)
32	integer lev_out ! niveau pour les print
33
34	REAL ztva(ngrid,klev)
35	REAL ztla(ngrid,klev)
36	REAL zqla(ngrid,klev)
37	REAL zqta(ngrid,klev)
38	REAL zha(ngrid,klev)
39
40	REAL detr_star(ngrid,klev)
41	REAL entr_star(ngrid,klev)
42	REAL detr(ngrid,klev)
43	REAL entr(ngrid,klev)
44
45	REAL zw2(ngrid,klev+1)
46	REAL w_est(ngrid,klev+1)
47	REAL f_star(ngrid,klev+1)
48	REAL wa_moy(ngrid,klev+1)
49
50	REAL ztva_est(ngrid,klev)
51	REAL zqla_est(ngrid,klev)
52	REAL zqsatth(ngrid,klev)
53
54	REAL linter(ngrid)
55	INTEGER lmix(ngrid)
56	REAL wmaxa(ngrid)
57
58	INTEGER ig,l,k
59
60	real zcor,zdelta,zcvm5,qlbef
61	real Tbef,qsatbef
62	real dqsat_dT,DT,num,denom
63	REAL REPS,RLvCp,DDT0
64	PARAMETER (DDT0=.01)
65	logical Zsat
66
67	Zsat=.false.
68	! Initialisation
69	RLvCp = RLVTT/RCPD
70
71	do l=1,klev
72	do ig=1,ngrid
73	zqla_est(ig,l)=0.
74	ztva_est(ig,l)=ztva(ig,l)
75	zqsatth(ig,l)=0.
76	enddo
77	enddo
78
79	!AM:initialisations du thermique
80	do k=1,klev
81	do ig=1,ngrid
82	ztva(ig,k)=ztv(ig,k)
83	ztla(ig,k)=zthl(ig,k)
84	zqla(ig,k)=0.
85	zqta(ig,k)=po(ig,k)
86	enddo
87	enddo
88	do k=1,klev
89	do ig=1,ngrid
90	detr_star(ig,k)=0.
91	entr_star(ig,k)=0.
92	detr(ig,k)=0.
93	entr(ig,k)=0.
94	enddo
95	enddo
96	if (lev_out.ge.1) print*,'7 OK convect8'
97	do k=1,klev+1
98	do ig=1,ngrid
99	zw2(ig,k)=0.
100	w_est(ig,k)=0.
101	f_star(ig,k)=0.
102	wa_moy(ig,k)=0.
103	enddo
104	enddo
105
106	if (lev_out.ge.1) print*,'8 OK convect8'
107	do ig=1,ngrid
108	linter(ig)=1.
109	lmix(ig)=1
110	wmaxa(ig)=0.
111	enddo
112
113	!-----------------------------------------------------------------------------------
114	!boucle de calcul de la vitesse verticale dans le thermique
115	!-----------------------------------------------------------------------------------
116	do l=1,klev-1
117	do ig=1,ngrid
118	if (ztv(ig,l).gt.ztv(ig,l+1) &
119	& .and.alim_star(ig,l).gt.1.e-10 &
120	& .and.zw2(ig,l).lt.1e-10) then
121	ztla(ig,l)=zthl(ig,l)
122	zqta(ig,l)=po(ig,l)
123	zqla(ig,l)=zl(ig,l)
124	f_star(ig,l+1)=alim_star(ig,l)
125	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) &
126	& *(zlev(ig,l+1)-zlev(ig,l)) &
127	& 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
128	w_est(ig,l+1)=zw2(ig,l+1)
129	!
130	else if ((zw2(ig,l).ge.1e-10).and. &
131	& (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
132	!estimation du detrainement a partir de la geometrie du pas precedent
133	!tests sur la definition du detr
134	!calcul de detr_star et entr_star
135
136
137
138	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
139	! FH le test miraculeux de Catherine ? Le bout du tunel ?
140	! w_est(ig,3)=zw2(ig,2)* &
141	! & ((f_star(ig,2))**2) &
142	! & /(f_star(ig,2)+alim_star(ig,2))**2+ &
143	! & 2.RG(ztva(ig,1)-ztv(ig,2))/ztv(ig,2) &
144	! & *(zlev(ig,3)-zlev(ig,2))
145	! if (l.gt.2) then
146	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
147
148	if (l.gt.2) then
149	w_est(ig,3)=zw2(ig,2)* &
150	& ((f_star(ig,2))**2) &
151	& /(f_star(ig,2)+alim_star(ig,2))**2+ &
152	& 2.RG(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
153	& *(zlev(ig,3)-zlev(ig,2))
154
155
156	!---------------------------------------------------------------------------
157	!calcul de l entrainement et du detrainement lateral
158	!---------------------------------------------------------------------------
159	!
160	!test:estimation de ztva_new_est sans entrainement
161	Tbef=ztla(ig,l-1)*zpspsk(ig,l)
162	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
163	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
164	qsatbef=MIN(0.5,qsatbef)
165	zcor=1./(1.-retv*qsatbef)
166	qsatbef=qsatbef*zcor
167	Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
168	if (Zsat) then
169	qlbef=max(0.,zqta(ig,l-1)-qsatbef)
170	DT = 0.5RLvCpqlbef
171	do while (abs(DT).gt.DDT0)
172	Tbef=Tbef+DT
173	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
174	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
175	qsatbef=MIN(0.5,qsatbef)
176	zcor=1./(1.-retv*qsatbef)
177	qsatbef=qsatbef*zcor
178	qlbef=zqta(ig,l-1)-qsatbef
179
180	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
181	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
182	zcor=1./(1.-retv*qsatbef)
183	dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
184	num=-Tbef+ztla(ig,l-1)zpspsk(ig,l)+RLvCpqlbef
185	denom=1.+RLvCp*dqsat_dT
186	DT=num/denom
187	enddo
188	zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef)
189	endif
190	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
191	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
192	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
193	& -zqla_est(ig,l))-zqla_est(ig,l))
194
195	w_est(ig,l+1)=zw2(ig,l)* &
196	& ((f_star(ig,l))**2) &
197	& /(f_star(ig,l)+alim_star(ig,l))**2+ &
198	& 2.RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
199	& *(zlev(ig,l+1)-zlev(ig,l))
200	if (w_est(ig,l+1).lt.0.) then
201	w_est(ig,l+1)=zw2(ig,l)
202	endif
203	!
204	!calcul du detrainement
205	if ((w_est(ig,l+1).gt.w_est(ig,l)).and. &
206	& (zlev(ig,l+1).lt.zmax_sec(ig)).and. &
207	& (zqla(ig,l-1).lt.1.e-10)) then
208	detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1) &
209	& sqrt(w_est(ig,l+1))sqrt(l_mix*zlev(ig,l+1)) &
210	& -rhobarz(ig,l)sqrt(w_est(ig,l))sqrt(l_mix*zlev(ig,l))) &
211	& /(r_aspect*zmax_sec(ig)))
212	if (lev_out.ge.20) print*,'coucou calcul detr 1'
213	else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and. &
214	& (zqla(ig,l-1).lt.1.e-10)) then
215	detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig)) &
216	& /(rhobarz(ig,lmix(ig))wmaxa(ig)) &
217	& (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1)) &
218	& *((zmax_sec(ig)-zlev(ig,l+1))/ &
219	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2. &
220	& -rhobarz(ig,l)*sqrt(w_est(ig,l)) &
221	& *((zmax_sec(ig)-zlev(ig,l))/ &
222	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
223	if (lev_out.ge.20) print*,'coucou calcul detr 2'
224	else
225	detr_star(ig,l)=0.002f0(ig)f_star(ig,l) &
226	& *(zlev(ig,l+1)-zlev(ig,l))
227	if (lev_out.ge.20) print*,'coucou calcul detr 3'
228
229	endif
230	detr_star(ig,l)=detr_star(ig,l)/f0(ig)
231	!
232	!calcul de entr_star
233	!
234	! Deplacement du calcul de entr_star pour eviter d'avoir aussi
235	! entr_star > fstar.
236	! Redeplacer suite a la transformation du cas detr>f
237	! FH
238	if(l.gt.lalim(ig)) then
239	entr_star(ig,l)=0.4*detr_star(ig,l)
240	else
241
242	! FH :
243	! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
244	! en haut de la couche d'alimentation.
245	! A remettre en questoin a la premiere occasion mais ca peut aider a
246	! ecrire un code robuste.
247	! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
248	! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
249	! d'alimentation, ce qui n'est pas forcement heureux.
250	entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
251	detr_star(ig,l)=entr_star(ig,l)
252	endif
253
254	!
255	if (detr_star(ig,l).gt.f_star(ig,l)) then
256
257	! Ce test est là entre autres parce qu'on passe par des valeurs
258	! delirantes de detr_star.
259	! ca vaut sans doute le coup de verifier pourquoi.
260
261	detr_star(ig,l)=f_star(ig,l)
262	if (l.gt.lalim(ig)+1) then
263	entr_star(ig,l)=0.
264	alim_star(ig,l)=0.
265	! FH ajout pour forcer a stoper le thermique juste sous le sommet
266	! de la couche (voir calcul de finter)
267	zw2(ig,l+1)=-1.e-10
268	linter(ig)=l+1
269	else
270	entr_star(ig,l)=detr_star(ig,l)
271	endif
272	endif
273
274	else
275	detr_star(ig,l)=0.
276	entr_star(ig,l)=0.
277	endif
278
279	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
280	! FH inutile si on conserve comme on l'a fait plus haut entr=detr
281	! dans la couche d'alimentation
282	!pas d entrainement dans la couche alim
283	! if ((l.le.lalim(ig))) then
284	! entr_star(ig,l)=0.
285	! endif
286	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
287	!
288	!prise en compte du detrainement et de l entrainement dans le calcul du flux
289	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
290	& -detr_star(ig,l)
291
292	!test sur le signe de f_star
293	if (f_star(ig,l+1).gt.1.e-10) then
294	!----------------------------------------------------------------------------
295	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
296	!---------------------------------------------------------------------------
297	!
298	Zsat=.false.
299	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
300	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
301	& /(f_star(ig,l+1)+detr_star(ig,l))
302	!
303	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
304	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
305	& /(f_star(ig,l+1)+detr_star(ig,l))
306	!
307	Tbef=ztla(ig,l)*zpspsk(ig,l)
308	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
309	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
310	qsatbef=MIN(0.5,qsatbef)
311	zcor=1./(1.-retv*qsatbef)
312	qsatbef=qsatbef*zcor
313	Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
314	if (Zsat) then
315	qlbef=max(0.,zqta(ig,l)-qsatbef)
316	DT = 0.5RLvCpqlbef
317	do while (abs(DT).gt.DDT0)
318	Tbef=Tbef+DT
319	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
320	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
321	qsatbef=MIN(0.5,qsatbef)
322	zcor=1./(1.-retv*qsatbef)
323	qsatbef=qsatbef*zcor
324	qlbef=zqta(ig,l)-qsatbef
325
326	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
327	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
328	zcor=1./(1.-retv*qsatbef)
329	dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
330	num=-Tbef+ztla(ig,l)zpspsk(ig,l)+RLvCpqlbef
331	denom=1.+RLvCp*dqsat_dT
332	DT=num/denom
333	enddo
334	zqla(ig,l) = max(0.,qlbef)
335	endif
336	!
337	! on ecrit de maniere conservative (sat ou non)
338	! T = Tl +Lv/Cp ql
339	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
340	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
341	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
342	zha(ig,l) = ztva(ig,l)
343	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
344	& -zqla(ig,l))-zqla(ig,l))
345
346	!on ecrit zqsat
347	zqsatth(ig,l)=qsatbef
348	!calcul de vitesse
349	zw2(ig,l+1)=zw2(ig,l)* &
350	& ((f_star(ig,l))**2) &
351	& /(f_star(ig,l+1)+detr_star(ig,l))**2+ &
352	& 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
353	& *(zlev(ig,l+1)-zlev(ig,l))
354
355	endif
356	endif
357	!
358	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
359
360	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
361	! stop'On tombe sur le cas particulier de thermcell_dry'
362	print*,'On tombe sur le cas particulier de thermcell_plume'
363	zw2(ig,l+1)=0.
364	linter(ig)=l+1
365	endif
366
367
368	if (zw2(ig,l+1).lt.0.) then
369	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
370	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
371	zw2(ig,l+1)=0.
372	endif
373
374	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
375
376	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
377	! lmix est le niveau de la couche ou w (wa_moy) est maximum
378	lmix(ig)=l+1
379	wmaxa(ig)=wa_moy(ig,l+1)
380	endif
381	enddo
382	enddo
383
384
385	return
386	end

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

Download in other formats:

Original Format