Context Navigation

lmdz_thermcell_plume_6A.f90 @ 5901

Last change on this file since 5901 was 5854, checked in by fhourdin, 8 weeks ago
Separation thermcell_plume 5B/6A
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: 26.4 KB

Line
1	MODULE lmdz_thermcell_plume_6A
2	!
3	! $Id: lmdz_thermcell_plume_6A.f90 5854 2025-10-10 14:00:56Z idelkadi $
4	!
5	CONTAINS
6
7	SUBROUTINE thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
8	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
9	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
10	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
11	& ,lev_out,lunout1,igout)
12	! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
13	!--------------------------------------------------------------------------
14	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
15	!--------------------------------------------------------------------------
16
17	USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
18	USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
19	USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
20	USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
21	USE lmdz_thermcell_alim, ONLY : thermcell_alim
22	USE lmdz_thermcell_qsat, ONLY : thermcell_qsat
23
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
101	zbetalpha=betalpha/(1.+betalpha)
102
103
104	! Initialisations des variables r?elles
105	if (1==1) then
106	ztva(:,:)=ztv(:,:)
107	ztva_est(:,:)=ztva(:,:)
108	ztv_est(:,:)=ztv(:,:)
109	ztla(:,:)=zthl(:,:)
110	zqta(:,:)=po(:,:)
111	zqla(:,:)=0.
112	zha(:,:) = ztva(:,:)
113	else
114	ztva(:,:)=0.
115	ztv_est(:,:)=0.
116	ztva_est(:,:)=0.
117	ztla(:,:)=0.
118	zqta(:,:)=0.
119	zha(:,:) =0.
120	endif
121
122	zqla_est(:,:)=0.
123	zqsatth(:,:)=0.
124	zqla(:,:)=0.
125	detr_star(:,:)=0.
126	entr_star(:,:)=0.
127	alim_star(:,:)=0.
128	alim_star_tot(:)=0.
129	csc(:,:)=0.
130	detr(:,:)=0.
131	entr(:,:)=0.
132	zw2(:,:)=0.
133	zbuoy(:,:)=0.
134	zbuoyjam(:,:)=0.
135	gamma(:,:)=0.
136	zeps(:,:)=0.
137	w_est(:,:)=0.
138	f_star(:,:)=0.
139	wa_moy(:,:)=0.
140	linter(:)=1.
141	! linter(:)=1.
142	! Initialisation des variables entieres
143	lmix(:)=1
144	lmix_bis(:)=2
145	wmaxa(:)=0.
146
147	! Initialisation a 0 en cas de sortie dans replay
148	zqsat(:)=0.
149	zta_est(:,:)=0.
150	zdqt(:,:)=0.
151	zdqtjam(:,:)=0.
152	c2(:,:)=0.
153
154
155	!-------------------------------------------------------------------------
156	! On ne considere comme actif que les colonnes dont les deux premieres
157	! couches sont instables.
158	!-------------------------------------------------------------------------
159
160	active(:)=ztv(:,1)>ztv(:,2)
161	d_temp(:)=0. ! Pour activer un contraste de temperature a la base
162	! du panache
163	! Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
164	CALL thermcell_alim(thermals_flag_alim,ngrid,nlay,ztv,d_temp,zlev,alim_star,lalim)
165
166	!------------------------------------------------------------------------------
167	! Calcul dans la premiere couche
168	! On decide dans cette version que le thermique n'est actif que si la premiere
169	! couche est instable.
170	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
171	! dans une couche l>1
172	!------------------------------------------------------------------------------
173	do ig=1,ngrid
174	! Le panache va prendre au debut les caracteristiques de l'air contenu
175	! dans cette couche.
176	if (active(ig)) then
177	ztla(ig,1)=zthl(ig,1)
178	zqta(ig,1)=po(ig,1)
179	zqla(ig,1)=zl(ig,1)
180	!cr: attention, prise en compte de f*(1)=1
181	f_star(ig,2)=alim_star(ig,1)
182	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
183	& *(zlev(ig,2)-zlev(ig,1)) &
184	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
185	w_est(ig,2)=zw2(ig,2)
186	endif
187	enddo
188	!
189
190	!==============================================================================
191	!boucle de calcul de la vitesse verticale dans le thermique
192	!==============================================================================
193	do l=2,nlay-1
194	!==============================================================================
195
196
197	! On decide si le thermique est encore actif ou non
198	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
199	do ig=1,ngrid
200	active(ig)=active(ig) &
201	& .and. zw2(ig,l)>1.e-10 &
202	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
203	enddo
204
205
206
207	!---------------------------------------------------------------------------
208	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
209	! sans tenir compte du detrainement et de l'entrainement dans cette
210	! couche
211	! C'est a dire qu'on suppose
212	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
213	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
214	! avant) a l'alimentation pour avoir un calcul plus propre
215	!---------------------------------------------------------------------------
216
217	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
218	call thermcell_qsat(ngrid, 1, active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
219	do ig=1,ngrid
220	! print*,'active',active(ig),ig,l
221	if(active(ig)) then
222	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
223	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
224	zta_est(ig,l)=ztva_est(ig,l)
225	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
226	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1)-zqla_est(ig,l)))
227
228
229	!Modif AJAM
230
231	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
232	zdz=zlev(ig,l+1)-zlev(ig,l)
233	lmel=fact_thermals_ed_dz*zlev(ig,l)
234	! lmel=0.09*zlev(ig,l)
235	zlmel=zlev(ig,l)+lmel
236	zlmelup=zlmel+(zdz/2)
237	zlmeldwn=zlmel-(zdz/2)
238
239	lt=l+1
240	zlt=zlev(ig,lt)
241	zdz3=zlev(ig,lt+1)-zlt
242	zltdwn=zlt-zdz3/2
243	zltup=zlt+zdz3/2
244
245	!=========================================================================
246	! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
247	!=========================================================================
248
249	!--------------------------------------------------
250	if (iflag_thermals_ed.lt.8) then
251	!--------------------------------------------------
252	!AJ052014: J'ai remplac?? la boucle do par un do while
253	! afin de faire moins de calcul dans la boucle
254	!--------------------------------------------------
255	do while (zlmelup.gt.zltup)
256	lt=lt+1
257	zlt=zlev(ig,lt)
258	zdz3=zlev(ig,lt+1)-zlt
259	zltdwn=zlt-zdz3/2
260	zltup=zlt+zdz3/2
261	enddo
262	!--------------------------------------------------
263	!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
264	! on cherche o?? se trouve l'altitude d'inversion
265	! en calculant ztv1 (interpolation de la valeur de
266	! theta au niveau lt en utilisant les niveaux lt-1 et
267	! lt-2) et ztv2 (interpolation avec les niveaux lt+1
268	! et lt+2). Si theta r??ellement calcul??e au niveau lt
269	! comprise entre ztv1 et ztv2, alors il y a inversion
270	! et on calcule son altitude zinv en supposant que ztv(lt)
271	! est une combinaison lineaire de ztv1 et ztv2.
272	! Ensuite, on calcule la flottabilite en comparant
273	! la temperature de la couche l a celle de l'air situe
274	! l+lmel plus haut, ce qui necessite de savoir quel fraction
275	! de cet air est au-dessus ou en-dessous de l'inversion
276	!--------------------------------------------------
277	atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2))
278	btv1=(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
279	& /(zlev(ig,lt-1)-zlev(ig,lt-2))
280	atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1))
281	btv2=(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
282	& /(zlev(ig,lt+2)-zlev(ig,lt+1))
283
284	ztv1=atv1*zlt+btv1
285	ztv2=atv2*zlt+btv2
286
287	if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
288
289	!--------------------------------------------------
290	!AJ052014: D??calage de zinv qui est entre le haut
291	! et le bas de la couche lt
292	!--------------------------------------------------
293	factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
294	zinv=zltdwn+zdz3*factinv
295
296
297	if (zlmeldwn.ge.zinv) then
298	ztv_est(ig,l)=atv2*zlmel+btv2
299	zbuoyjam(ig,l)=fact_shellRG(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
300	& +(1.-fact_shell)*zbuoy(ig,l)
301	elseif (zlmelup.ge.zinv) then
302	ztv_est2=atv20.5(zlmelup+zinv)+btv2
303	ztv_est1=atv10.5(zinv+zlmeldwn)+btv1
304	ztv_est(ig,l)=((zlmelup-zinv)/zdz)ztv_est2+((zinv-zlmeldwn)/zdz)ztv_est1
305
306	zbuoyjam(ig,l)=fact_shellRG(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- &
307	& ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
308	& ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
309
310	else
311	ztv_est(ig,l)=atv1*zlmel+btv1
312	zbuoyjam(ig,l)=fact_shellRG(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
313	& +(1.-fact_shell)*zbuoy(ig,l)
314	endif
315
316	else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
317
318	if (zlmeldwn.gt.zltdwn) then
319	zbuoyjam(ig,l)=fact_shellRG((ztva_est(ig,l)- &
320	& ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l)
321	else
322	zbuoyjam(ig,l)=fact_shellRG(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
323	& ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
324	& ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
325
326	endif
327
328	! zbuoyjam(ig,l)=fact_shellRG(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
329	! & ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
330	! & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
331	! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
332	! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
333	! & po(ig,lt-1))/po(ig,lt-1))
334	endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
335
336	else ! if (iflag_thermals_ed.lt.8) then
337	lt=l+1
338	zlt=zlev(ig,lt)
339	zdz2=zlev(ig,lt)-zlev(ig,l)
340
341	do while (lmel.gt.zdz2)
342	lt=lt+1
343	zlt=zlev(ig,lt)
344	zdz2=zlt-zlev(ig,l)
345	enddo
346	zdz3=zlev(ig,lt+1)-zlt
347	zltdwn=zlev(ig,lt)-zdz3/2
348	zlmelup=zlmel+(zdz/2)
349	coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
350	zbuoyjam(ig,l)=1.RG(coefzlmel*(ztva_est(ig,l)- &
351	& ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
352	& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
353	endif ! if (iflag_thermals_ed.lt.8) then
354
355	!------------------------------------------------
356	!AJAM:nouveau calcul de w?
357	!------------------------------------------------
358	zdz=zlev(ig,l+1)-zlev(ig,l)
359	zdzbis=zlev(ig,l)-zlev(ig,l-1)
360	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
361
362	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
363	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
364	zdw2=afact*zbuoy(ig,l)/fact_epsilon
365	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
366	! zdw2bis=0.5*(zdw2+zdw2bis)
367	lm=Max(1,l-2)
368	! zdw2=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l) &
369	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
370	! zdw2bis=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l-1) &
371	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
372	! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
373	! w_est(ig,l+1)=(zdz/zdzbis)Max(0.0001,exp(-zw2fact) &
374	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
375	! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
376	! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))zdw2+w_est(ig,l)exp(-zw2fact))
377
378	!--------------------------------------------------
379	!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
380	!--------------------------------------------------
381	if (iflag_thermals_ed==8) then
382	! Ancienne version
383	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
384	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
385	! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
386
387	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
388
389	! Nouvelle version Arnaud
390	else
391	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
392	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
393	! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
394
395	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
396
397	! w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))(exp(-zw2fact) &
398	! & (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
399	! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
400
401
402
403	! w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2biszw2fact)exp(-zw2fact))
404
405	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(zw2(ig,l)+zdw2zw2fact)*exp(-zw2fact)+ &
406	! & (zdzbis-zdz)/zdzbis(zw2(ig,l-1)+zdw2biszw2factbis)*exp(-zw2factbis))
407
408	! w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
409
410	endif
411
412
413	if (iflag_thermals_ed<6) then
414	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
415	! fact_epsilon=0.0005/(zalpha+0.025)**0.5
416	! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
417	fact_epsilon=0.0002/(zalpha+0.1)
418	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
419	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
420	zdw2=afact*zbuoy(ig,l)/fact_epsilon
421	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
422	! w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2zw2fact)exp(-zw2fact))
423
424	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
425	! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
426	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
427
428	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
429
430
431	endif
432	!--------------------------------------------------
433	!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
434	!on fait max(0.0001,.....)
435	!--------------------------------------------------
436
437	! if (w_est(ig,l+1).lt.0.) then
438	! w_est(ig,l+1)=zw2(ig,l)
439	! w_est(ig,l+1)=0.0001
440	! endif
441
442	endif
443	enddo
444
445
446	!-------------------------------------------------
447	!calcul des taux d'entrainement et de detrainement
448	!-------------------------------------------------
449
450	do ig=1,ngrid
451	if (active(ig)) then
452
453	! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
454	zw2m=w_est(ig,l+1)
455	! zw2m=zw2(ig,l)
456	zdz=zlev(ig,l+1)-zlev(ig,l)
457	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
458	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
459	zbuoybis=zbuoy(ig,l)
460	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
461	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
462
463
464	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
465	! & afact*zbuoybis/zw2m - fact_epsilon )
466
467	! entr_star(ig,l)=MAX(0.,f_star(ig,l)zdzzbetalpha* &
468	! & afact*zbuoybis/zw2m - fact_epsilon )
469
470
471
472	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
473
474	!=========================================================================
475	! 4. Calcul de l'entrainement et du detrainement
476	!=========================================================================
477
478	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
479	! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
480	! entrbis=entr_star(ig,l)
481
482	if (iflag_thermals_ed.lt.6) then
483	fact_epsilon=0.0002/(zalpha+0.1)
484	endif
485
486
487
488	detr_star(ig,l)=f_star(ig,l)*zdz &
489	& ( mix0 0.1 / (zalpha+0.001) &
490	& + MAX(detr_min, -afactzbetalphazbuoyjam(ig,l)/zw2m &
491	& + detr_q_coef(zdqt(ig,l)/zw2m)*detr_q_power))
492
493	! detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
494	! & ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
495
496	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
497
498	entr_star(ig,l)=f_star(ig,l)zdz ( &
499	& mix0 * 0.1 / (zalpha+0.001) &
500	& + zbetalpha*MAX(entr_min, &
501	& afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
502
503
504	! entr_star(ig,l)=f_star(ig,l)zdz ( &
505	! & mix0 * 0.1 / (zalpha+0.001) &
506	! & + MAX(entr_min, &
507	! & zbetalphaafactzbuoyjam(ig,l)/zw2m - fact_epsilon + &
508	! & detr_q_coef(zdqt(ig,l)/zw2m)*detr_q_power))
509
510
511	! entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
512	! & ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
513
514	! entr_star(ig,l)=Max(0.,f_star(ig,l)zdzzbetalpha* &
515	! & afact*zbuoy(ig,l)/zw2m &
516	! & - 1.*fact_epsilon)
517
518
519	! En dessous de lalim, on prend le max de alim_star et entr_star pour
520	! alim_star et 0 sinon
521	if (l.lt.lalim(ig)) then
522	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
523	entr_star(ig,l)=0.
524	endif
525	! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
526	! alim_star(ig,l)=entrbis
527	! endif
528
529	! print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
530	! Calcul du flux montant normalise
531	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
532	& -detr_star(ig,l)
533
534	endif
535	enddo
536
537
538	!============================================================================
539	! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
540	!===========================================================================
541
542	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
543	do ig=1,ngrid
544	if (activetmp(ig)) then
545	Zsat=.false.
546	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
547	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
548	& /(f_star(ig,l+1)+detr_star(ig,l))
549	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
550	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
551	& /(f_star(ig,l+1)+detr_star(ig,l))
552
553	endif
554	enddo
555
556	ztemp(:)=zpspsk(:,l)*ztla(:,l)
557	call thermcell_qsat(ngrid, 1, activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
558	do ig=1,ngrid
559	if (activetmp(ig)) then
560	! on ecrit de maniere conservative (sat ou non)
561	! T = Tl +Lv/Cp ql
562	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
563	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
564	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
565	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
566	zha(ig,l) = ztva(ig,l)
567	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l)-zqla(ig,l)))
568	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
569	zdz=zlev(ig,l+1)-zlev(ig,l)
570	zdzbis=zlev(ig,l)-zlev(ig,l-1)
571	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
572	!!!!!!! fact_epsilon=0.002
573	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
574	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
575	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
576	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
577	! zdw2=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l) &
578	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
579	! lm=Max(1,l-2)
580	! zdw2bis=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l-1) &
581	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
582	! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
583	! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(zw2(ig,l)+zdw2zw2fact)*exp(-zw2fact)+ &
584	! & (zdzbis-zdz)/zdzbis(zw2(ig,l-1)+zdw2biszw2factbis)*exp(-zw2factbis))
585	! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2zw2fact)exp(-zw2fact))
586	! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
587	! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
588	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
589	if (iflag_thermals_ed==8) then
590	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
591	else
592	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
593	endif
594	! zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))(exp(-zw2fact) &
595	! & (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
596	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
597
598
599	if (iflag_thermals_ed.lt.6) then
600	zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
601	! fact_epsilon=0.0005/(zalpha+0.025)**0.5
602	! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
603	fact_epsilon=0.0002/(zalpha+0.1)**1
604	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
605	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
606	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
607	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
608
609	! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
610	! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
611	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
612	! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2zw2fact)exp(-zw2fact))
613	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
614
615	endif
616
617
618	endif
619	enddo
620
621	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
622	!
623	!===========================================================================
624	! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
625	!===========================================================================
626
627	nbpb=0
628	do ig=1,ngrid
629	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
630	! stop'On tombe sur le cas particulier de thermcell_dry'
631	! print*,'On tombe sur le cas particulier de thermcell_plume'
632	nbpb=nbpb+1
633	zw2(ig,l+1)=0.
634	linter(ig)=l+1
635	endif
636
637	if (zw2(ig,l+1).lt.0.) then
638	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
639	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
640	zw2(ig,l+1)=0.
641	!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
642	elseif (f_star(ig,l+1).lt.0.) then
643	linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
644	& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
645	zw2(ig,l+1)=0.
646	!fin CR:04/05/12
647	endif
648
649	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
650
651	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
652	! lmix est le niveau de la couche ou w (wa_moy) est maximum
653	!on rajoute le calcul de lmix_bis
654	if (zqla(ig,l).lt.1.e-10) then
655	lmix_bis(ig)=l+1
656	endif
657	lmix(ig)=l+1
658	wmaxa(ig)=wa_moy(ig,l+1)
659	endif
660	enddo
661
662	if (nbpb>0) then
663	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
664	endif
665
666	!=========================================================================
667	! FIN DE LA BOUCLE VERTICALE
668	enddo
669	!=========================================================================
670
671	!on recalcule alim_star_tot
672	do ig=1,ngrid
673	alim_star_tot(ig)=0.
674	enddo
675	do ig=1,ngrid
676	do l=1,lalim(ig)-1
677	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
678	enddo
679	enddo
680
681
682	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
683
684	RETURN
685	END SUBROUTINE thermcell_plume_6A
686	END MODULE lmdz_thermcell_plume_6A

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

Download in other formats:

Original Format