Context Navigation

thermcell_plume.F90 @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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: 27.4 KB

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

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

Download in other formats:

Original Format