Context Navigation

lmdz_thermcell_plume_6A.f90 @ 5105

Last change on this file since 5105 was 5105, checked in by abarral, 8 weeks ago
Replace 1DUTILS.h by module lmdz_1dutils.f90 Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files) Convert .F to .f90 Fix gradsdef.h formatting Remove unnecessary "RETURN" at the end of functions/subroutines
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: 43.1 KB

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

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

Download in other formats:

Original Format