Context Navigation

lmdz_thermcell_plume_6A.F90 @ 5434

Last change on this file since 5434 was 5202, checked in by Laurent Fairhead, 3 months ago
Updating cirrus branch to trunk revision 5171
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: 40.7 KB

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

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

Download in other formats:

Original Format