Context Navigation

thermcell_plume.F90 @ 2387

Last change on this file since 2387 was 2387, checked in by fhourdin, 9 years ago
Poursuite du travail sur les thermiques. Extraction d'une routine calculant l'alimentation laterale à la base. + une correction de Jean-Yves sur des initialisations maquantes dans thermcell_main (mais qui sont maintenant dans thermcell_alp).
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: 41.1 KB

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

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

Download in other formats:

Original Format