Context Navigation

thermcell_flux2.F90 @ 5306

Last change on this file since 5306 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 16.4 KB

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

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

Download in other formats:

Original Format