Context Navigation

thermcell_plume.F90 @ 2106

Last change on this file since 2106 was 2106, checked in by fhourdin, 10 years ago
Version moidifiÃe des thermiques pour une meilleure reprÃ©snation des stratocumulus + contrÃ´le par lecture de .def d'un ertain nombre de paramÃ¨tres internes. Modified version of thermals for stratocumulus + reading new free parameters. Modification by Arnaud Jam and Frederic Hourdin
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.2 KB

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

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

Download in other formats:

Original Format