Context Navigation

thermcell_plume.F90 @ 1968

Last change on this file since 1968 was 1968, checked in by fhourdin, 10 years ago

Nouvelle version de thermcell_plume, permettant d'avoir des stratocumulus

avec 40 couches dans le cas fire.

New version of thermcell_plume, giving stratocumulus in the fire case

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: 29.2 KB

Line
1	!
2	! $Id: thermcell_plume.F90 1968 2014-02-11 15:01:12Z 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	!--------------------------------------------------------------------------
11	! thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
12	! Last modified : Arnaud Jam 2014/02/11
13	! Better representation of stratocumulus
14	!--------------------------------------------------------------------------
15
16	IMPLICIT NONE
17
18	#include "YOMCST.h"
19	#include "YOETHF.h"
20	#include "FCTTRE.h"
21	#include "iniprint.h"
22	#include "thermcell.h"
23
24	INTEGER itap
25	INTEGER lunout1,igout
26	INTEGER ngrid,klev
27	REAL ptimestep
28	REAL ztv(ngrid,klev)
29	REAL zthl(ngrid,klev)
30	REAL po(ngrid,klev)
31	REAL zl(ngrid,klev)
32	REAL rhobarz(ngrid,klev)
33	REAL zlev(ngrid,klev+1)
34	REAL pplev(ngrid,klev+1)
35	REAL pphi(ngrid,klev)
36	REAL zpspsk(ngrid,klev)
37	REAL alim_star(ngrid,klev)
38	REAL f0(ngrid)
39	INTEGER lalim(ngrid)
40	integer lev_out ! niveau pour les print
41	integer nbpb
42
43	real alim_star_tot(ngrid)
44
45	REAL ztva(ngrid,klev)
46	REAL ztla(ngrid,klev)
47	REAL zqla(ngrid,klev)
48	REAL zqta(ngrid,klev)
49	REAL zha(ngrid,klev)
50
51	REAL detr_star(ngrid,klev)
52	REAL coefc
53	REAL entr_star(ngrid,klev)
54	REAL detr(ngrid,klev)
55	REAL entr(ngrid,klev)
56
57	REAL csc(ngrid,klev)
58
59	REAL zw2(ngrid,klev+1)
60	REAL w_est(ngrid,klev+1)
61	REAL f_star(ngrid,klev+1)
62	REAL wa_moy(ngrid,klev+1)
63
64	REAL ztva_est(ngrid,klev)
65	REAL ztv_est(ngrid,klev)
66	REAL zqla_est(ngrid,klev)
67	REAL zqsatth(ngrid,klev)
68	REAL zta_est(ngrid,klev)
69	REAL ztemp(ngrid),zqsat(ngrid)
70	REAL zdw2,zdw2bis
71	REAL zw2modif
72	REAL zw2fact,zw2factbis
73	REAL zeps(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,lt,it
81
82	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
83	real zbuoyjam(ngrid,klev)
84	real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
85	real zcor,zdelta,zcvm5,qlbef
86	real betalpha,zbetalpha
87	real eps, afact
88	REAL REPS,RLvCp,DDT0
89	PARAMETER (DDT0=.01)
90	logical Zsat
91	LOGICAL active(ngrid),activetmp(ngrid)
92	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
93	REAL c2(ngrid,klev)
94	Zsat=.false.
95	! Initialisation
96
97	RLvCp = RLVTT/RCPD
98	fact_epsilon=0.002
99	betalpha=0.9
100	afact=2./3.
101
102	zbetalpha=betalpha/(1.+betalpha)
103
104
105	! Initialisations des variables r?elles
106	if (1==1) then
107	ztva(:,:)=ztv(:,:)
108	ztva_est(:,:)=ztva(:,:)
109	ztv_est(:,:)=ztv(:,:)
110	ztla(:,:)=zthl(:,:)
111	zqta(:,:)=po(:,:)
112	zqla(:,:)=0.
113	zha(:,:) = ztva(:,:)
114	else
115	ztva(:,:)=0.
116	ztv_est(:,:)=0.
117	ztva_est(:,:)=0.
118	ztla(:,:)=0.
119	zqta(:,:)=0.
120	zha(:,:) =0.
121	endif
122
123	zqla_est(:,:)=0.
124	zqsatth(:,:)=0.
125	zqla(:,:)=0.
126	detr_star(:,:)=0.
127	entr_star(:,:)=0.
128	alim_star(:,:)=0.
129	alim_star_tot(:)=0.
130	csc(:,:)=0.
131	detr(:,:)=0.
132	entr(:,:)=0.
133	zw2(:,:)=0.
134	zbuoy(:,:)=0.
135	zbuoyjam(:,:)=0.
136	gamma(:,:)=0.
137	zeps(:,:)=0.
138	w_est(:,:)=0.
139	f_star(:,:)=0.
140	wa_moy(:,:)=0.
141	linter(:)=1.
142	! linter(:)=1.
143	! Initialisation des variables entieres
144	lmix(:)=1
145	lmix_bis(:)=2
146	wmaxa(:)=0.
147	lalim(:)=1
148
149
150	!-------------------------------------------------------------------------
151	! On ne considere comme actif que les colonnes dont les deux premieres
152	! couches sont instables.
153	!-------------------------------------------------------------------------
154	active(:)=ztv(:,1)>ztv(:,2)
155
156	!-------------------------------------------------------------------------
157	! Definition de l'alimentation a l'origine dans thermcell_init
158	!-------------------------------------------------------------------------
159	do l=1,klev-1
160	do ig=1,ngrid
161	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
162	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
163	& *sqrt(zlev(ig,l+1))
164	lalim(ig)=l+1
165	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
166	print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)
167	endif
168	enddo
169	enddo
170	do l=1,klev
171	do ig=1,ngrid
172	if (alim_star_tot(ig) > 1.e-10 ) then
173	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
174	endif
175	enddo
176	enddo
177	alim_star_tot(:)=1.
178
179
180
181
182	!------------------------------------------------------------------------------
183	! Calcul dans la premiere couche
184	! On decide dans cette version que le thermique n'est actif que si la premiere
185	! couche est instable.
186	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
187	! dans une couche l>1
188	!------------------------------------------------------------------------------
189	do ig=1,ngrid
190	! Le panache va prendre au debut les caracteristiques de l'air contenu
191	! dans cette couche.
192	if (active(ig)) then
193	ztla(ig,1)=zthl(ig,1)
194	zqta(ig,1)=po(ig,1)
195	zqla(ig,1)=zl(ig,1)
196	!cr: attention, prise en compte de f*(1)=1
197	f_star(ig,2)=alim_star(ig,1)
198	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
199	& *(zlev(ig,2)-zlev(ig,1)) &
200	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
201	w_est(ig,2)=zw2(ig,2)
202	endif
203	enddo
204	!
205
206	!==============================================================================
207	!boucle de calcul de la vitesse verticale dans le thermique
208	!==============================================================================
209	do l=2,klev-1
210	!==============================================================================
211
212
213	! On decide si le thermique est encore actif ou non
214	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
215	do ig=1,ngrid
216	active(ig)=active(ig) &
217	& .and. zw2(ig,l)>1.e-10 &
218	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
219	enddo
220
221
222
223	!---------------------------------------------------------------------------
224	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
225	! sans tenir compte du detrainement et de l'entrainement dans cette
226	! couche
227	! C'est a dire qu'on suppose
228	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
229	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
230	! avant) a l'alimentation pour avoir un calcul plus propre
231	!---------------------------------------------------------------------------
232
233	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
234	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
235	do ig=1,ngrid
236	print*,'active',active(ig),ig,l
237	if(active(ig)) then
238	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
239	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
240	zta_est(ig,l)=ztva_est(ig,l)
241	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
242	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
243	& -zqla_est(ig,l))-zqla_est(ig,l))
244
245
246	!------------------------------------------------
247	!AJAM:nouveau calcul de w?
248	!------------------------------------------------
249	zdz=zlev(ig,l+1)-zlev(ig,l)
250	zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
251	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
252
253	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
254	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
255	zdw2=afact*zbuoy(ig,l)/fact_epsilon
256	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
257	! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
258	! w_est(ig,l+1)=(zdz/zdzbis)Max(0.0001,exp(-zw2fact) &
259	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
260	! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
261	! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))zdw2+w_est(ig,l)exp(-zw2fact))
262	w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
263	& (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
264	& (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
265	if (w_est(ig,l+1).lt.0.) then
266	! w_est(ig,l+1)=zw2(ig,l)
267	w_est(ig,l+1)=0.0001
268	endif
269	endif
270	enddo
271
272
273	!-------------------------------------------------
274	!calcul des taux d'entrainement et de detrainement
275	!-------------------------------------------------
276
277	do ig=1,ngrid
278	if (active(ig)) then
279
280	! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
281	zw2m=w_est(ig,l+1)
282	! zw2m=zw2(ig,l)
283	zdz=zlev(ig,l+1)-zlev(ig,l)
284	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
285	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
286	zbuoybis=zbuoy(ig,l)
287	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
288	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
289
290
291	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
292	! & afact*zbuoybis/zw2m - fact_epsilon )
293
294	! entr_star(ig,l)=MAX(0.,f_star(ig,l)zdzzbetalpha* &
295	! & afact*zbuoybis/zw2m - fact_epsilon )
296
297	!Modif AJAM
298
299	lmel=0.1*zlev(ig,l)
300	! lmel=2.5*(zlev(ig,l)-zlev(ig,l-1))
301	lt=l+1
302	do it=1,klev-(l+1)
303	zdz2=zlev(ig,lt)-zlev(ig,l)
304	if (zdz2.gt.lmel) then
305	zdz3=zlev(ig,lt)-zlev(ig,lt-1)
306	! ztv_est(ig,l)=(lmel/zdz2)*(ztv(ig,lt)-ztv(ig,l))+ztv(ig,l)
307	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)
308
309	zbuoyjam(ig,l)=1.RG(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
310	& ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
311	& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
312
313	! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
314	! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
315	! & po(ig,lt-1))/po(ig,lt-1))
316
317	else
318	lt=lt+1
319	endif
320	enddo
321
322	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
323
324	entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
325	& afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
326
327	entrbis=entr_star(ig,l)
328
329
330	detr_star(ig,l)=f_star(ig,l)*zdz &
331	& MAX(1.e-4, -afactzbetalpha*zbuoyjam(ig,l)/zw2m &
332	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
333
334
335	! zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
336	!
337	! entr_star(ig,l)=Max(0.,f_star(ig,l)zdzzbetalpha* &
338	! & afact*zbuoy(ig,l)/zw2m &
339	! & - 1.*fact_epsilon)
340
341
342	! En dessous de lalim, on prend le max de alim_star et entr_star pour
343	! alim_star et 0 sinon
344	if (l.lt.lalim(ig)) then
345	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
346	entr_star(ig,l)=0.
347	endif
348	! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
349	! alim_star(ig,l)=entrbis
350	! endif
351
352	print*,'alim0',l,lalim(ig),alim_star(ig,l),entrbis,f_star(ig,l)
353	! Calcul du flux montant normalise
354	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
355	& -detr_star(ig,l)
356
357	endif
358	enddo
359
360
361	!----------------------------------------------------------------------------
362	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
363	!---------------------------------------------------------------------------
364	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
365	do ig=1,ngrid
366	if (activetmp(ig)) then
367	Zsat=.false.
368	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
369	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
370	& /(f_star(ig,l+1)+detr_star(ig,l))
371	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
372	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
373	& /(f_star(ig,l+1)+detr_star(ig,l))
374
375	endif
376	enddo
377
378	ztemp(:)=zpspsk(:,l)*ztla(:,l)
379	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
380	do ig=1,ngrid
381	if (activetmp(ig)) then
382	! on ecrit de maniere conservative (sat ou non)
383	! T = Tl +Lv/Cp ql
384	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
385	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
386	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
387	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
388	zha(ig,l) = ztva(ig,l)
389	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
390	& -zqla(ig,l))-zqla(ig,l))
391	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
392	zdz=zlev(ig,l+1)-zlev(ig,l)
393	zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
394	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
395
396	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
397	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
398	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
399	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
400	! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
401	zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
402	& (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
403	& (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
404	endif
405	enddo
406
407	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
408	!
409	!---------------------------------------------------------------------------
410	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
411	!---------------------------------------------------------------------------
412
413	nbpb=0
414	do ig=1,ngrid
415	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
416	! stop'On tombe sur le cas particulier de thermcell_dry'
417	! print*,'On tombe sur le cas particulier de thermcell_plume'
418	nbpb=nbpb+1
419	zw2(ig,l+1)=0.
420	linter(ig)=l+1
421	endif
422
423	if (zw2(ig,l+1).lt.0.) then
424	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
425	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
426	zw2(ig,l+1)=0.
427	endif
428
429	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
430
431	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
432	! lmix est le niveau de la couche ou w (wa_moy) est maximum
433	!on rajoute le calcul de lmix_bis
434	if (zqla(ig,l).lt.1.e-10) then
435	lmix_bis(ig)=l+1
436	endif
437	lmix(ig)=l+1
438	wmaxa(ig)=wa_moy(ig,l+1)
439	endif
440	enddo
441
442	if (nbpb>0) then
443	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
444	endif
445
446	!=========================================================================
447	! FIN DE LA BOUCLE VERTICALE
448	enddo
449	!=========================================================================
450
451	!on recalcule alim_star_tot
452	do ig=1,ngrid
453	alim_star_tot(ig)=0.
454	enddo
455	do ig=1,ngrid
456	do l=1,lalim(ig)-1
457	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
458	enddo
459	enddo
460
461
462	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
463
464	#undef wrgrads_thermcell
465	#ifdef wrgrads_thermcell
466	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
467	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
468	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
469	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
470	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
471	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
472	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
473	#endif
474
475
476	return
477	end
478
479
480	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
481	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
482	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
483	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
484	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
485	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
486	SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
487	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
488	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
489	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
490	& ,lev_out,lunout1,igout)
491
492	!--------------------------------------------------------------------------
493	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
494	! Version conforme a l'article de Rio et al. 2010.
495	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
496	!--------------------------------------------------------------------------
497
498	IMPLICIT NONE
499
500	#include "YOMCST.h"
501	#include "YOETHF.h"
502	#include "FCTTRE.h"
503	#include "iniprint.h"
504	#include "thermcell.h"
505
506	INTEGER itap
507	INTEGER lunout1,igout
508	INTEGER ngrid,klev
509	REAL ptimestep
510	REAL ztv(ngrid,klev)
511	REAL zthl(ngrid,klev)
512	REAL po(ngrid,klev)
513	REAL zl(ngrid,klev)
514	REAL rhobarz(ngrid,klev)
515	REAL zlev(ngrid,klev+1)
516	REAL pplev(ngrid,klev+1)
517	REAL pphi(ngrid,klev)
518	REAL zpspsk(ngrid,klev)
519	REAL alim_star(ngrid,klev)
520	REAL f0(ngrid)
521	INTEGER lalim(ngrid)
522	integer lev_out ! niveau pour les print
523	integer nbpb
524
525	real alim_star_tot(ngrid)
526
527	REAL ztva(ngrid,klev)
528	REAL ztla(ngrid,klev)
529	REAL zqla(ngrid,klev)
530	REAL zqta(ngrid,klev)
531	REAL zha(ngrid,klev)
532
533	REAL detr_star(ngrid,klev)
534	REAL coefc
535	REAL entr_star(ngrid,klev)
536	REAL detr(ngrid,klev)
537	REAL entr(ngrid,klev)
538
539	REAL csc(ngrid,klev)
540
541	REAL zw2(ngrid,klev+1)
542	REAL w_est(ngrid,klev+1)
543	REAL f_star(ngrid,klev+1)
544	REAL wa_moy(ngrid,klev+1)
545
546	REAL ztva_est(ngrid,klev)
547	REAL zqla_est(ngrid,klev)
548	REAL zqsatth(ngrid,klev)
549	REAL zta_est(ngrid,klev)
550	REAL ztemp(ngrid),zqsat(ngrid)
551	REAL zdw2
552	REAL zw2modif
553	REAL zw2fact
554	REAL zeps(ngrid,klev)
555
556	REAL linter(ngrid)
557	INTEGER lmix(ngrid)
558	INTEGER lmix_bis(ngrid)
559	REAL wmaxa(ngrid)
560
561	INTEGER ig,l,k
562
563	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
564	real zbuoybis
565	real zcor,zdelta,zcvm5,qlbef,zdz2
566	real betalpha,zbetalpha
567	real eps, afact
568	REAL REPS,RLvCp,DDT0
569	PARAMETER (DDT0=.01)
570	logical Zsat
571	LOGICAL active(ngrid),activetmp(ngrid)
572	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
573	REAL c2(ngrid,klev)
574	Zsat=.false.
575	! Initialisation
576
577	RLvCp = RLVTT/RCPD
578	fact_epsilon=0.002
579	betalpha=0.9
580	afact=2./3.
581
582	zbetalpha=betalpha/(1.+betalpha)
583
584
585	! Initialisations des variables reeles
586	if (1==0) then
587	ztva(:,:)=ztv(:,:)
588	ztva_est(:,:)=ztva(:,:)
589	ztla(:,:)=zthl(:,:)
590	zqta(:,:)=po(:,:)
591	zha(:,:) = ztva(:,:)
592	else
593	ztva(:,:)=0.
594	ztva_est(:,:)=0.
595	ztla(:,:)=0.
596	zqta(:,:)=0.
597	zha(:,:) =0.
598	endif
599
600	zqla_est(:,:)=0.
601	zqsatth(:,:)=0.
602	zqla(:,:)=0.
603	detr_star(:,:)=0.
604	entr_star(:,:)=0.
605	alim_star(:,:)=0.
606	alim_star_tot(:)=0.
607	csc(:,:)=0.
608	detr(:,:)=0.
609	entr(:,:)=0.
610	zw2(:,:)=0.
611	zbuoy(:,:)=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	endif
834
835	if (zw2(ig,l+1).lt.0.) then
836	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
837	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
838	zw2(ig,l+1)=0.
839	endif
840
841	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
842
843	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
844	! lmix est le niveau de la couche ou w (wa_moy) est maximum
845	!on rajoute le calcul de lmix_bis
846	if (zqla(ig,l).lt.1.e-10) then
847	lmix_bis(ig)=l+1
848	endif
849	lmix(ig)=l+1
850	wmaxa(ig)=wa_moy(ig,l+1)
851	endif
852	enddo
853
854	if (nbpb>0) then
855	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
856	endif
857
858	!=========================================================================
859	! FIN DE LA BOUCLE VERTICALE
860	enddo
861	!=========================================================================
862
863	!on recalcule alim_star_tot
864	do ig=1,ngrid
865	alim_star_tot(ig)=0.
866	enddo
867	do ig=1,ngrid
868	do l=1,lalim(ig)-1
869	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
870	enddo
871	enddo
872
873
874	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
875
876	#undef wrgrads_thermcell
877	#ifdef wrgrads_thermcell
878	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
879	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
880	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
881	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
882	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
883	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
884	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
885	#endif
886
887
888	return
889	end

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

Download in other formats:

Original Format