Context Navigation

lmdz_thermcell_plume_6A.F90 @ 5407

Last change on this file since 5407 was 5390, checked in by yann meurdesoif, 9 days ago

Remove UTF8 character that inihibit fortran parsing with GPU morphosis
Add missing END SUBROUTINE instead of simple END, that inhibit correct parsing with regulat expression parser (quick and dirty parsing)

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.8 KB

Line
1	MODULE lmdz_thermcell_plume_6A
2	!
3	! $Id: lmdz_thermcell_plume_6A.F90 5390 2024-12-05 16:09:25Z abarral $
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 SUBROUTINE thermcell_plume_6A
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 SUBROUTINE thermcell_plume_5B
1118	END MODULE lmdz_thermcell_plume_6A

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

Download in other formats:

Original Format