Context Navigation

thermcell_plume.F90 @ 2093

Last change on this file since 2093 was 2046, checked in by fhourdin, 11 years ago

Nouvelles options des thermiques en route vers les stratocumulus

+ sorties de la discretisation verticale dans le vieux disvert_old

New options in thermals for stratocumulus

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

Line
1	!
2	! $Id: thermcell_plume.F90 2046 2014-05-18 13:18:19Z 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 fact_shell
83	real ztv1,ztv2,factinv,zinv,zlmel
84	real ztv_est1,ztv_est2
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	fact_shell=0.85
102
103	zbetalpha=betalpha/(1.+betalpha)
104
105
106	! Initialisations des variables r?elles
107	if (1==1) then
108	ztva(:,:)=ztv(:,:)
109	ztva_est(:,:)=ztva(:,:)
110	ztv_est(:,:)=ztv(:,:)
111	ztla(:,:)=zthl(:,:)
112	zqta(:,:)=po(:,:)
113	zqla(:,:)=0.
114	zha(:,:) = ztva(:,:)
115	else
116	ztva(:,:)=0.
117	ztv_est(:,:)=0.
118	ztva_est(:,:)=0.
119	ztla(:,:)=0.
120	zqta(:,:)=0.
121	zha(:,:) =0.
122	endif
123
124	zqla_est(:,:)=0.
125	zqsatth(:,:)=0.
126	zqla(:,:)=0.
127	detr_star(:,:)=0.
128	entr_star(:,:)=0.
129	alim_star(:,:)=0.
130	alim_star_tot(:)=0.
131	csc(:,:)=0.
132	detr(:,:)=0.
133	entr(:,:)=0.
134	zw2(:,:)=0.
135	zbuoy(:,:)=0.
136	zbuoyjam(:,:)=0.
137	gamma(:,:)=0.
138	zeps(:,:)=0.
139	w_est(:,:)=0.
140	f_star(:,:)=0.
141	wa_moy(:,:)=0.
142	linter(:)=1.
143	! linter(:)=1.
144	! Initialisation des variables entieres
145	lmix(:)=1
146	lmix_bis(:)=2
147	wmaxa(:)=0.
148	lalim(:)=1
149
150
151	!-------------------------------------------------------------------------
152	! On ne considere comme actif que les colonnes dont les deux premieres
153	! couches sont instables.
154	!-------------------------------------------------------------------------
155	active(:)=ztv(:,1)>ztv(:,2)
156
157	!-------------------------------------------------------------------------
158	! Definition de l'alimentation a l'origine dans thermcell_init
159	!-------------------------------------------------------------------------
160	do l=1,klev-1
161	do ig=1,ngrid
162	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
163	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
164	& *sqrt(zlev(ig,l+1))
165	lalim(ig)=l+1
166	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
167	! print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)
168	endif
169	enddo
170	enddo
171	do l=1,klev
172	do ig=1,ngrid
173	if (alim_star_tot(ig) > 1.e-10 ) then
174	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
175	endif
176	enddo
177	enddo
178	alim_star_tot(:)=1.
179
180
181
182
183	!------------------------------------------------------------------------------
184	! Calcul dans la premiere couche
185	! On decide dans cette version que le thermique n'est actif que si la premiere
186	! couche est instable.
187	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
188	! dans une couche l>1
189	!------------------------------------------------------------------------------
190	do ig=1,ngrid
191	! Le panache va prendre au debut les caracteristiques de l'air contenu
192	! dans cette couche.
193	if (active(ig)) then
194	ztla(ig,1)=zthl(ig,1)
195	zqta(ig,1)=po(ig,1)
196	zqla(ig,1)=zl(ig,1)
197	!cr: attention, prise en compte de f*(1)=1
198	f_star(ig,2)=alim_star(ig,1)
199	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
200	& *(zlev(ig,2)-zlev(ig,1)) &
201	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
202	w_est(ig,2)=zw2(ig,2)
203	endif
204	enddo
205	!
206
207	!==============================================================================
208	!boucle de calcul de la vitesse verticale dans le thermique
209	!==============================================================================
210	do l=2,klev-1
211	!==============================================================================
212
213
214	! On decide si le thermique est encore actif ou non
215	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
216	do ig=1,ngrid
217	active(ig)=active(ig) &
218	& .and. zw2(ig,l)>1.e-10 &
219	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
220	enddo
221
222
223
224	!---------------------------------------------------------------------------
225	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
226	! sans tenir compte du detrainement et de l'entrainement dans cette
227	! couche
228	! C'est a dire qu'on suppose
229	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
230	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
231	! avant) a l'alimentation pour avoir un calcul plus propre
232	!---------------------------------------------------------------------------
233
234	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
235	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
236	do ig=1,ngrid
237	! print*,'active',active(ig),ig,l
238	if(active(ig)) then
239	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
240	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
241	zta_est(ig,l)=ztva_est(ig,l)
242	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
243	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
244	& -zqla_est(ig,l))-zqla_est(ig,l))
245
246
247	!------------------------------------------------
248	!AJAM:nouveau calcul de w?
249	!------------------------------------------------
250	zdz=zlev(ig,l+1)-zlev(ig,l)
251	zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
252	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
253
254	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
255	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
256	zdw2=afact*zbuoy(ig,l)/fact_epsilon
257	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
258	! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
259	! w_est(ig,l+1)=(zdz/zdzbis)Max(0.0001,exp(-zw2fact) &
260	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
261	! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
262	! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))zdw2+w_est(ig,l)exp(-zw2fact))
263
264	!--------------------------------------------------
265	!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
266	!--------------------------------------------------
267	if (iflag_thermals_ed==8) then
268	! Ancienne version
269	w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
270	& (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
271	& (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
272
273	! Nouvelle version Arnaud
274	else
275	w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
276	& (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
277	& (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
278	endif
279
280
281	if (iflag_thermals_ed<6) then
282	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
283	! fact_epsilon=0.0005/(zalpha+0.025)**0.5
284	! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
285	fact_epsilon=0.0002/(zalpha+0.1)
286	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
287	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
288	zdw2=afact*zbuoy(ig,l)/fact_epsilon
289	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
290	w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
291	& (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
292	& (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
293
294	endif
295	!--------------------------------------------------
296	!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
297	!on fait max(0.0001,.....)
298	!--------------------------------------------------
299
300	! if (w_est(ig,l+1).lt.0.) then
301	! w_est(ig,l+1)=zw2(ig,l)
302	! w_est(ig,l+1)=0.0001
303	! endif
304
305	endif
306	enddo
307
308
309	!-------------------------------------------------
310	!calcul des taux d'entrainement et de detrainement
311	!-------------------------------------------------
312
313	do ig=1,ngrid
314	if (active(ig)) then
315
316	! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
317	zw2m=w_est(ig,l+1)
318	! zw2m=zw2(ig,l)
319	zdz=zlev(ig,l+1)-zlev(ig,l)
320	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
321	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
322	zbuoybis=zbuoy(ig,l)
323	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
324	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
325
326
327	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
328	! & afact*zbuoybis/zw2m - fact_epsilon )
329
330	! entr_star(ig,l)=MAX(0.,f_star(ig,l)zdzzbetalpha* &
331	! & afact*zbuoybis/zw2m - fact_epsilon )
332
333	!Modif AJAM
334
335	lmel=fact_thermals_ed_dz*zlev(ig,l)
336	zlmel=zlev(ig,l)+lmel
337	! lmel=2.5*(zlev(ig,l)-zlev(ig,l-1))
338	lt=l+1
339	zdz2=zlev(ig,lt)-zlev(ig,l)
340	!--------------------------------------------------
341	!AJ052014: J'ai remplac? la boucle do par un do while
342	! afin de faire moins de calcul dans la boucle
343	!--------------------------------------------------
344	do while (zdz2.lt.lmel)
345	lt=lt+1
346	zdz2=zlev(ig,lt)-zlev(ig,l)
347	end do
348
349	zdz3=zlev(ig,lt)-zlev(ig,lt-1)
350
351	!--------------------------------------------------
352	!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
353	! on cherche o? se trouve l'altitude d'inversion
354	! en calculant ztv1 (interpolation de la valeur de
355	! theta au niveau lt en utilisant les niveaux lt-1 et
356	! lt-2) et ztv2 (interpolation avec les niveaux lt+1
357	! et lt+2). Si theta r?ellement calcul?e au niveau lt
358	! comprise entre ztv1 et ztv2, alors il y a inversion
359	! et on calcule son altitude zinv en supposant que ztv(lt)
360	! est une combinaison lin?aire de ztv1 et ztv2.
361	! Ensuite, on calcule la flottabilit? en comparant
362	! la temp?rature de la couche l ? celle de l'air situ?
363	! l+lmel plus haut, ce qui n?cessite de savoir quel fraction
364	! de cet air est au-dessus ou en-dessous de l'inversion
365	!--------------------------------------------------
366
367
368	if (iflag_thermals_ed.lt.8) then
369
370	ztv1=(ztv(ig,lt-1)-ztv(ig,lt-2))*zlev(ig,lt)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
371	& +(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
372	& /(zlev(ig,lt-1)-zlev(ig,lt-2))
373
374	ztv2=(ztv(ig,lt+2)-ztv(ig,lt+1))*zlev(ig,lt)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
375	& +(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
376	& /(zlev(ig,lt+2)-zlev(ig,lt+1))
377
378	if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
379
380	factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
381	zinv=zlev(ig,lt-1)+factinv*(zlev(ig,lt)-zlev(ig,lt-1))
382
383	if (zlmel+0.5*zdz.ge.zinv) then
384	if (zlmel-0.5*zdz.ge.zinv) then
385
386	ztv_est(ig,l)=(ztv(ig,lt+2)-ztv(ig,lt+1))(zlmel-0.zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
387	& +(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
388	& /(zlev(ig,lt+2)-zlev(ig,lt+1))
389
390	zbuoyjam(ig,l)=fact_shellRG(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l)
391
392	else
393
394	ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))0.5(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
395	& +(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
396	& /(zlev(ig,lt+2)-zlev(ig,lt+1))
397	ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))0.5(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
398	& +(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
399	& /(zlev(ig,lt-1)-zlev(ig,lt-2))
400	zbuoyjam(ig,l)=fact_shellRG(((zlmel+0.5zdz-zinv)/zdz)(ztva_est(ig,l)- &
401	& ztv_est1)/ztv_est1+((zinv-zlmel+0.5zdz)/zdz)(ztva_est(ig,l)- &
402	& ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l)
403
404	endif
405
406	else
407
408	ztv_est(ig,l)=(ztv(ig,lt-1)-ztv(ig,lt-2))(zlmel-0.zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
409	& +(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
410	& /(zlev(ig,lt-1)-zlev(ig,lt-2))
411
412	zbuoyjam(ig,l)=fact_shellRG(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l)
413	! ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))0.5(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
414	! & +(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
415	! & /(zlev(ig,lt+2)-zlev(ig,lt+1))
416	! ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))0.5(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
417	! & +(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
418	! & /(zlev(ig,lt-1)-zlev(ig,lt-2))
419	! zbuoyjam(ig,l)=fact_shellRG(((zlmel+0.5zdz-zinv)/zdz)(ztva_est(ig,l)- &
420	! & ztv_est1)/ztv_est1+((zinv-zlmel+0.5zdz)/zdz)(ztva_est(ig,l)- &
421	! & ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l)
422
423
424
425
426	! print*,'on est pass? par l?',l,lt,zbuoyjam(ig,l),zbuoy(ig,l)
427	endif
428
429
430	else
431
432	! ztv_est(ig,l)=(lmel/zdz2)*(ztv(ig,lt)-ztv(ig,l))+ztv(ig,l)
433	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)
434
435	zbuoyjam(ig,l)=fact_shellRG(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
436	& ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
437	& ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
438
439	! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
440	! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
441	! & po(ig,lt-1))/po(ig,lt-1))
442
443	endif
444
445	else
446
447	zbuoyjam(ig,l)=1.RG(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
448	& ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
449	& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
450
451	endif
452
453
454	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
455
456	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
457	! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
458
459	! entrbis=entr_star(ig,l)
460
461	if (iflag_thermals_ed.lt.6) then
462	fact_epsilon=0.0002/(zalpha+0.1)
463	endif
464
465	detr_star(ig,l)=f_star(ig,l)*zdz &
466	& MAX(1.e-4, -afactzbetalpha*zbuoyjam(ig,l)/zw2m &
467	& + 0.012(zdqt(ig,l)/zw2m)*0.5)
468
469	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
470
471	entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
472	& afact*zbuoy(ig,l)/zw2m - fact_epsilon)
473
474	! entr_star(ig,l)=Max(0.,f_star(ig,l)zdzzbetalpha* &
475	! & afact*zbuoy(ig,l)/zw2m &
476	! & - 1.*fact_epsilon)
477
478
479	! En dessous de lalim, on prend le max de alim_star et entr_star pour
480	! alim_star et 0 sinon
481	if (l.lt.lalim(ig)) then
482	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
483	entr_star(ig,l)=0.
484	endif
485	! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
486	! alim_star(ig,l)=entrbis
487	! endif
488
489	! print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
490	! Calcul du flux montant normalise
491	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
492	& -detr_star(ig,l)
493
494	endif
495	enddo
496
497
498	!----------------------------------------------------------------------------
499	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
500	!---------------------------------------------------------------------------
501	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
502	do ig=1,ngrid
503	if (activetmp(ig)) then
504	Zsat=.false.
505	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
506	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
507	& /(f_star(ig,l+1)+detr_star(ig,l))
508	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
509	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
510	& /(f_star(ig,l+1)+detr_star(ig,l))
511
512	endif
513	enddo
514
515	ztemp(:)=zpspsk(:,l)*ztla(:,l)
516	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
517	do ig=1,ngrid
518	if (activetmp(ig)) then
519	! on ecrit de maniere conservative (sat ou non)
520	! T = Tl +Lv/Cp ql
521	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
522	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
523	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
524	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
525	zha(ig,l) = ztva(ig,l)
526	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
527	& -zqla(ig,l))-zqla(ig,l))
528	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
529	zdz=zlev(ig,l+1)-zlev(ig,l)
530	zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
531	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
532	fact_epsilon=0.002
533	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
534	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
535	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
536	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
537	! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
538	zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
539	& (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
540	& (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
541
542	if (iflag_thermals_ed.lt.6) then
543	zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
544	! fact_epsilon=0.0005/(zalpha+0.025)**0.5
545	! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
546	fact_epsilon=0.0002/(zalpha+0.1)**1
547	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
548	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
549	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
550	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
551
552	zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
553	& (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
554	& (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
555
556	endif
557
558
559	endif
560	enddo
561
562	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
563	!
564	!---------------------------------------------------------------------------
565	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
566	!---------------------------------------------------------------------------
567
568	nbpb=0
569	do ig=1,ngrid
570	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
571	! stop'On tombe sur le cas particulier de thermcell_dry'
572	! print*,'On tombe sur le cas particulier de thermcell_plume'
573	nbpb=nbpb+1
574	zw2(ig,l+1)=0.
575	linter(ig)=l+1
576	endif
577
578	if (zw2(ig,l+1).lt.0.) then
579	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
580	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
581	zw2(ig,l+1)=0.
582	!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
583	elseif (f_star(ig,l+1).lt.0.) then
584	linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
585	& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
586	zw2(ig,l+1)=0.
587	!fin CR:04/05/12
588	endif
589
590	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
591
592	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
593	! lmix est le niveau de la couche ou w (wa_moy) est maximum
594	!on rajoute le calcul de lmix_bis
595	if (zqla(ig,l).lt.1.e-10) then
596	lmix_bis(ig)=l+1
597	endif
598	lmix(ig)=l+1
599	wmaxa(ig)=wa_moy(ig,l+1)
600	endif
601	enddo
602
603	if (nbpb>0) then
604	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
605	endif
606
607	!=========================================================================
608	! FIN DE LA BOUCLE VERTICALE
609	enddo
610	!=========================================================================
611
612	!on recalcule alim_star_tot
613	do ig=1,ngrid
614	alim_star_tot(ig)=0.
615	enddo
616	do ig=1,ngrid
617	do l=1,lalim(ig)-1
618	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
619	enddo
620	enddo
621
622
623	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
624
625	#undef wrgrads_thermcell
626	#ifdef wrgrads_thermcell
627	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
628	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
629	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
630	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
631	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
632	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
633	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
634	#endif
635
636
637	return
638	end
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
675	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
676	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
677	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
678	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
679	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
680	SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
681	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
682	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
683	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
684	& ,lev_out,lunout1,igout)
685
686	!--------------------------------------------------------------------------
687	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
688	! Version conforme a l'article de Rio et al. 2010.
689	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
690	!--------------------------------------------------------------------------
691
692	IMPLICIT NONE
693
694	#include "YOMCST.h"
695	#include "YOETHF.h"
696	#include "FCTTRE.h"
697	#include "iniprint.h"
698	#include "thermcell.h"
699
700	INTEGER itap
701	INTEGER lunout1,igout
702	INTEGER ngrid,klev
703	REAL ptimestep
704	REAL ztv(ngrid,klev)
705	REAL zthl(ngrid,klev)
706	REAL po(ngrid,klev)
707	REAL zl(ngrid,klev)
708	REAL rhobarz(ngrid,klev)
709	REAL zlev(ngrid,klev+1)
710	REAL pplev(ngrid,klev+1)
711	REAL pphi(ngrid,klev)
712	REAL zpspsk(ngrid,klev)
713	REAL alim_star(ngrid,klev)
714	REAL f0(ngrid)
715	INTEGER lalim(ngrid)
716	integer lev_out ! niveau pour les print
717	integer nbpb
718
719	real alim_star_tot(ngrid)
720
721	REAL ztva(ngrid,klev)
722	REAL ztla(ngrid,klev)
723	REAL zqla(ngrid,klev)
724	REAL zqta(ngrid,klev)
725	REAL zha(ngrid,klev)
726
727	REAL detr_star(ngrid,klev)
728	REAL coefc
729	REAL entr_star(ngrid,klev)
730	REAL detr(ngrid,klev)
731	REAL entr(ngrid,klev)
732
733	REAL csc(ngrid,klev)
734
735	REAL zw2(ngrid,klev+1)
736	REAL w_est(ngrid,klev+1)
737	REAL f_star(ngrid,klev+1)
738	REAL wa_moy(ngrid,klev+1)
739
740	REAL ztva_est(ngrid,klev)
741	REAL zqla_est(ngrid,klev)
742	REAL zqsatth(ngrid,klev)
743	REAL zta_est(ngrid,klev)
744	REAL zbuoyjam(ngrid,klev)
745	REAL ztemp(ngrid),zqsat(ngrid)
746	REAL zdw2
747	REAL zw2modif
748	REAL zw2fact
749	REAL zeps(ngrid,klev)
750
751	REAL linter(ngrid)
752	INTEGER lmix(ngrid)
753	INTEGER lmix_bis(ngrid)
754	REAL wmaxa(ngrid)
755
756	INTEGER ig,l,k
757
758	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
759	real zbuoybis
760	real zcor,zdelta,zcvm5,qlbef,zdz2
761	real betalpha,zbetalpha
762	real eps, afact
763	REAL REPS,RLvCp,DDT0
764	PARAMETER (DDT0=.01)
765	logical Zsat
766	LOGICAL active(ngrid),activetmp(ngrid)
767	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
768	REAL c2(ngrid,klev)
769	Zsat=.false.
770	! Initialisation
771
772	RLvCp = RLVTT/RCPD
773	fact_epsilon=0.002
774	betalpha=0.9
775	afact=2./3.
776
777	zbetalpha=betalpha/(1.+betalpha)
778
779
780	! Initialisations des variables reeles
781	if (1==1) then
782	ztva(:,:)=ztv(:,:)
783	ztva_est(:,:)=ztva(:,:)
784	ztla(:,:)=zthl(:,:)
785	zqta(:,:)=po(:,:)
786	zha(:,:) = ztva(:,:)
787	else
788	ztva(:,:)=0.
789	ztva_est(:,:)=0.
790	ztla(:,:)=0.
791	zqta(:,:)=0.
792	zha(:,:) =0.
793	endif
794
795	zqla_est(:,:)=0.
796	zqsatth(:,:)=0.
797	zqla(:,:)=0.
798	detr_star(:,:)=0.
799	entr_star(:,:)=0.
800	alim_star(:,:)=0.
801	alim_star_tot(:)=0.
802	csc(:,:)=0.
803	detr(:,:)=0.
804	entr(:,:)=0.
805	zw2(:,:)=0.
806	zbuoy(:,:)=0.
807	zbuoyjam(:,:)=0.
808	gamma(:,:)=0.
809	zeps(:,:)=0.
810	w_est(:,:)=0.
811	f_star(:,:)=0.
812	wa_moy(:,:)=0.
813	linter(:)=1.
814	! linter(:)=1.
815	! Initialisation des variables entieres
816	lmix(:)=1
817	lmix_bis(:)=2
818	wmaxa(:)=0.
819	lalim(:)=1
820
821
822	!-------------------------------------------------------------------------
823	! On ne considere comme actif que les colonnes dont les deux premieres
824	! couches sont instables.
825	!-------------------------------------------------------------------------
826	active(:)=ztv(:,1)>ztv(:,2)
827
828	!-------------------------------------------------------------------------
829	! Definition de l'alimentation a l'origine dans thermcell_init
830	!-------------------------------------------------------------------------
831	do l=1,klev-1
832	do ig=1,ngrid
833	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
834	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
835	& *sqrt(zlev(ig,l+1))
836	lalim(ig)=l+1
837	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
838	endif
839	enddo
840	enddo
841	do l=1,klev
842	do ig=1,ngrid
843	if (alim_star_tot(ig) > 1.e-10 ) then
844	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
845	endif
846	enddo
847	enddo
848	alim_star_tot(:)=1.
849
850
851
852	!------------------------------------------------------------------------------
853	! Calcul dans la premiere couche
854	! On decide dans cette version que le thermique n'est actif que si la premiere
855	! couche est instable.
856	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
857	! dans une couche l>1
858	!------------------------------------------------------------------------------
859	do ig=1,ngrid
860	! Le panache va prendre au debut les caracteristiques de l'air contenu
861	! dans cette couche.
862	if (active(ig)) then
863	ztla(ig,1)=zthl(ig,1)
864	zqta(ig,1)=po(ig,1)
865	zqla(ig,1)=zl(ig,1)
866	!cr: attention, prise en compte de f*(1)=1
867	f_star(ig,2)=alim_star(ig,1)
868	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
869	& *(zlev(ig,2)-zlev(ig,1)) &
870	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
871	w_est(ig,2)=zw2(ig,2)
872	endif
873	enddo
874	!
875
876	!==============================================================================
877	!boucle de calcul de la vitesse verticale dans le thermique
878	!==============================================================================
879	do l=2,klev-1
880	!==============================================================================
881
882
883	! On decide si le thermique est encore actif ou non
884	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
885	do ig=1,ngrid
886	active(ig)=active(ig) &
887	& .and. zw2(ig,l)>1.e-10 &
888	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
889	enddo
890
891
892
893	!---------------------------------------------------------------------------
894	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
895	! sans tenir compte du detrainement et de l'entrainement dans cette
896	! couche
897	! C'est a dire qu'on suppose
898	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
899	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
900	! avant) a l'alimentation pour avoir un calcul plus propre
901	!---------------------------------------------------------------------------
902
903	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
904	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
905
906	do ig=1,ngrid
907	! print*,'active',active(ig),ig,l
908	if(active(ig)) then
909	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
910	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
911	zta_est(ig,l)=ztva_est(ig,l)
912	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
913	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
914	& -zqla_est(ig,l))-zqla_est(ig,l))
915
916	!------------------------------------------------
917	!AJAM:nouveau calcul de w?
918	!------------------------------------------------
919	zdz=zlev(ig,l+1)-zlev(ig,l)
920	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
921
922	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
923	zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
924	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
925
926
927	if (w_est(ig,l+1).lt.0.) then
928	w_est(ig,l+1)=zw2(ig,l)
929	endif
930	endif
931	enddo
932
933
934	!-------------------------------------------------
935	!calcul des taux d'entrainement et de detrainement
936	!-------------------------------------------------
937
938	do ig=1,ngrid
939	if (active(ig)) then
940
941	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
942	zw2m=w_est(ig,l+1)
943	zdz=zlev(ig,l+1)-zlev(ig,l)
944	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
945	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
946	zbuoybis=zbuoy(ig,l)
947	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
948	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
949
950
951	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
952	& afact*zbuoybis/zw2m - fact_epsilon )
953
954
955	detr_star(ig,l)=f_star(ig,l)*zdz &
956	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
957	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
958
959	! En dessous de lalim, on prend le max de alim_star et entr_star pour
960	! alim_star et 0 sinon
961	if (l.lt.lalim(ig)) then
962	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
963	entr_star(ig,l)=0.
964	endif
965
966	! Calcul du flux montant normalise
967	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
968	& -detr_star(ig,l)
969
970	endif
971	enddo
972
973
974	!----------------------------------------------------------------------------
975	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
976	!---------------------------------------------------------------------------
977	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
978	do ig=1,ngrid
979	if (activetmp(ig)) then
980	Zsat=.false.
981	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
982	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
983	& /(f_star(ig,l+1)+detr_star(ig,l))
984	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
985	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
986	& /(f_star(ig,l+1)+detr_star(ig,l))
987
988	endif
989	enddo
990
991	ztemp(:)=zpspsk(:,l)*ztla(:,l)
992	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
993
994	do ig=1,ngrid
995	if (activetmp(ig)) then
996	! on ecrit de maniere conservative (sat ou non)
997	! T = Tl +Lv/Cp ql
998	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
999	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
1000	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1001	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1002	zha(ig,l) = ztva(ig,l)
1003	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
1004	& -zqla(ig,l))-zqla(ig,l))
1005	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1006	zdz=zlev(ig,l+1)-zlev(ig,l)
1007	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
1008
1009	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
1010	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
1011	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
1012	endif
1013	enddo
1014
1015	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
1016	!
1017	!---------------------------------------------------------------------------
1018	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
1019	!---------------------------------------------------------------------------
1020
1021	nbpb=0
1022	do ig=1,ngrid
1023	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
1024	! stop'On tombe sur le cas particulier de thermcell_dry'
1025	! print*,'On tombe sur le cas particulier de thermcell_plume'
1026	nbpb=nbpb+1
1027	zw2(ig,l+1)=0.
1028	linter(ig)=l+1
1029	endif
1030
1031	if (zw2(ig,l+1).lt.0.) then
1032	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
1033	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1034	zw2(ig,l+1)=0.
1035	endif
1036
1037	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1038
1039	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1040	! lmix est le niveau de la couche ou w (wa_moy) est maximum
1041	!on rajoute le calcul de lmix_bis
1042	if (zqla(ig,l).lt.1.e-10) then
1043	lmix_bis(ig)=l+1
1044	endif
1045	lmix(ig)=l+1
1046	wmaxa(ig)=wa_moy(ig,l+1)
1047	endif
1048	enddo
1049
1050	if (nbpb>0) then
1051	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
1052	endif
1053
1054	!=========================================================================
1055	! FIN DE LA BOUCLE VERTICALE
1056	enddo
1057	!=========================================================================
1058
1059	!on recalcule alim_star_tot
1060	do ig=1,ngrid
1061	alim_star_tot(ig)=0.
1062	enddo
1063	do ig=1,ngrid
1064	do l=1,lalim(ig)-1
1065	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
1066	enddo
1067	enddo
1068
1069
1070	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
1071
1072	#undef wrgrads_thermcell
1073	#ifdef wrgrads_thermcell
1074	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
1075	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
1076	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
1077	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
1078	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
1079	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
1080	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
1081	#endif
1082
1083
1084	return
1085	end

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

Download in other formats:

Original Format