Context Navigation

lmdz_thermcell_plume_6A.f90 @ 5715

Last change on this file since 5715 was 5620, checked in by evignon, 2 months ago
correction du calcul de la temperature virtuelle dans thermcell_plume_6A. Attention! Cette commission conduit une perte de convergence numérique
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: 39.5 KB

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

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

Download in other formats:

Original Format