Context Navigation

thermcell_plume.F90 @ 1311

Last change on this file since 1311 was 1311, checked in by Laurent Fairhead, 14 years ago

Modifications to thermals for TKE transport

Modifications aux thermiques pour le transport de la TKE

pbl_surface_mode.F90 : ok_flux_surf=.false. seulement pour klon>1
physiq.F : option iflag_pbl=10 pour transporter la TKE avec les thermiques.
calltherm.F90 : passage de iflag_thermals_ed en argument pour thermcell_plume
thermcell_main.F90 : Appel a plusieurs version de thermcell_plume en option
thermcell_plume.F90 : plusieurs versions dans le meme fichier (temporaire)
thermcell_height.F90 : verrue pour les cas ou les thermiques montent tout

en haut

yamada4 : inclusion de la diffusion verticale en option iflag_pbl=9

+ variables anciennement common, puis save/allocatable, remises en local

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 28.9 KB

Line
1	!
2	! $Id: thermcell_plume.F90 1311 2010-02-18 13:14:02Z fairhead $
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	!--------------------------------------------------------------------------
11	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
12	!--------------------------------------------------------------------------
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	real zcon2(ngrid)
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 zqla_est(ngrid,klev)
64	REAL zqsatth(ngrid,klev)
65	REAL zta_est(ngrid,klev)
66	REAL zdw2
67	REAL zw2modif
68	REAL zeps
69
70	REAL linter(ngrid)
71	INTEGER lmix(ngrid)
72	INTEGER lmix_bis(ngrid)
73	REAL wmaxa(ngrid)
74
75	INTEGER ig,l,k
76
77	real zdz,zfact,zbuoy,zalpha
78	real zcor,zdelta,zcvm5,qlbef
79	real Tbef,qsatbef
80	real dqsat_dT,DT,num,denom
81	REAL REPS,RLvCp,DDT0
82	PARAMETER (DDT0=.01)
83	logical Zsat
84	LOGICAL active(ngrid),activetmp(ngrid)
85	REAL fact_gamma,fact_epsilon,fact_gamma2
86	REAL c2(ngrid,klev)
87
88	REAL zw2fact,expa
89	Zsat=.false.
90	! Initialisation
91	RLvCp = RLVTT/RCPD
92
93	if (iflag_thermals_ed==0) then
94	fact_gamma=1.
95	fact_epsilon=1.
96	else if (iflag_thermals_ed==1) then
97	! Valeurs au moment de la premiere soumission des papiers
98	fact_gamma=1.
99	fact_epsilon=0.002
100	fact_gamma2=0.6
101	! Suggestions des Fleurs, Septembre 2009
102	fact_epsilon=0.015
103	!test cr
104	! fact_epsilon=0.002
105	fact_gamma=0.9
106	fact_gamma2=0.7
107
108	else if (iflag_thermals_ed==2) then
109	fact_gamma=1.
110	fact_epsilon=2.
111	else if (iflag_thermals_ed==3) then
112	fact_gamma=3./4.
113	fact_epsilon=3.
114	endif
115
116	! write(lunout,*)'THERM 31H '
117
118	! Initialisations des variables reeles
119	if (1==0) then
120	ztva(:,:)=ztv(:,:)
121	ztva_est(:,:)=ztva(:,:)
122	ztla(:,:)=zthl(:,:)
123	zqta(:,:)=po(:,:)
124	zha(:,:) = ztva(:,:)
125	else
126	ztva(:,:)=0.
127	ztva_est(:,:)=0.
128	ztla(:,:)=0.
129	zqta(:,:)=0.
130	zha(:,:) =0.
131	endif
132
133	zqla_est(:,:)=0.
134	zqsatth(:,:)=0.
135	zqla(:,:)=0.
136	detr_star(:,:)=0.
137	entr_star(:,:)=0.
138	alim_star(:,:)=0.
139	alim_star_tot(:)=0.
140	csc(:,:)=0.
141	detr(:,:)=0.
142	entr(:,:)=0.
143	zw2(:,:)=0.
144	w_est(:,:)=0.
145	f_star(:,:)=0.
146	wa_moy(:,:)=0.
147	linter(:)=1.
148	linter(:)=1.
149
150	! Initialisation des variables entieres
151	lmix(:)=1
152	lmix_bis(:)=2
153	wmaxa(:)=0.
154	lalim(:)=1
155
156	!-------------------------------------------------------------------------
157	! On ne considere comme actif que les colonnes dont les deux premieres
158	! couches sont instables.
159	!-------------------------------------------------------------------------
160	active(:)=ztv(:,1)>ztv(:,2)
161
162	!-------------------------------------------------------------------------
163	! Definition de l'alimentation a l'origine dans thermcell_init
164	!-------------------------------------------------------------------------
165	do l=1,klev-1
166	do ig=1,ngrid
167	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
168	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
169	& *sqrt(zlev(ig,l+1))
170	lalim(:)=l+1
171	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
172	endif
173	enddo
174	enddo
175	do l=1,klev
176	do ig=1,ngrid
177	if (alim_star_tot(ig) > 1.e-10 ) then
178	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
179	endif
180	enddo
181	enddo
182	alim_star_tot(:)=1.
183
184
185	!------------------------------------------------------------------------------
186	! Calcul dans la premiere couche
187	! On decide dans cette version que le thermique n'est actif que si la premiere
188	! couche est instable.
189	! Pourrait etre change si on veut que le thermiques puisse se déclencher
190	! dans une couche l>1
191	!------------------------------------------------------------------------------
192	do ig=1,ngrid
193	! Le panache va prendre au debut les caracteristiques de l'air contenu
194	! dans cette couche.
195	if (active(ig)) then
196	ztla(ig,1)=zthl(ig,1)
197	zqta(ig,1)=po(ig,1)
198	zqla(ig,1)=zl(ig,1)
199	!cr: attention, prise en compte de f*(1)=1
200	f_star(ig,2)=alim_star(ig,1)
201	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
202	& *(zlev(ig,2)-zlev(ig,1)) &
203	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
204	w_est(ig,2)=zw2(ig,2)
205	endif
206	enddo
207	!
208
209	!==============================================================================
210	!boucle de calcul de la vitesse verticale dans le thermique
211	!==============================================================================
212	do l=2,klev-1
213	!==============================================================================
214
215
216	! On decide si le thermique est encore actif ou non
217	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
218	do ig=1,ngrid
219	active(ig)=active(ig) &
220	& .and. zw2(ig,l)>1.e-10 &
221	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
222	enddo
223
224
225
226	! Premier calcul de la vitesse verticale a partir de la temperature
227	! potentielle virtuelle
228	! if (1.eq.1) then
229	! w_est(ig,3)=zw2(ig,2)* &
230	! & ((f_star(ig,2))**2) &
231	! & /(f_star(ig,2)+alim_star(ig,2))**2+ &
232	! & 2.RG(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
233	! & *(zlev(ig,3)-zlev(ig,2))
234	! endif
235
236
237	!---------------------------------------------------------------------------
238	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
239	! sans tenir compte du detrainement et de l'entrainement dans cette
240	! couche
241	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
242	! avant) a l'alimentation pour avoir un calcul plus propre
243	!---------------------------------------------------------------------------
244
245	call thermcell_condens(ngrid,active, &
246	& zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
247
248	do ig=1,ngrid
249	if(active(ig)) then
250	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
251	zta_est(ig,l)=ztva_est(ig,l)
252	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
253	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
254	& -zqla_est(ig,l))-zqla_est(ig,l))
255
256	w_est(ig,l+1)=zw2(ig,l)* &
257	& ((f_star(ig,l))**2) &
258	& /(f_star(ig,l)+alim_star(ig,l))**2+ &
259	& 2.RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
260	& *(zlev(ig,l+1)-zlev(ig,l))
261	if (w_est(ig,l+1).lt.0.) then
262	w_est(ig,l+1)=zw2(ig,l)
263	endif
264	endif
265	enddo
266
267	!-------------------------------------------------
268	!calcul des taux d'entrainement et de detrainement
269	!-------------------------------------------------
270
271	do ig=1,ngrid
272	if (active(ig)) then
273	zdz=zlev(ig,l+1)-zlev(ig,l)
274	zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
275	zfact=fact_gamma/(1.+fact_gamma)
276
277	! estimation de la fraction couverte par les thermiques
278	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
279
280	!calcul de la soumission papier
281	if (1.eq.1) then
282	fact_epsilon=0.0007
283	! zfact=0.9/(1.+0.9)
284	zfact=0.3
285	fact_gamma=0.7
286	fact_gamma2=0.6
287	expa=0.25
288	! Calcul du taux d'entrainement entr_star (epsilon)
289	entr_star(ig,l)=f_star(ig,l)zdz ( zfact * MAX(0., &
290	& zbuoy/w_est(ig,l+1) )&
291	!- fact_epsilon/zalpha**0.25 ) &
292	& +0.000 )
293
294	! entr_star(ig,l)=f_star(ig,l)zdz ( 1./3 * MAX(0., &
295	! & zbuoy/w_est(ig,l+1) - 1./zalpha**0.25 ) &
296	! & +0.000 )
297	! Calcul du taux de detrainement detr_star (delta)
298	! if (zqla_est(ig,l).lt.1.e-10) then
299	! detr_star(ig,l)=f_star(ig,l)zdz ( &
300	! & fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) ) &
301	! & +0.0006 )
302	! else
303	! detr_star(ig,l)=f_star(ig,l)zdz ( &
304	! & fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) ) &
305	! & +0.002 )
306	! endif
307	detr_star(ig,l)=f_star(ig,l)zdz ( &
308	& fact_gamma2 * MAX(0., &
309	!fact_epsilon/zalpha**0.25
310	& -zbuoy/w_est(ig,l+1) ) &
311	! & + 0.002(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))*0.6 &
312	!test
313	& + 0.006(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))*0.6 &
314	& +0.0000 )
315	else
316
317	! Calcul du taux d'entrainement entr_star (epsilon)
318	entr_star(ig,l)=f_star(ig,l)zdz ( zfact * MAX(0., &
319	& zbuoy/w_est(ig,l+1) - fact_epsilon ) &
320	& +0.0000 )
321
322	! Calcul du taux de detrainement detr_star (delta)
323	detr_star(ig,l)=f_star(ig,l)zdz ( &
324	& fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) ) &
325	& + 0.002(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))*0.6 &
326	& +0.0000 )
327
328	endif
329	!endif choix du calcul de E* et D*
330
331	!cr test
332	! entr_star(ig,l)=entr_star(ig,l)+0.2*detr_star(ig,l)
333
334	! Prise en compte de la fraction
335	! detr_star(ig,l)=detr_star(ig,l)*sqrt(0.01/max(zalpha,1.e-5))
336
337	! En dessous de lalim, on prend le max de alim_star et entr_star pour
338	! alim_star et 0 sinon
339	if (l.lt.lalim(ig)) then
340	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
341	entr_star(ig,l)=0.
342	endif
343
344	! Calcul du flux montant normalise
345	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
346	& -detr_star(ig,l)
347
348	endif
349	enddo
350
351	!----------------------------------------------------------------------------
352	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
353	!---------------------------------------------------------------------------
354	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
355	do ig=1,ngrid
356	if (activetmp(ig)) then
357	Zsat=.false.
358	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
359	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
360	& /(f_star(ig,l+1)+detr_star(ig,l))
361	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
362	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
363	& /(f_star(ig,l+1)+detr_star(ig,l))
364
365	endif
366	enddo
367
368	call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
369
370
371	do ig=1,ngrid
372	if (activetmp(ig)) then
373	if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
374	! on ecrit de maniere conservative (sat ou non)
375	! T = Tl +Lv/Cp ql
376	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
377	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
378	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
379	zha(ig,l) = ztva(ig,l)
380	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
381	& -zqla(ig,l))-zqla(ig,l))
382
383	!on ecrit zqsat
384	zqsatth(ig,l)=qsatbef
385
386	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
387	! zw2(ig,l+1)=&
388	! & zw2(ig,l)(1-fact_epsilon/(1.+fact_gamma)2.*(zlev(ig,l+1)-zlev(ig,l))) &
389	! & +2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
390	! & *1./(1.+fact_gamma) &
391	! & *(zlev(ig,l+1)-zlev(ig,l))
392	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
393	! La meme en plus modulaire :
394	zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
395	zdz=zlev(ig,l+1)-zlev(ig,l)
396
397
398	zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
399
400	if (1==0) then
401	zw2modif=zw2(ig,l)(1-fact_epsilon/(1.+fact_gamma)2.*zdz)
402	zdw2=2.zbuoy/(1.+fact_gamma)zdz
403	zw2(ig,l+1)=zw2modif+zdw2
404	else
405	! Tentative de reecriture de l'equation de w2. A reprendre ...
406	! zdw2=2zdzzbuoy
407	! zw2modif=zw2(ig,l)(1.-2.zdz*(zeps+fact_epsilon))
408	!!!!! zw2(ig,l+1)=(zw2(ig,l)+zdw2)/(1.+2.*zfactw2(ig,l))
409	!formulation Arnaud
410	! zw2fact=zbuoyzalpha*expa/fact_epsilon
411	! zw2(ig,l+1)=(zw2(ig,l)-zw2fact)exp(-2.fact_epsilon/(zalpha*expa(1+fact_gamma))*zdz) &
412	! & +zw2fact
413	if (zbuoy.gt.1.e-10) then
414	zw2fact=zbuoyzalphaexpa/fact_epsilon(fact_gamma-zfact)
415	zw2(ig,l+1)=(zw2(ig,l)-zw2fact)exp(-2.fact_epsilon/(zalpha*expa)zdz) &
416	& +zw2fact
417	else
418	zw2fact=zbuoyzalphaexpa/fact_epsilon(fact_gamma)
419	zw2(ig,l+1)=(zw2(ig,l)-zw2fact)exp(-2.fact_epsilon/(zalpha*expa)zdz) &
420	& +zw2fact
421
422	endif
423
424	endif
425	! zw2(ig,l+1)=zw2modif+zdw2
426	endif
427	enddo
428
429	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
430	!
431	!---------------------------------------------------------------------------
432	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
433	!---------------------------------------------------------------------------
434
435	do ig=1,ngrid
436	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
437	! stop'On tombe sur le cas particulier de thermcell_dry'
438	write(lunout,*) &
439	& 'On tombe sur le cas particulier de thermcell_plume'
440	zw2(ig,l+1)=0.
441	linter(ig)=l+1
442	endif
443
444	if (zw2(ig,l+1).lt.0.) then
445	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
446	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
447	zw2(ig,l+1)=0.
448	endif
449
450	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
451
452	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
453	! lmix est le niveau de la couche ou w (wa_moy) est maximum
454	!on rajoute le calcul de lmix_bis
455	if (zqla(ig,l).lt.1.e-10) then
456	lmix_bis(ig)=l+1
457	endif
458	lmix(ig)=l+1
459	wmaxa(ig)=wa_moy(ig,l+1)
460	endif
461	enddo
462
463	!=========================================================================
464	! FIN DE LA BOUCLE VERTICALE
465	enddo
466	!=========================================================================
467
468	!on recalcule alim_star_tot
469	do ig=1,ngrid
470	alim_star_tot(ig)=0.
471	enddo
472	do ig=1,ngrid
473	do l=1,lalim(ig)-1
474	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
475	enddo
476	enddo
477
478
479	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
480
481	! print*,'thermcell_plume OK'
482
483	return
484	end
485
486
487	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
488	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
489	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
490	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
491	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
492	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
493	SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
494	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
495	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
496	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
497	& ,lev_out,lunout1,igout)
498
499	!--------------------------------------------------------------------------
500	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
501	! Version conforme a l'article de Rio et al. 2010.
502	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
503	!--------------------------------------------------------------------------
504
505	IMPLICIT NONE
506
507	#include "YOMCST.h"
508	#include "YOETHF.h"
509	#include "FCTTRE.h"
510	#include "iniprint.h"
511	#include "thermcell.h"
512
513	INTEGER itap
514	INTEGER lunout1,igout
515	INTEGER ngrid,klev
516	REAL ptimestep
517	REAL ztv(ngrid,klev)
518	REAL zthl(ngrid,klev)
519	REAL po(ngrid,klev)
520	REAL zl(ngrid,klev)
521	REAL rhobarz(ngrid,klev)
522	REAL zlev(ngrid,klev+1)
523	REAL pplev(ngrid,klev+1)
524	REAL pphi(ngrid,klev)
525	REAL zpspsk(ngrid,klev)
526	REAL alim_star(ngrid,klev)
527	REAL f0(ngrid)
528	INTEGER lalim(ngrid)
529	integer lev_out ! niveau pour les print
530
531	real alim_star_tot(ngrid)
532
533	REAL ztva(ngrid,klev)
534	REAL ztla(ngrid,klev)
535	REAL zqla(ngrid,klev)
536	REAL zqta(ngrid,klev)
537	REAL zha(ngrid,klev)
538
539	REAL detr_star(ngrid,klev)
540	REAL coefc
541	REAL entr_star(ngrid,klev)
542	REAL detr(ngrid,klev)
543	REAL entr(ngrid,klev)
544
545	REAL csc(ngrid,klev)
546
547	REAL zw2(ngrid,klev+1)
548	REAL w_est(ngrid,klev+1)
549	REAL f_star(ngrid,klev+1)
550	REAL wa_moy(ngrid,klev+1)
551
552	REAL ztva_est(ngrid,klev)
553	REAL zqla_est(ngrid,klev)
554	REAL zqsatth(ngrid,klev)
555	REAL zta_est(ngrid,klev)
556	REAL zdw2
557	REAL zw2modif
558	REAL zw2fact
559	REAL zeps(ngrid,klev)
560
561	REAL linter(ngrid)
562	INTEGER lmix(ngrid)
563	INTEGER lmix_bis(ngrid)
564	REAL wmaxa(ngrid)
565
566	INTEGER ig,l,k
567
568	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
569	real zbuoybis
570	real zcor,zdelta,zcvm5,qlbef,zdz2
571	real betalpha,zbetalpha
572	real Tbef,qsatbef,b1,eps, afact
573	real dqsat_dT,DT,num,denom,m
574	REAL REPS,RLvCp,DDT0
575	PARAMETER (DDT0=.01)
576	logical Zsat
577	LOGICAL active(ngrid),activetmp(ngrid)
578	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
579	REAL c2(ngrid,klev)
580	Zsat=.false.
581	! Initialisation
582
583	RLvCp = RLVTT/RCPD
584	fact_epsilon=0.002
585	betalpha=0.9
586	afact=2./3.
587
588	zbetalpha=betalpha/(1.+betalpha)
589
590	! print*,'THERM 31B'
591
592	! Initialisations des variables reeles
593	if (1==0) then
594	ztva(:,:)=ztv(:,:)
595	ztva_est(:,:)=ztva(:,:)
596	ztla(:,:)=zthl(:,:)
597	zqta(:,:)=po(:,:)
598	zha(:,:) = ztva(:,:)
599	else
600	ztva(:,:)=0.
601	ztva_est(:,:)=0.
602	ztla(:,:)=0.
603	zqta(:,:)=0.
604	zha(:,:) =0.
605	endif
606
607	zqla_est(:,:)=0.
608	zqsatth(:,:)=0.
609	zqla(:,:)=0.
610	detr_star(:,:)=0.
611	entr_star(:,:)=0.
612	alim_star(:,:)=0.
613	alim_star_tot(:)=0.
614	csc(:,:)=0.
615	detr(:,:)=0.
616	entr(:,:)=0.
617	zw2(:,:)=0.
618	zbuoy(:,:)=0.
619	gamma(:,:)=0.
620	zeps(:,:)=0.
621	w_est(:,:)=0.
622	f_star(:,:)=0.
623	wa_moy(:,:)=0.
624	linter(:)=1.
625	! linter(:)=1.
626	b1=2.
627	! Initialisation des variables entieres
628	lmix(:)=1
629	lmix_bis(:)=2
630	wmaxa(:)=0.
631	lalim(:)=1
632
633	! print*,'THERMCELL PLUME ARNAUD DEDANS'
634
635	!-------------------------------------------------------------------------
636	! On ne considere comme actif que les colonnes dont les deux premieres
637	! couches sont instables.
638	!-------------------------------------------------------------------------
639	active(:)=ztv(:,1)>ztv(:,2)
640
641	!-------------------------------------------------------------------------
642	! Definition de l'alimentation a l'origine dans thermcell_init
643	!-------------------------------------------------------------------------
644	do l=1,klev-1
645	do ig=1,ngrid
646	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
647	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
648	& *sqrt(zlev(ig,l+1))
649	lalim(:)=l+1
650	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
651	endif
652	enddo
653	enddo
654	do l=1,klev
655	do ig=1,ngrid
656	if (alim_star_tot(ig) > 1.e-10 ) then
657	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
658	endif
659	enddo
660	enddo
661	alim_star_tot(:)=1.
662
663
664
665	!------------------------------------------------------------------------------
666	! Calcul dans la premiere couche
667	! On decide dans cette version que le thermique n'est actif que si la premiere
668	! couche est instable.
669	! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
670	! dans une couche l>1
671	!------------------------------------------------------------------------------
672	do ig=1,ngrid
673	! Le panache va prendre au debut les caracteristiques de l'air contenu
674	! dans cette couche.
675	if (active(ig)) then
676	ztla(ig,1)=zthl(ig,1)
677	zqta(ig,1)=po(ig,1)
678	zqla(ig,1)=zl(ig,1)
679	!cr: attention, prise en compte de f*(1)=1
680	f_star(ig,2)=alim_star(ig,1)
681	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
682	& *(zlev(ig,2)-zlev(ig,1)) &
683	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
684	w_est(ig,2)=zw2(ig,2)
685	endif
686	enddo
687	!
688
689	!==============================================================================
690	!boucle de calcul de la vitesse verticale dans le thermique
691	!==============================================================================
692	do l=2,klev-1
693	!==============================================================================
694
695
696	! On decide si le thermique est encore actif ou non
697	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
698	do ig=1,ngrid
699	active(ig)=active(ig) &
700	& .and. zw2(ig,l)>1.e-10 &
701	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
702	enddo
703
704
705
706	!---------------------------------------------------------------------------
707	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
708	! sans tenir compte du detrainement et de l'entrainement dans cette
709	! couche
710	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
711	! avant) a l'alimentation pour avoir un calcul plus propre
712	!---------------------------------------------------------------------------
713
714	call thermcell_condens(ngrid,active, &
715	& zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
716
717
718	do ig=1,ngrid
719	! print*,'active',active(ig),ig,l
720	if(active(ig)) then
721	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
722	zta_est(ig,l)=ztva_est(ig,l)
723	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
724	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
725	& -zqla_est(ig,l))-zqla_est(ig,l))
726
727	!------------------------------------------------
728	!AJAM:nouveau calcul de w²
729	!------------------------------------------------
730	zdz=zlev(ig,l+1)-zlev(ig,l)
731	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
732
733	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
734	zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
735	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
736
737
738	if (w_est(ig,l+1).lt.0.) then
739	w_est(ig,l+1)=zw2(ig,l)
740	endif
741	endif
742	enddo
743
744
745	!-------------------------------------------------
746	!calcul des taux d'entrainement et de detrainement
747	!-------------------------------------------------
748
749	do ig=1,ngrid
750	if (active(ig)) then
751
752	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
753	zw2m=w_est(ig,l+1)
754	zdz=zlev(ig,l+1)-zlev(ig,l)
755	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
756	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
757	zbuoybis=zbuoy(ig,l)
758	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
759	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
760
761
762	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
763	& afact*zbuoybis/zw2m - fact_epsilon )
764
765
766	detr_star(ig,l)=f_star(ig,l)*zdz &
767	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
768	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
769
770	! En dessous de lalim, on prend le max de alim_star et entr_star pour
771	! alim_star et 0 sinon
772	if (l.lt.lalim(ig)) then
773	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
774	entr_star(ig,l)=0.
775	endif
776
777	! Calcul du flux montant normalise
778	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
779	& -detr_star(ig,l)
780
781	endif
782	enddo
783
784
785	!----------------------------------------------------------------------------
786	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
787	!---------------------------------------------------------------------------
788	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
789	do ig=1,ngrid
790	if (activetmp(ig)) then
791	Zsat=.false.
792	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
793	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
794	& /(f_star(ig,l+1)+detr_star(ig,l))
795	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
796	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
797	& /(f_star(ig,l+1)+detr_star(ig,l))
798
799	endif
800	enddo
801
802	call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
803
804
805	do ig=1,ngrid
806	if (activetmp(ig)) then
807	if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
808	! on ecrit de maniere conservative (sat ou non)
809	! T = Tl +Lv/Cp ql
810	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
811	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
812	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
813	zha(ig,l) = ztva(ig,l)
814	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
815	& -zqla(ig,l))-zqla(ig,l))
816
817	!on ecrit zqsat
818	zqsatth(ig,l)=qsatbef
819
820	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
821	zdz=zlev(ig,l+1)-zlev(ig,l)
822	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
823
824	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
825	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
826	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
827	endif
828	enddo
829
830	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
831	!
832	!---------------------------------------------------------------------------
833	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
834	!---------------------------------------------------------------------------
835
836	do ig=1,ngrid
837	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
838	! stop'On tombe sur le cas particulier de thermcell_dry'
839	print*,'On tombe sur le cas particulier de thermcell_plume'
840	zw2(ig,l+1)=0.
841	linter(ig)=l+1
842	endif
843
844	if (zw2(ig,l+1).lt.0.) then
845	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
846	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
847	zw2(ig,l+1)=0.
848	endif
849
850	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
851
852	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
853	! lmix est le niveau de la couche ou w (wa_moy) est maximum
854	!on rajoute le calcul de lmix_bis
855	if (zqla(ig,l).lt.1.e-10) then
856	lmix_bis(ig)=l+1
857	endif
858	lmix(ig)=l+1
859	wmaxa(ig)=wa_moy(ig,l+1)
860	endif
861	enddo
862
863	!=========================================================================
864	! FIN DE LA BOUCLE VERTICALE
865	enddo
866	!=========================================================================
867
868	! print*,'THERMCELL PLUME ARNAUD DEDANS 7'
869	!on recalcule alim_star_tot
870	do ig=1,ngrid
871	alim_star_tot(ig)=0.
872	enddo
873	do ig=1,ngrid
874	do l=1,lalim(ig)-1
875	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
876	enddo
877	enddo
878
879
880	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
881
882	#undef wrgrads_thermcell
883	#ifdef wrgrads_thermcell
884	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
885	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
886	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
887	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
888	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
889	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
890	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
891	#endif
892
893
894	! print*,'THERMCELL PLUME ARNAUD DEDANS 8'
895	return
896	end
897

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

Download in other formats:

Original Format