Context Navigation

thermcell_plume.F90 @ 5434

Last change on this file since 5434 was 2000, checked in by fhourdin, 11 years ago
Introduction d'un parametre pour controler l'altitude a laquelle on va chercher l'air pour l'entrainement pour obtenir des stratocumulus dz = fact_thermals_ed_dz * z (these Arnaud Jam)
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 Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.8 KB

Line
1	!
2	! $Id: thermcell_plume.F90 2000 2014-03-26 18:05:47Z fhourdin $
3	!
4	SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
5	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
6	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
7	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
8	& ,lev_out,lunout1,igout)
9	!--------------------------------------------------------------------------
10	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
11	!--------------------------------------------------------------------------
12
13	IMPLICIT NONE
14
15	#include "YOMCST.h"
16	#include "YOETHF.h"
17	#include "FCTTRE.h"
18	#include "iniprint.h"
19	#include "thermcell.h"
20
21	INTEGER itap
22	INTEGER lunout1,igout
23	INTEGER ngrid,klev
24	REAL ptimestep
25	REAL ztv(ngrid,klev)
26	REAL zthl(ngrid,klev)
27	REAL po(ngrid,klev)
28	REAL zl(ngrid,klev)
29	REAL rhobarz(ngrid,klev)
30	REAL zlev(ngrid,klev+1)
31	REAL pplev(ngrid,klev+1)
32	REAL pphi(ngrid,klev)
33	REAL zpspsk(ngrid,klev)
34	REAL alim_star(ngrid,klev)
35	REAL f0(ngrid)
36	INTEGER lalim(ngrid)
37	integer lev_out ! niveau pour les print
38	integer nbpb
39
40	real alim_star_tot(ngrid)
41
42	REAL ztva(ngrid,klev)
43	REAL ztla(ngrid,klev)
44	REAL zqla(ngrid,klev)
45	REAL zqta(ngrid,klev)
46	REAL zha(ngrid,klev)
47
48	REAL detr_star(ngrid,klev)
49	REAL coefc
50	REAL entr_star(ngrid,klev)
51	REAL detr(ngrid,klev)
52	REAL entr(ngrid,klev)
53
54	REAL csc(ngrid,klev)
55
56	REAL zw2(ngrid,klev+1)
57	REAL w_est(ngrid,klev+1)
58	REAL f_star(ngrid,klev+1)
59	REAL wa_moy(ngrid,klev+1)
60
61	REAL ztva_est(ngrid,klev)
62	REAL ztv_est(ngrid,klev)
63	REAL zqla_est(ngrid,klev)
64	REAL zqsatth(ngrid,klev)
65	REAL zta_est(ngrid,klev)
66	REAL ztemp(ngrid),zqsat(ngrid)
67	REAL zdw2,zdw2bis
68	REAL zw2modif
69	REAL zw2fact,zw2factbis
70	REAL zeps(ngrid,klev)
71
72	REAL linter(ngrid)
73	INTEGER lmix(ngrid)
74	INTEGER lmix_bis(ngrid)
75	REAL wmaxa(ngrid)
76
77	INTEGER ig,l,k,lt,it
78
79	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
80	real zbuoyjam(ngrid,klev)
81	real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
82	real zcor,zdelta,zcvm5,qlbef
83	real betalpha,zbetalpha
84	real eps, afact
85	REAL REPS,RLvCp,DDT0
86	PARAMETER (DDT0=.01)
87	logical Zsat
88	LOGICAL active(ngrid),activetmp(ngrid)
89	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
90	REAL c2(ngrid,klev)
91	Zsat=.false.
92	! Initialisation
93
94	! print*,'THERMCELL PLUME OK'
95	RLvCp = RLVTT/RCPD
96	fact_epsilon=0.002
97	betalpha=0.9
98	afact=2./3.
99
100	zbetalpha=betalpha/(1.+betalpha)
101
102
103	! Initialisations des variables r?elles
104	if (1==1) then
105	ztva(:,:)=ztv(:,:)
106	ztva_est(:,:)=ztva(:,:)
107	ztv_est(:,:)=ztv(:,:)
108	ztla(:,:)=zthl(:,:)
109	zqta(:,:)=po(:,:)
110	zqla(:,:)=0.
111	zha(:,:) = ztva(:,:)
112	else
113	ztva(:,:)=0.
114	ztv_est(:,:)=0.
115	ztva_est(:,:)=0.
116	ztla(:,:)=0.
117	zqta(:,:)=0.
118	zha(:,:) =0.
119	endif
120
121	zqla_est(:,:)=0.
122	zqsatth(:,:)=0.
123	zqla(:,:)=0.
124	detr_star(:,:)=0.
125	entr_star(:,:)=0.
126	alim_star(:,:)=0.
127	alim_star_tot(:)=0.
128	csc(:,:)=0.
129	detr(:,:)=0.
130	entr(:,:)=0.
131	zw2(:,:)=0.
132	zbuoy(:,:)=0.
133	zbuoyjam(:,:)=0.
134	gamma(:,:)=0.
135	zeps(:,:)=0.
136	w_est(:,:)=0.
137	f_star(:,:)=0.
138	wa_moy(:,:)=0.
139	linter(:)=1.
140	! linter(:)=1.
141	! Initialisation des variables entieres
142	lmix(:)=1
143	lmix_bis(:)=2
144	wmaxa(:)=0.
145	lalim(:)=1
146
147
148	!-------------------------------------------------------------------------
149	! On ne considere comme actif que les colonnes dont les deux premieres
150	! couches sont instables.
151	!-------------------------------------------------------------------------
152	active(:)=ztv(:,1)>ztv(:,2)
153
154	!-------------------------------------------------------------------------
155	! Definition de l'alimentation a l'origine dans thermcell_init
156	!-------------------------------------------------------------------------
157	do l=1,klev-1
158	do ig=1,ngrid
159	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
160	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
161	& *sqrt(zlev(ig,l+1))
162	lalim(ig)=l+1
163	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
164	! print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)
165	endif
166	enddo
167	enddo
168	do l=1,klev
169	do ig=1,ngrid
170	if (alim_star_tot(ig) > 1.e-10 ) then
171	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
172	endif
173	enddo
174	enddo
175	alim_star_tot(:)=1.
176
177
178
179
180	!------------------------------------------------------------------------------
181	! Calcul dans la premiere couche
182	! On decide dans cette version que le thermique n'est actif que si la premiere
183	! couche est instable.
184	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
185	! dans une couche l>1
186	!------------------------------------------------------------------------------
187	do ig=1,ngrid
188	! Le panache va prendre au debut les caracteristiques de l'air contenu
189	! dans cette couche.
190	if (active(ig)) then
191	ztla(ig,1)=zthl(ig,1)
192	zqta(ig,1)=po(ig,1)
193	zqla(ig,1)=zl(ig,1)
194	!cr: attention, prise en compte de f*(1)=1
195	f_star(ig,2)=alim_star(ig,1)
196	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
197	& *(zlev(ig,2)-zlev(ig,1)) &
198	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
199	w_est(ig,2)=zw2(ig,2)
200	endif
201	enddo
202	!
203
204	!==============================================================================
205	!boucle de calcul de la vitesse verticale dans le thermique
206	!==============================================================================
207	do l=2,klev-1
208	!==============================================================================
209
210
211	! On decide si le thermique est encore actif ou non
212	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
213	do ig=1,ngrid
214	active(ig)=active(ig) &
215	& .and. zw2(ig,l)>1.e-10 &
216	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
217	enddo
218
219
220
221	!---------------------------------------------------------------------------
222	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
223	! sans tenir compte du detrainement et de l'entrainement dans cette
224	! couche
225	! C'est a dire qu'on suppose
226	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
227	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
228	! avant) a l'alimentation pour avoir un calcul plus propre
229	!---------------------------------------------------------------------------
230
231	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
232	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
233	do ig=1,ngrid
234	! print*,'active',active(ig),ig,l
235	if(active(ig)) then
236	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
237	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
238	zta_est(ig,l)=ztva_est(ig,l)
239	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
240	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
241	& -zqla_est(ig,l))-zqla_est(ig,l))
242
243
244	!------------------------------------------------
245	!AJAM:nouveau calcul de w?
246	!------------------------------------------------
247	zdz=zlev(ig,l+1)-zlev(ig,l)
248	zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
249	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
250
251	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
252	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
253	zdw2=afact*zbuoy(ig,l)/fact_epsilon
254	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
255	! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
256	! w_est(ig,l+1)=(zdz/zdzbis)Max(0.0001,exp(-zw2fact) &
257	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
258	! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
259	! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))zdw2+w_est(ig,l)exp(-zw2fact))
260	w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
261	& (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
262	& (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
263	if (w_est(ig,l+1).lt.0.) then
264	! w_est(ig,l+1)=zw2(ig,l)
265	w_est(ig,l+1)=0.0001
266	endif
267	endif
268	enddo
269
270
271	!-------------------------------------------------
272	!calcul des taux d'entrainement et de detrainement
273	!-------------------------------------------------
274
275	do ig=1,ngrid
276	if (active(ig)) then
277
278	! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
279	zw2m=w_est(ig,l+1)
280	! zw2m=zw2(ig,l)
281	zdz=zlev(ig,l+1)-zlev(ig,l)
282	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
283	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
284	zbuoybis=zbuoy(ig,l)
285	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
286	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
287
288
289	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
290	! & afact*zbuoybis/zw2m - fact_epsilon )
291
292	! entr_star(ig,l)=MAX(0.,f_star(ig,l)zdzzbetalpha* &
293	! & afact*zbuoybis/zw2m - fact_epsilon )
294
295	!Modif AJAM
296
297	lmel=fact_thermals_ed_dz*zlev(ig,l)
298	! lmel=2.5*(zlev(ig,l)-zlev(ig,l-1))
299	lt=l+1
300	do it=1,klev-(l+1)
301	zdz2=zlev(ig,lt)-zlev(ig,l)
302	if (zdz2.gt.lmel) then
303	zdz3=zlev(ig,lt)-zlev(ig,lt-1)
304	! ztv_est(ig,l)=(lmel/zdz2)*(ztv(ig,lt)-ztv(ig,l))+ztv(ig,l)
305	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)
306
307	zbuoyjam(ig,l)=1.RG(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
308	& ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
309	& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
310
311	! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
312	! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
313	! & po(ig,lt-1))/po(ig,lt-1))
314
315	else
316	lt=lt+1
317	endif
318	enddo
319
320	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
321
322	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
323	! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
324
325	! entrbis=entr_star(ig,l)
326
327
328	detr_star(ig,l)=f_star(ig,l)*zdz &
329	& MAX(1.e-4, -afactzbetalpha*zbuoyjam(ig,l)/zw2m &
330	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
331
332	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
333
334	entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
335	& afact*zbuoy(ig,l)/zw2m - fact_epsilon )
336	! & afactzbuoy(ig,l)/zw2m - fact_epsilon+ 0.012(zdqt(ig,l)/zw2m)**0.5)
337
338
339	! entr_star(ig,l)=Max(0.,f_star(ig,l)zdzzbetalpha* &
340	! & afact*zbuoy(ig,l)/zw2m &
341	! & - 1.*fact_epsilon)
342
343
344	! En dessous de lalim, on prend le max de alim_star et entr_star pour
345	! alim_star et 0 sinon
346	if (l.lt.lalim(ig)) then
347	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
348	entr_star(ig,l)=0.
349	endif
350	! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
351	! alim_star(ig,l)=entrbis
352	! endif
353
354	!print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
355	! Calcul du flux montant normalise
356	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
357	& -detr_star(ig,l)
358
359	endif
360	enddo
361
362
363	!----------------------------------------------------------------------------
364	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
365	!---------------------------------------------------------------------------
366	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
367	do ig=1,ngrid
368	if (activetmp(ig)) then
369	Zsat=.false.
370	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
371	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
372	& /(f_star(ig,l+1)+detr_star(ig,l))
373	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
374	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
375	& /(f_star(ig,l+1)+detr_star(ig,l))
376
377	endif
378	enddo
379
380	ztemp(:)=zpspsk(:,l)*ztla(:,l)
381	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
382	do ig=1,ngrid
383	if (activetmp(ig)) then
384	! on ecrit de maniere conservative (sat ou non)
385	! T = Tl +Lv/Cp ql
386	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
387	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
388	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
389	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
390	zha(ig,l) = ztva(ig,l)
391	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
392	& -zqla(ig,l))-zqla(ig,l))
393	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
394	zdz=zlev(ig,l+1)-zlev(ig,l)
395	zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
396	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
397
398	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
399	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
400	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
401	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
402	! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
403	zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
404	& (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
405	& (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
406	endif
407	enddo
408
409	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
410	!
411	!---------------------------------------------------------------------------
412	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
413	!---------------------------------------------------------------------------
414
415	nbpb=0
416	do ig=1,ngrid
417	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
418	! stop'On tombe sur le cas particulier de thermcell_dry'
419	! print*,'On tombe sur le cas particulier de thermcell_plume'
420	nbpb=nbpb+1
421	zw2(ig,l+1)=0.
422	linter(ig)=l+1
423	endif
424
425	if (zw2(ig,l+1).lt.0.) then
426	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
427	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
428	zw2(ig,l+1)=0.
429	!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
430	elseif (f_star(ig,l+1).lt.0.) then
431	linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
432	& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
433	zw2(ig,l+1)=0.
434	!fin CR:04/05/12
435	endif
436
437	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
438
439	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
440	! lmix est le niveau de la couche ou w (wa_moy) est maximum
441	!on rajoute le calcul de lmix_bis
442	if (zqla(ig,l).lt.1.e-10) then
443	lmix_bis(ig)=l+1
444	endif
445	lmix(ig)=l+1
446	wmaxa(ig)=wa_moy(ig,l+1)
447	endif
448	enddo
449
450	if (nbpb>0) then
451	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
452	endif
453
454	!=========================================================================
455	! FIN DE LA BOUCLE VERTICALE
456	enddo
457	!=========================================================================
458
459	!on recalcule alim_star_tot
460	do ig=1,ngrid
461	alim_star_tot(ig)=0.
462	enddo
463	do ig=1,ngrid
464	do l=1,lalim(ig)-1
465	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
466	enddo
467	enddo
468
469
470	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
471
472	return
473	end
474
475
476
477
478	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
479	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
480	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
481	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
482	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
483	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
484	SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
485	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
486	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
487	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
488	& ,lev_out,lunout1,igout)
489
490	!--------------------------------------------------------------------------
491	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
492	! Version conforme a l'article de Rio et al. 2010.
493	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
494	!--------------------------------------------------------------------------
495
496	IMPLICIT NONE
497
498	#include "YOMCST.h"
499	#include "YOETHF.h"
500	#include "FCTTRE.h"
501	#include "iniprint.h"
502	#include "thermcell.h"
503
504	INTEGER itap
505	INTEGER lunout1,igout
506	INTEGER ngrid,klev
507	REAL ptimestep
508	REAL ztv(ngrid,klev)
509	REAL zthl(ngrid,klev)
510	REAL po(ngrid,klev)
511	REAL zl(ngrid,klev)
512	REAL rhobarz(ngrid,klev)
513	REAL zlev(ngrid,klev+1)
514	REAL pplev(ngrid,klev+1)
515	REAL pphi(ngrid,klev)
516	REAL zpspsk(ngrid,klev)
517	REAL alim_star(ngrid,klev)
518	REAL f0(ngrid)
519	INTEGER lalim(ngrid)
520	integer lev_out ! niveau pour les print
521	integer nbpb
522
523	real alim_star_tot(ngrid)
524
525	REAL ztva(ngrid,klev)
526	REAL ztla(ngrid,klev)
527	REAL zqla(ngrid,klev)
528	REAL zqta(ngrid,klev)
529	REAL zha(ngrid,klev)
530
531	REAL detr_star(ngrid,klev)
532	REAL coefc
533	REAL entr_star(ngrid,klev)
534	REAL detr(ngrid,klev)
535	REAL entr(ngrid,klev)
536
537	REAL csc(ngrid,klev)
538
539	REAL zw2(ngrid,klev+1)
540	REAL w_est(ngrid,klev+1)
541	REAL f_star(ngrid,klev+1)
542	REAL wa_moy(ngrid,klev+1)
543
544	REAL ztva_est(ngrid,klev)
545	REAL zqla_est(ngrid,klev)
546	REAL zqsatth(ngrid,klev)
547	REAL zta_est(ngrid,klev)
548	REAL zbuoyjam(ngrid,klev)
549	REAL ztemp(ngrid),zqsat(ngrid)
550	REAL zdw2
551	REAL zw2modif
552	REAL zw2fact
553	REAL zeps(ngrid,klev)
554
555	REAL linter(ngrid)
556	INTEGER lmix(ngrid)
557	INTEGER lmix_bis(ngrid)
558	REAL wmaxa(ngrid)
559
560	INTEGER ig,l,k
561
562	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
563	real zbuoybis
564	real zcor,zdelta,zcvm5,qlbef,zdz2
565	real betalpha,zbetalpha
566	real eps, afact
567	REAL REPS,RLvCp,DDT0
568	PARAMETER (DDT0=.01)
569	logical Zsat
570	LOGICAL active(ngrid),activetmp(ngrid)
571	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
572	REAL c2(ngrid,klev)
573	Zsat=.false.
574	! Initialisation
575
576	RLvCp = RLVTT/RCPD
577	fact_epsilon=0.002
578	betalpha=0.9
579	afact=2./3.
580
581	zbetalpha=betalpha/(1.+betalpha)
582
583
584	! Initialisations des variables reeles
585	if (1==1) then
586	ztva(:,:)=ztv(:,:)
587	ztva_est(:,:)=ztva(:,:)
588	ztla(:,:)=zthl(:,:)
589	zqta(:,:)=po(:,:)
590	zha(:,:) = ztva(:,:)
591	else
592	ztva(:,:)=0.
593	ztva_est(:,:)=0.
594	ztla(:,:)=0.
595	zqta(:,:)=0.
596	zha(:,:) =0.
597	endif
598
599	zqla_est(:,:)=0.
600	zqsatth(:,:)=0.
601	zqla(:,:)=0.
602	detr_star(:,:)=0.
603	entr_star(:,:)=0.
604	alim_star(:,:)=0.
605	alim_star_tot(:)=0.
606	csc(:,:)=0.
607	detr(:,:)=0.
608	entr(:,:)=0.
609	zw2(:,:)=0.
610	zbuoy(:,:)=0.
611	zbuoyjam(:,:)=0.
612	gamma(:,:)=0.
613	zeps(:,:)=0.
614	w_est(:,:)=0.
615	f_star(:,:)=0.
616	wa_moy(:,:)=0.
617	linter(:)=1.
618	! linter(:)=1.
619	! Initialisation des variables entieres
620	lmix(:)=1
621	lmix_bis(:)=2
622	wmaxa(:)=0.
623	lalim(:)=1
624
625
626	!-------------------------------------------------------------------------
627	! On ne considere comme actif que les colonnes dont les deux premieres
628	! couches sont instables.
629	!-------------------------------------------------------------------------
630	active(:)=ztv(:,1)>ztv(:,2)
631
632	!-------------------------------------------------------------------------
633	! Definition de l'alimentation a l'origine dans thermcell_init
634	!-------------------------------------------------------------------------
635	do l=1,klev-1
636	do ig=1,ngrid
637	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
638	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
639	& *sqrt(zlev(ig,l+1))
640	lalim(ig)=l+1
641	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
642	endif
643	enddo
644	enddo
645	do l=1,klev
646	do ig=1,ngrid
647	if (alim_star_tot(ig) > 1.e-10 ) then
648	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
649	endif
650	enddo
651	enddo
652	alim_star_tot(:)=1.
653
654
655
656	!------------------------------------------------------------------------------
657	! Calcul dans la premiere couche
658	! On decide dans cette version que le thermique n'est actif que si la premiere
659	! couche est instable.
660	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
661	! dans une couche l>1
662	!------------------------------------------------------------------------------
663	do ig=1,ngrid
664	! Le panache va prendre au debut les caracteristiques de l'air contenu
665	! dans cette couche.
666	if (active(ig)) then
667	ztla(ig,1)=zthl(ig,1)
668	zqta(ig,1)=po(ig,1)
669	zqla(ig,1)=zl(ig,1)
670	!cr: attention, prise en compte de f*(1)=1
671	f_star(ig,2)=alim_star(ig,1)
672	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
673	& *(zlev(ig,2)-zlev(ig,1)) &
674	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
675	w_est(ig,2)=zw2(ig,2)
676	endif
677	enddo
678	!
679
680	!==============================================================================
681	!boucle de calcul de la vitesse verticale dans le thermique
682	!==============================================================================
683	do l=2,klev-1
684	!==============================================================================
685
686
687	! On decide si le thermique est encore actif ou non
688	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
689	do ig=1,ngrid
690	active(ig)=active(ig) &
691	& .and. zw2(ig,l)>1.e-10 &
692	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
693	enddo
694
695
696
697	!---------------------------------------------------------------------------
698	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
699	! sans tenir compte du detrainement et de l'entrainement dans cette
700	! couche
701	! C'est a dire qu'on suppose
702	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
703	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
704	! avant) a l'alimentation pour avoir un calcul plus propre
705	!---------------------------------------------------------------------------
706
707	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
708	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
709
710	do ig=1,ngrid
711	! print*,'active',active(ig),ig,l
712	if(active(ig)) then
713	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
714	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
715	zta_est(ig,l)=ztva_est(ig,l)
716	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
717	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
718	& -zqla_est(ig,l))-zqla_est(ig,l))
719
720	!------------------------------------------------
721	!AJAM:nouveau calcul de w?
722	!------------------------------------------------
723	zdz=zlev(ig,l+1)-zlev(ig,l)
724	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
725
726	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
727	zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
728	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
729
730
731	if (w_est(ig,l+1).lt.0.) then
732	w_est(ig,l+1)=zw2(ig,l)
733	endif
734	endif
735	enddo
736
737
738	!-------------------------------------------------
739	!calcul des taux d'entrainement et de detrainement
740	!-------------------------------------------------
741
742	do ig=1,ngrid
743	if (active(ig)) then
744
745	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
746	zw2m=w_est(ig,l+1)
747	zdz=zlev(ig,l+1)-zlev(ig,l)
748	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
749	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
750	zbuoybis=zbuoy(ig,l)
751	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
752	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
753
754
755	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
756	& afact*zbuoybis/zw2m - fact_epsilon )
757
758
759	detr_star(ig,l)=f_star(ig,l)*zdz &
760	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
761	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
762
763	! En dessous de lalim, on prend le max de alim_star et entr_star pour
764	! alim_star et 0 sinon
765	if (l.lt.lalim(ig)) then
766	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
767	entr_star(ig,l)=0.
768	endif
769
770	! Calcul du flux montant normalise
771	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
772	& -detr_star(ig,l)
773
774	endif
775	enddo
776
777
778	!----------------------------------------------------------------------------
779	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
780	!---------------------------------------------------------------------------
781	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
782	do ig=1,ngrid
783	if (activetmp(ig)) then
784	Zsat=.false.
785	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
786	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
787	& /(f_star(ig,l+1)+detr_star(ig,l))
788	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
789	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
790	& /(f_star(ig,l+1)+detr_star(ig,l))
791
792	endif
793	enddo
794
795	ztemp(:)=zpspsk(:,l)*ztla(:,l)
796	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
797
798	do ig=1,ngrid
799	if (activetmp(ig)) then
800	! on ecrit de maniere conservative (sat ou non)
801	! T = Tl +Lv/Cp ql
802	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
803	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
804	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
805	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
806	zha(ig,l) = ztva(ig,l)
807	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
808	& -zqla(ig,l))-zqla(ig,l))
809	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
810	zdz=zlev(ig,l+1)-zlev(ig,l)
811	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
812
813	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
814	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
815	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
816	endif
817	enddo
818
819	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
820	!
821	!---------------------------------------------------------------------------
822	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
823	!---------------------------------------------------------------------------
824
825	nbpb=0
826	do ig=1,ngrid
827	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
828	! stop'On tombe sur le cas particulier de thermcell_dry'
829	! print*,'On tombe sur le cas particulier de thermcell_plume'
830	nbpb=nbpb+1
831	zw2(ig,l+1)=0.
832	linter(ig)=l+1
833	!CR:04/05/12:calcul linter
834	elseif (f_star(ig,l+1).lt.0.) then
835	linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
836	& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
837	zw2(ig,l+1)=0.
838	!fin CR:04/05/12
839	endif
840
841
842	if (zw2(ig,l+1).lt.0.) then
843	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
844	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
845	zw2(ig,l+1)=0.
846	endif
847
848	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
849
850	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
851	! lmix est le niveau de la couche ou w (wa_moy) est maximum
852	!on rajoute le calcul de lmix_bis
853	if (zqla(ig,l).lt.1.e-10) then
854	lmix_bis(ig)=l+1
855	endif
856	lmix(ig)=l+1
857	wmaxa(ig)=wa_moy(ig,l+1)
858	endif
859	enddo
860
861	if (nbpb>0) then
862	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
863	endif
864
865	!=========================================================================
866	! FIN DE LA BOUCLE VERTICALE
867	enddo
868	!=========================================================================
869
870	!on recalcule alim_star_tot
871	do ig=1,ngrid
872	alim_star_tot(ig)=0.
873	enddo
874	do ig=1,ngrid
875	do l=1,lalim(ig)-1
876	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
877	enddo
878	enddo
879
880
881	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
882
883
884
885	return
886	end

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

Download in other formats:

Original Format