Context Navigation

thermcell_plume_6A.F90 @ 5420

Last change on this file since 5420 was 3451, checked in by fhourdin, 6 years ago

Saving old version of thermcell_plume in thermcell_plume_6A.F90
It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
to the 5B and 6A versions used for CMIP5 and CMIP6.
The latest was previously named thermcellV1_plume.
The new thermcell_plume is a clean version (removing obsolete
options) of thermcell_plume_6A.
The 3 versions are controled by
flag_thermals_ed <= 9 thermcell_plume_6A

<= 19 thermcell_plume_5B
else thermcell_plume (default 20 for convergence with 6A)

Fredho

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

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

Download in other formats:

Original Format