Context Navigation

thermcell_flux.F90

Last change on this file was 1907, checked in by lguez, 11 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: 16.7 KB

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

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

Download in other formats:

Original Format