Context Navigation

thermcell_plume.F90 @ 409

Last change on this file since 409 was 179, checked in by lfita, 10 years ago
No Ale, Alp, ... values because pphi is at the intermediate level!!! Removing raw prints of checkig for Ale, Alp
File size: 27.7 KB

Line
1	!
2	! $Id: thermcell_plume.F90 1503 2011-03-23 11:57:52Z idelkadi $
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	integer nbpb
40	real zcon2(ngrid)
41
42	real alim_star_tot(ngrid)
43
44	REAL ztva(ngrid,klev)
45	REAL ztla(ngrid,klev)
46	REAL zqla(ngrid,klev)
47	REAL zqta(ngrid,klev)
48	REAL zha(ngrid,klev)
49
50	REAL detr_star(ngrid,klev)
51	REAL coefc
52	REAL entr_star(ngrid,klev)
53	REAL detr(ngrid,klev)
54	REAL entr(ngrid,klev)
55
56	REAL csc(ngrid,klev)
57
58	REAL zw2(ngrid,klev+1)
59	REAL w_est(ngrid,klev+1)
60	REAL f_star(ngrid,klev+1)
61	REAL wa_moy(ngrid,klev+1)
62
63	REAL ztva_est(ngrid,klev)
64	REAL zqla_est(ngrid,klev)
65	REAL zqsatth(ngrid,klev)
66	REAL zta_est(ngrid,klev)
67	REAL zdw2
68	REAL zw2modif
69	REAL zeps
70
71	REAL linter(ngrid)
72	INTEGER lmix(ngrid)
73	INTEGER lmix_bis(ngrid)
74	REAL wmaxa(ngrid)
75
76	INTEGER ig,l,k
77
78	real zdz,zfact,zbuoy,zalpha,zdrag
79	real zcor,zdelta,zcvm5,qlbef
80	real Tbef,qsatbef
81	real dqsat_dT,DT,num,denom
82	REAL REPS,RLvCp,DDT0
83	PARAMETER (DDT0=.01)
84	logical Zsat
85	LOGICAL active(ngrid),activetmp(ngrid)
86	REAL fact_gamma,fact_epsilon,fact_gamma2
87	REAL c2(ngrid,klev)
88	REAL a1,m
89
90	REAL zw2fact,expa
91	Zsat=.false.
92
93	! Initialisation
94	RLvCp = RLVTT/RCPD
95
96
97	fact_epsilon=0.002
98	a1=2./3.
99	fact_gamma=0.9
100	zfact=fact_gamma/(1+fact_gamma)
101	fact_gamma2=zfact
102	expa=0.
103
104
105	! Initialisations des variables reeles
106	if (1==1) then
107	ztva(:,:)=ztv(:,:)
108	ztva_est(:,:)=ztva(:,:)
109	ztla(:,:)=zthl(:,:)
110	zqta(:,:)=po(:,:)
111	zha(:,:) = ztva(:,:)
112	else
113	ztva(:,:)=0.
114	ztva_est(:,:)=0.
115	ztla(:,:)=0.
116	zqta(:,:)=0.
117	zha(:,:) =0.
118	endif
119
120	zqla_est(:,:)=0.
121	zqsatth(:,:)=0.
122	zqla(:,:)=0.
123	detr_star(:,:)=0.
124	entr_star(:,:)=0.
125	alim_star(:,:)=0.
126	alim_star_tot(:)=0.
127	csc(:,:)=0.
128	detr(:,:)=0.
129	entr(:,:)=0.
130	zw2(:,:)=0.
131	w_est(:,:)=0.
132	f_star(:,:)=0.
133	wa_moy(:,:)=0.
134	linter(:)=1.
135	linter(:)=1.
136
137	! Initialisation des variables entieres
138	lmix(:)=1
139	lmix_bis(:)=2
140	wmaxa(:)=0.
141	lalim(:)=1
142
143	!-------------------------------------------------------------------------
144	! On ne considere comme actif que les colonnes dont les deux premieres
145	! couches sont instables.
146	!-------------------------------------------------------------------------
147	active(:)=ztv(:,1)>ztv(:,2)
148
149	!-------------------------------------------------------------------------
150	! Definition de l'alimentation a l'origine dans thermcell_init
151	!-------------------------------------------------------------------------
152	do l=1,klev-1
153	do ig=1,ngrid
154	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
155	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
156	& *sqrt(zlev(ig,l+1))
157	lalim(ig)=l+1
158	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
159	endif
160	enddo
161	enddo
162	do l=1,klev
163	do ig=1,ngrid
164	if (alim_star_tot(ig) > 1.e-10 ) then
165	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
166	endif
167	enddo
168	enddo
169	alim_star_tot(:)=1.
170
171	!------------------------------------------------------------------------------
172	! Calcul dans la premiere couche
173	! On decide dans cette version que le thermique n'est actif que si la premiere
174	! couche est instable.
175	! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
176	! dans une couche l>1
177	!------------------------------------------------------------------------------
178	do ig=1,ngrid
179	! Le panache va prendre au debut les caracteristiques de l'air contenu
180	! dans cette couche.
181	if (active(ig)) then
182	ztla(ig,1)=zthl(ig,1)
183	zqta(ig,1)=po(ig,1)
184	zqla(ig,1)=zl(ig,1)
185	!cr: attention, prise en compte de f*(1)=1
186	f_star(ig,2)=alim_star(ig,1)
187	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
188	& *(zlev(ig,2)-zlev(ig,1)) &
189	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
190	w_est(ig,2)=zw2(ig,2)
191	endif
192	enddo
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	! on ecrit de maniere conservative (sat ou non)
346	! T = Tl +Lv/Cp ql
347	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
348	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
349	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
350	zha(ig,l) = ztva(ig,l)
351	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
352	& -zqla(ig,l))-zqla(ig,l))
353
354	!on ecrit zqsat
355	zqsatth(ig,l)=qsatbef
356
357	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
358	! zw2(ig,l+1)=&
359	! & zw2(ig,l)(1-fact_epsilon/(1.+fact_gamma)2.*(zlev(ig,l+1)-zlev(ig,l))) &
360	! & +2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
361	! & *1./(1.+fact_gamma) &
362	! & *(zlev(ig,l+1)-zlev(ig,l))
363	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
364	! La meme en plus modulaire :
365	zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
366	zdz=zlev(ig,l+1)-zlev(ig,l)
367
368
369	zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
370
371	!if (1==0) then
372	! zw2modif=zw2(ig,l)(1-fact_epsilon/(1.+fact_gamma)2.*zdz)
373	! zdw2=2.zbuoy/(1.+fact_gamma)zdz
374	! zw2(ig,l+1)=zw2modif+zdw2
375	!else
376	zdrag=fact_epsilon/(zalpha**expa)
377	zw2fact=zbuoy/zdrag*a1
378	zw2(ig,l+1)=(zw2(ig,l)-zw2fact)exp(-2.zdrag/(1+fact_gamma)*zdz) &
379	& +zw2fact
380	!endif
381
382	endif
383	enddo
384
385	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
386	!
387	!---------------------------------------------------------------------------
388	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
389	!---------------------------------------------------------------------------
390
391	nbpb=0
392	do ig=1,ngrid
393	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
394	! stop'On tombe sur le cas particulier de thermcell_dry'
395	! print*,'On tombe sur le cas particulier de thermcell_plume'
396	nbpb=nbpb+1
397	zw2(ig,l+1)=0.
398	linter(ig)=l+1
399	endif
400
401	if (zw2(ig,l+1).lt.0.) then
402	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
403	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
404	zw2(ig,l+1)=0.
405	endif
406
407	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
408
409	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
410	! lmix est le niveau de la couche ou w (wa_moy) est maximum
411	!on rajoute le calcul de lmix_bis
412	if (zqla(ig,l).lt.1.e-10) then
413	lmix_bis(ig)=l+1
414	endif
415	lmix(ig)=l+1
416	wmaxa(ig)=wa_moy(ig,l+1)
417	endif
418	enddo
419
420	if (nbpb>0) then
421	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
422	endif
423
424	!=========================================================================
425	! FIN DE LA BOUCLE VERTICALE
426	enddo
427	!=========================================================================
428
429	!on recalcule alim_star_tot
430	do ig=1,ngrid
431	alim_star_tot(ig)=0.
432	enddo
433	do ig=1,ngrid
434	do l=1,lalim(ig)-1
435	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
436	enddo
437	enddo
438
439
440	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
441
442
443	return
444	end
445
446	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
447	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
448	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
449	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
450	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
451	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
452	SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
453	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
454	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
455	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
456	& ,lev_out,lunout1,igout)
457
458	!--------------------------------------------------------------------------
459	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
460	! Version conforme a l'article de Rio et al. 2010.
461	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
462	!--------------------------------------------------------------------------
463
464	IMPLICIT NONE
465
466	#include "YOMCST.h"
467	#include "YOETHF.h"
468	#include "FCTTRE.h"
469	#include "iniprint.h"
470	#include "thermcell.h"
471
472	INTEGER itap
473	INTEGER lunout1,igout
474	INTEGER ngrid,klev
475	REAL ptimestep
476	REAL ztv(ngrid,klev)
477	REAL zthl(ngrid,klev)
478	REAL po(ngrid,klev)
479	REAL zl(ngrid,klev)
480	REAL rhobarz(ngrid,klev)
481	REAL zlev(ngrid,klev+1)
482	REAL pplev(ngrid,klev+1)
483	REAL pphi(ngrid,klev)
484	REAL zpspsk(ngrid,klev)
485	REAL alim_star(ngrid,klev)
486	REAL f0(ngrid)
487	INTEGER lalim(ngrid)
488	integer lev_out ! niveau pour les print
489	integer nbpb
490
491	real alim_star_tot(ngrid)
492
493	REAL ztva(ngrid,klev)
494	REAL ztla(ngrid,klev)
495	REAL zqla(ngrid,klev)
496	REAL zqta(ngrid,klev)
497	REAL zha(ngrid,klev)
498
499	REAL detr_star(ngrid,klev)
500	REAL coefc
501	REAL entr_star(ngrid,klev)
502	REAL detr(ngrid,klev)
503	REAL entr(ngrid,klev)
504
505	REAL csc(ngrid,klev)
506
507	REAL zw2(ngrid,klev+1)
508	REAL w_est(ngrid,klev+1)
509	REAL f_star(ngrid,klev+1)
510	REAL wa_moy(ngrid,klev+1)
511
512	REAL ztva_est(ngrid,klev)
513	REAL zqla_est(ngrid,klev)
514	REAL zqsatth(ngrid,klev)
515	REAL zta_est(ngrid,klev)
516	REAL ztemp(ngrid),zqsat(ngrid)
517	REAL zdw2
518	REAL zw2modif
519	REAL zw2fact
520	REAL zeps(ngrid,klev)
521
522	REAL linter(ngrid)
523	INTEGER lmix(ngrid)
524	INTEGER lmix_bis(ngrid)
525	REAL wmaxa(ngrid)
526
527	INTEGER ig,l,k
528
529	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
530	real zbuoybis
531	real zcor,zdelta,zcvm5,qlbef,zdz2
532	real betalpha,zbetalpha
533	real eps, afact
534	REAL REPS,RLvCp,DDT0
535	PARAMETER (DDT0=.01)
536	logical Zsat
537	LOGICAL active(ngrid),activetmp(ngrid)
538	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
539	REAL c2(ngrid,klev)
540
541	Zsat=.false.
542	! Initialisation
543
544	RLvCp = RLVTT/RCPD
545	fact_epsilon=0.002
546	betalpha=0.9
547	afact=2./3.
548
549	zbetalpha=betalpha/(1.+betalpha)
550
551
552	! Initialisations des variables reeles
553	if (1==0) then
554	ztva(:,:)=ztv(:,:)
555	ztva_est(:,:)=ztva(:,:)
556	ztla(:,:)=zthl(:,:)
557	zqta(:,:)=po(:,:)
558	zha(:,:) = ztva(:,:)
559	else
560	ztva(:,:)=0.
561	ztva_est(:,:)=0.
562	ztla(:,:)=0.
563	zqta(:,:)=0.
564	zha(:,:) =0.
565	endif
566
567	zqla_est(:,:)=0.
568	zqsatth(:,:)=0.
569	zqla(:,:)=0.
570	detr_star(:,:)=0.
571	entr_star(:,:)=0.
572	alim_star(:,:)=0.
573	alim_star_tot(:)=0.
574	csc(:,:)=0.
575	detr(:,:)=0.
576	entr(:,:)=0.
577	zw2(:,:)=0.
578	zbuoy(:,:)=0.
579	gamma(:,:)=0.
580	zeps(:,:)=0.
581	w_est(:,:)=0.
582	f_star(:,:)=0.
583	wa_moy(:,:)=0.
584	linter(:)=1.
585	! linter(:)=1.
586	! Initialisation des variables entieres
587	lmix(:)=1
588	lmix_bis(:)=2
589	wmaxa(:)=0.
590	lalim(:)=1
591
592
593	!-------------------------------------------------------------------------
594	! On ne considere comme actif que les colonnes dont les deux premieres
595	! couches sont instables.
596	!-------------------------------------------------------------------------
597	active(:)=ztv(:,1)>ztv(:,2)
598
599	!-------------------------------------------------------------------------
600	! Definition de l'alimentation a l'origine dans thermcell_init
601	!-------------------------------------------------------------------------
602	do l=1,klev-1
603	do ig=1,ngrid
604	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
605	! L. Fita, LMD July 2014. Adding IF for zlev == 0. NOT NEEDED!!
606	!! IF ( zlev(ig,l+1) < 0.) THEN
607	!! alim_star(ig,l) = 0.
608	!! ELSE
609	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
610	& *sqrt(zlev(ig,l+1))
611	lalim(ig)=l+1
612	!! END IF
613	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
614	endif
615	enddo
616	enddo
617	do l=1,klev
618	do ig=1,ngrid
619	if (alim_star_tot(ig) > 1.e-10 ) then
620	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
621	endif
622	enddo
623	enddo
624	alim_star_tot(:)=1.
625	!------------------------------------------------------------------------------
626	! Calcul dans la premiere couche
627	! On decide dans cette version que le thermique n'est actif que si la premiere
628	! couche est instable.
629	! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
630	! dans une couche l>1
631	!------------------------------------------------------------------------------
632	do ig=1,ngrid
633	! Le panache va prendre au debut les caracteristiques de l'air contenu
634	! dans cette couche.
635	if (active(ig)) then
636	ztla(ig,1)=zthl(ig,1)
637	zqta(ig,1)=po(ig,1)
638	zqla(ig,1)=zl(ig,1)
639	!cr: attention, prise en compte de f*(1)=1
640	f_star(ig,2)=alim_star(ig,1)
641	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
642	& *(zlev(ig,2)-zlev(ig,1)) &
643	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
644	w_est(ig,2)=zw2(ig,2)
645	endif
646	enddo
647	!
648
649	!==============================================================================
650	!boucle de calcul de la vitesse verticale dans le thermique
651	!==============================================================================
652	do l=2,klev-1
653	!==============================================================================
654
655
656	! On decide si le thermique est encore actif ou non
657	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
658	do ig=1,ngrid
659	active(ig)=active(ig) &
660	& .and. zw2(ig,l)>1.e-10 &
661	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
662	enddo
663
664
665
666	!---------------------------------------------------------------------------
667	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
668	! sans tenir compte du detrainement et de l'entrainement dans cette
669	! couche
670	! C'est a dire qu'on suppose
671	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
672	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
673	! avant) a l'alimentation pour avoir un calcul plus propre
674	!---------------------------------------------------------------------------
675
676	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
677	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
678
679	do ig=1,ngrid
680	! print*,'active',active(ig),ig,l
681	if(active(ig)) then
682	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
683	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
684	zta_est(ig,l)=ztva_est(ig,l)
685	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
686	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
687	& -zqla_est(ig,l))-zqla_est(ig,l))
688
689	!------------------------------------------------
690	!AJAM:nouveau calcul de w\B2
691	!------------------------------------------------
692	zdz=zlev(ig,l+1)-zlev(ig,l)
693	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
694
695	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
696	zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
697	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
698
699
700	if (w_est(ig,l+1).lt.0.) then
701	w_est(ig,l+1)=zw2(ig,l)
702	endif
703	endif
704	enddo
705
706
707	!-------------------------------------------------
708	!calcul des taux d'entrainement et de detrainement
709	!-------------------------------------------------
710
711	do ig=1,ngrid
712	if (active(ig)) then
713
714	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
715	zw2m=w_est(ig,l+1)
716	zdz=zlev(ig,l+1)-zlev(ig,l)
717	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
718	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
719	zbuoybis=zbuoy(ig,l)
720	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
721	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
722
723
724	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
725	& afact*zbuoybis/zw2m - fact_epsilon )
726
727
728	detr_star(ig,l)=f_star(ig,l)*zdz &
729	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
730	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
731
732	! En dessous de lalim, on prend le max de alim_star et entr_star pour
733	! alim_star et 0 sinon
734	if (l.lt.lalim(ig)) then
735	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
736	entr_star(ig,l)=0.
737	endif
738
739	! Calcul du flux montant normalise
740	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
741	& -detr_star(ig,l)
742
743	endif
744	enddo
745
746	!----------------------------------------------------------------------------
747	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
748	!---------------------------------------------------------------------------
749	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
750	do ig=1,ngrid
751	if (activetmp(ig)) then
752	Zsat=.false.
753	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
754	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
755	& /(f_star(ig,l+1)+detr_star(ig,l))
756	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
757	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
758	& /(f_star(ig,l+1)+detr_star(ig,l))
759
760	endif
761	enddo
762
763	ztemp(:)=zpspsk(:,l)*ztla(:,l)
764	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
765
766	do ig=1,ngrid
767	if (activetmp(ig)) then
768	! on ecrit de maniere conservative (sat ou non)
769	! T = Tl +Lv/Cp ql
770	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
771	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
772	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
773	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
774	zha(ig,l) = ztva(ig,l)
775	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
776	& -zqla(ig,l))-zqla(ig,l))
777	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
778	zdz=zlev(ig,l+1)-zlev(ig,l)
779	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
780
781	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
782	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
783	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
784	endif
785	enddo
786
787	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
788	!
789	!---------------------------------------------------------------------------
790	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
791	!---------------------------------------------------------------------------
792
793	nbpb=0
794	do ig=1,ngrid
795	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
796	! stop'On tombe sur le cas particulier de thermcell_dry'
797	! print*,'On tombe sur le cas particulier de thermcell_plume'
798	nbpb=nbpb+1
799	zw2(ig,l+1)=0.
800	linter(ig)=l+1
801	endif
802
803	if (zw2(ig,l+1).lt.0.) then
804	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
805	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
806	zw2(ig,l+1)=0.
807	endif
808
809	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
810
811	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
812	! lmix est le niveau de la couche ou w (wa_moy) est maximum
813	!on rajoute le calcul de lmix_bis
814	if (zqla(ig,l).lt.1.e-10) then
815	lmix_bis(ig)=l+1
816	endif
817	lmix(ig)=l+1
818	wmaxa(ig)=wa_moy(ig,l+1)
819	endif
820	enddo
821
822	if (nbpb>0) then
823	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
824	endif
825
826	!=========================================================================
827	! FIN DE LA BOUCLE VERTICALE
828	enddo
829	!=========================================================================
830
831	!on recalcule alim_star_tot
832	do ig=1,ngrid
833	alim_star_tot(ig)=0.
834	enddo
835	do ig=1,ngrid
836	do l=1,lalim(ig)-1
837	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
838	enddo
839	enddo
840
841
842	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
843
844	#undef wrgrads_thermcell
845	#ifdef wrgrads_thermcell
846	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
847	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
848	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
849	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
850	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
851	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
852	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
853	#endif
854
855
856	return
857	end
858

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

Download in other formats:

Original Format