Context Navigation

thermcell_plume.F90 @ 2925

Last change on this file since 2925 was 2160, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
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: 42.1 KB

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

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

Download in other formats:

Original Format