Context Navigation

thermcell_plume.F90 @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 27.3 KB

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

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

Download in other formats:

Original Format