Context Navigation

thermcell_flux2.F90 @ 5373

Last change on this file since 5373 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 14.9 KB

Line
1	SUBROUTINE thermcell_flux2(ngrid,klev,ptimestep,masse, &
2	& lalim,lmax,alim_star, &
3	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
4	& detr,zqla,lev_out,lunout1,igout)
5	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
6
7
8	!---------------------------------------------------------------------------
9	!thermcell_flux: deduction des flux
10	!---------------------------------------------------------------------------
11
12	IMPLICIT NONE
13	#include "iniprint.h"
14	#include "thermcell.h"
15
16	INTEGER ig,l
17	INTEGER ngrid,klev
18
19	REAL alim_star(ngrid,klev),entr_star(ngrid,klev)
20	REAL detr_star(ngrid,klev)
21	REAL zw2(ngrid,klev+1)
22	REAL zlev(ngrid,klev+1)
23	REAL masse(ngrid,klev)
24	REAL ptimestep
25	REAL rhobarz(ngrid,klev)
26	REAL f(ngrid)
27	INTEGER lmax(ngrid)
28	INTEGER lalim(ngrid)
29	REAL zqla(ngrid,klev)
30	REAL zmax(ngrid)
31
32	integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
33	integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
34
35
36	REAL entr(ngrid,klev),detr(ngrid,klev)
37	REAL fm(ngrid,klev+1)
38	REAL zfm
39
40	integer igout
41	integer lev_out
42	integer lunout1
43
44	REAL f_old,ddd0,eee0,ddd,eee,zzz
45
46	REAL fomass_max,alphamax
47	save fomass_max,alphamax
48
49	fomass_max=0.5
50	alphamax=0.7
51
52	ncorecfm1=0
53	ncorecfm2=0
54	ncorecfm3=0
55	ncorecfm4=0
56	ncorecfm5=0
57	ncorecfm6=0
58	ncorecfm7=0
59	ncorecfm8=0
60	ncorecalpha=0
61
62	!initialisation
63	fm(:,:)=0.
64
65	if (prt_level.ge.10) then
66	write(lunout1,*) 'Dans thermcell_flux 0'
67	write(lunout1,*) 'flux base ',f(igout)
68	write(lunout1,*) 'lmax ',lmax(igout)
69	write(lunout1,*) 'lalim ',lalim(igout)
70	write(lunout1,*) 'ig= ',igout
71	write(lunout1,) ' l E A* D* '
72	write(lunout1,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) &
73	& ,l=1,lmax(igout))
74	endif
75
76
77	!-------------------------------------------------------------------------
78	! Verification de la nullite des entrainement et detrainement au dessus
79	! de lmax(ig)
80	!-------------------------------------------------------------------------
81
82	do l=1,klev
83	do ig=1,ngrid
84	if (l.le.lmax(ig)) then
85	if (entr_star(ig,l).gt.1.) then
86	print*,'WARNING thermcell_flux 1 ig,l,lmax(ig)',ig,l,lmax(ig)
87	print*,'entr_star(ig,l)',entr_star(ig,l)
88	print*,'alim_star(ig,l)',alim_star(ig,l)
89	print*,'detr_star(ig,l)',detr_star(ig,l)
90	! stop
91	endif
92	else
93	if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0.) then
94	print*,'cas 1 : ig,l,lmax(ig)',ig,l,lmax(ig)
95	print*,'entr_star(ig,l)',entr_star(ig,l)
96	print*,'alim_star(ig,l)',alim_star(ig,l)
97	print*,'detr_star(ig,l)',detr_star(ig,l)
98	stop
99	endif
100	endif
101	enddo
102	enddo
103
104	!-------------------------------------------------------------------------
105	! Multiplication par le flux de masse issu de la femreture
106	!-------------------------------------------------------------------------
107
108	do l=1,klev
109	entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
110	detr(:,l)=f(:)*detr_star(:,l)
111	enddo
112
113	if (prt_level.ge.10) then
114	write(lunout1,*) 'Dans thermcell_flux 1'
115	write(lunout1,*) 'flux base ',f(igout)
116	write(lunout1,*) 'lmax ',lmax(igout)
117	write(lunout1,*) 'lalim ',lalim(igout)
118	write(lunout1,*) 'ig= ',igout
119	write(lunout1,*) ' l E D W2'
120	write(lunout1,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) &
121	& ,zw2(igout,l+1),l=1,lmax(igout))
122	endif
123
124	fm(:,1)=0.
125	do l=1,klev
126	do ig=1,ngrid
127	if (l.lt.lmax(ig)) then
128	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
129	elseif(l.eq.lmax(ig)) then
130	fm(ig,l+1)=0.
131	detr(ig,l)=fm(ig,l)+entr(ig,l)
132	else
133	fm(ig,l+1)=0.
134	endif
135	enddo
136	enddo
137
138
139
140	! Test provisoire : pour comprendre pourquoi on corrige plein de fois
141	! le cas fm6, on commence par regarder une premiere fois avant les
142	! autres corrections.
143
144	do l=1,klev
145	do ig=1,ngrid
146	if (detr(ig,l).gt.fm(ig,l)) then
147	ncorecfm8=ncorecfm8+1
148	! igout=ig
149	endif
150	enddo
151	enddo
152
153	! if (prt_level.ge.10) &
154	! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
155	! & ptimestep,masse,entr,detr,fm,'2 ')
156
157
158
159	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
160	! FH Version en cours de test;
161	! par rapport a thermcell_flux, on fait une grande boucle sur "l"
162	! et on modifie le flux avec tous les contr�les appliques d'affilee
163	! pour la meme couche
164	! Momentanement, on duplique le calcule du flux pour pouvoir comparer
165	! les flux avant et apres modif
166	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
167
168	do l=1,klev
169
170	do ig=1,ngrid
171	if (l.lt.lmax(ig)) then
172	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
173	elseif(l.eq.lmax(ig)) then
174	fm(ig,l+1)=0.
175	detr(ig,l)=fm(ig,l)+entr(ig,l)
176	else
177	fm(ig,l+1)=0.
178	endif
179	enddo
180
181
182	!-------------------------------------------------------------------------
183	! Verification de la positivite des flux de masse
184	!-------------------------------------------------------------------------
185
186	! do l=1,klev
187	do ig=1,ngrid
188	if (fm(ig,l+1).lt.0.) then
189	! print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
190	ncorecfm1=ncorecfm1+1
191	fm(ig,l+1)=fm(ig,l)
192	detr(ig,l)=entr(ig,l)
193	endif
194	enddo
195	! enddo
196
197	if (prt_level.ge.10) &
198	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
199	& entr(igout,l),detr(igout,l),fm(igout,l+1)
200
201	!-------------------------------------------------------------------------
202	!Test sur fraca croissant
203	!-------------------------------------------------------------------------
204	if (iflag_thermals_optflux==0) then
205	! do l=1,klev
206	do ig=1,ngrid
207	if (l.ge.lalim(ig).and.l.le.lmax(ig) &
208	& .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
209	! zzz est le flux en l+1 a frac constant
210	zzz=fm(ig,l)rhobarz(ig,l+1)zw2(ig,l+1) &
211	& /(rhobarz(ig,l)*zw2(ig,l))
212	if (fm(ig,l+1).gt.zzz) then
213	detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz
214	fm(ig,l+1)=zzz
215	ncorecfm4=ncorecfm4+1
216	endif
217	endif
218	enddo
219	! enddo
220	endif
221
222	if (prt_level.ge.10) &
223	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
224	& entr(igout,l),detr(igout,l),fm(igout,l+1)
225
226
227	!-------------------------------------------------------------------------
228	!test sur flux de masse croissant
229	!-------------------------------------------------------------------------
230	if (iflag_thermals_optflux==0) then
231	! do l=1,klev
232	do ig=1,ngrid
233	if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
234	f_old=fm(ig,l+1)
235	fm(ig,l+1)=fm(ig,l)
236	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
237	ncorecfm5=ncorecfm5+1
238	endif
239	enddo
240	! enddo
241	endif
242
243	if (prt_level.ge.10) &
244	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
245	& entr(igout,l),detr(igout,l),fm(igout,l+1)
246
247	!fin 1.eq.0
248	!-------------------------------------------------------------------------
249	!detr ne peut pas etre superieur a fm
250	!-------------------------------------------------------------------------
251
252	if(1.eq.1) then
253
254	! do l=1,klev
255	do ig=1,ngrid
256	if (entr(ig,l)<0.) then
257	print*,'N1 ig,l,entr',ig,l,entr(ig,l)
258	stop 'entr negatif'
259	endif
260	if (detr(ig,l).gt.fm(ig,l)) then
261	ncorecfm6=ncorecfm6+1
262	detr(ig,l)=fm(ig,l)
263	entr(ig,l)=fm(ig,l+1)
264
265	! Dans le cas ou on est au dessus de la couche d'alimentation et que le
266	! detrainement est plus fort que le flux de masse, on stope le thermique.
267	!test:on commente
268	! if (l.gt.lalim(ig)) then
269	! lmax(ig)=l
270	! fm(ig,l+1)=0.
271	! entr(ig,l)=0.
272	! else
273	! ncorecfm7=ncorecfm7+1
274	! endif
275	endif
276
277	if(l.gt.lmax(ig)) then
278	detr(ig,l)=0.
279	fm(ig,l+1)=0.
280	entr(ig,l)=0.
281	endif
282
283	if (entr(ig,l).lt.0.) then
284	print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
285	print*,'entr(ig,l)',entr(ig,l)
286	print*,'fm(ig,l)',fm(ig,l)
287	stop 'probleme dans thermcell flux'
288	endif
289	enddo
290	! enddo
291	endif
292
293
294	if (prt_level.ge.10) &
295	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
296	& entr(igout,l),detr(igout,l),fm(igout,l+1)
297
298	!-------------------------------------------------------------------------
299	!fm ne peut pas etre negatif
300	!-------------------------------------------------------------------------
301
302	! do l=1,klev
303	do ig=1,ngrid
304	if (fm(ig,l+1).lt.0.) then
305	detr(ig,l)=detr(ig,l)+fm(ig,l+1)
306	fm(ig,l+1)=0.
307	! print*,'fm2<0',l+1,lmax(ig)
308	ncorecfm2=ncorecfm2+1
309	endif
310	if (detr(ig,l).lt.0.) then
311	print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig)
312	print*,'detr(ig,l)',detr(ig,l)
313	print*,'fm(ig,l)',fm(ig,l)
314	stop 'probleme dans thermcell flux'
315	endif
316	enddo
317	! enddo
318
319	if (prt_level.ge.10) &
320	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
321	& entr(igout,l),detr(igout,l),fm(igout,l+1)
322
323	!-----------------------------------------------------------------------
324	!la fraction couverte ne peut pas etre superieure a 1
325	!-----------------------------------------------------------------------
326
327
328	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
329	! FH Partie a revisiter.
330	! Il semble qu'etaient codees ici deux optiques dans le cas
331	! F/ (rho *w) > 1
332	! soit limiter la hauteur du thermique en considerant que c'est
333	! la derniere chouche, soit limiter F a rho w.
334	! Dans le second cas, il faut en fait limiter a un peu moins
335	! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
336	! dans thermcell_main et qu'il semble de toutes facons deraisonable
337	! d'avoir des fractions de 1..
338	! Ci dessous, et dans l'etat actuel, le premier des deux if est
339	! sans doute inutile.
340	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
341
342	! do l=1,klev
343	do ig=1,ngrid
344	if (zw2(ig,l+1).gt.1.e-10) then
345	zfm=rhobarz(ig,l+1)zw2(ig,l+1)alphamax
346	if ( fm(ig,l+1) .gt. zfm) then
347	f_old=fm(ig,l+1)
348	fm(ig,l+1)=zfm
349	! zw2(ig,l+1)=0.
350	! zqla(ig,l+1)=0.
351	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
352	! lmax(ig)=l+1
353	! zmax(ig)=zlev(ig,lmax(ig))
354	! print*,'alpha>1',l+1,lmax(ig)
355	ncorecalpha=ncorecalpha+1
356	endif
357	endif
358	enddo
359	! enddo
360	!
361
362
363	if (prt_level.ge.10) &
364	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
365	& entr(igout,l),detr(igout,l),fm(igout,l+1)
366
367	! Fin de la grande boucle sur les niveaux verticaux
368	enddo
369
370	! if (prt_level.ge.10) &
371	! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
372	! & ptimestep,masse,entr,detr,fm,'8 ')
373
374
375	!-----------------------------------------------------------------------
376	! On fait en sorte que la quantite totale d'air entraine dans le
377	! panache ne soit pas trop grande comparee a la masse de la maille
378	!-----------------------------------------------------------------------
379
380	if (1.eq.1) then
381	do l=1,klev-1
382	do ig=1,ngrid
383	eee0=entr(ig,l)
384	ddd0=detr(ig,l)
385	eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
386	ddd=detr(ig,l)-eee
387	if (eee.gt.0.) then
388	ncorecfm3=ncorecfm3+1
389	entr(ig,l)=entr(ig,l)-eee
390	if ( ddd.gt.0.) then
391	! l'entrainement est trop fort mais l'exces peut etre compense par une
392	! diminution du detrainement)
393	detr(ig,l)=ddd
394	else
395	! l'entrainement est trop fort mais l'exces doit etre compense en partie
396	! par un entrainement plus fort dans la couche superieure
397	if(l.eq.lmax(ig)) then
398	detr(ig,l)=fm(ig,l)+entr(ig,l)
399	else
400	if(l.ge.lmax(ig).and.0.eq.1) then
401	print*,'ig,l',ig,l
402	print*,'eee0',eee0
403	print*,'ddd0',ddd0
404	print*,'eee',eee
405	print*,'ddd',ddd
406	print*,'entr',entr(ig,l)
407	print*,'detr',detr(ig,l)
408	print*,'masse',masse(ig,l)
409	print*,'fomass_max',fomass_max
410	print,'masse(ig,l)fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
411	print*,'ptimestep',ptimestep
412	print*,'lmax(ig)',lmax(ig)
413	print*,'fm(ig,l+1)',fm(ig,l+1)
414	print*,'fm(ig,l)',fm(ig,l)
415	stop 'probleme dans thermcell_flux'
416	endif
417	entr(ig,l+1)=entr(ig,l+1)-ddd
418	detr(ig,l)=0.
419	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
420	detr(ig,l)=0.
421	endif
422	endif
423	endif
424	enddo
425	enddo
426	endif
427	!
428	! ddd=detr(ig)-entre
429	!on s assure que tout s annule bien en zmax
430	do ig=1,ngrid
431	fm(ig,lmax(ig)+1)=0.
432	entr(ig,lmax(ig))=0.
433	detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
434	enddo
435
436	!-----------------------------------------------------------------------
437	! Impression du nombre de bidouilles qui ont ete necessaires
438	!-----------------------------------------------------------------------
439
440	!IM 090508 beg
441	! if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) then
442	!
443	! print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
444	! & ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
445	! & ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
446	! & ncorecfm6,'x fm6', &
447	! & ncorecfm7,'x fm7', &
448	! & ncorecfm8,'x fm8', &
449	! & ncorecalpha,'x alpha'
450	! endif
451	!IM 090508 end
452
453	! if (prt_level.ge.10) &
454	! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
455	! & ptimestep,masse,entr,detr,fm,'fin')
456
457
458	return
459	end

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

Download in other formats:

Original Format