Context Navigation

thermcell_flux.F90 @ 4040

Last change on this file since 4040 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:keywords set to `Author Date Id Revision`
File size: 16.2 KB

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

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

Download in other formats:

Original Format