Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lmdz_thermcell_flux2.F90 @ 5006

Last change on this file since 5006 was 4590, checked in by fhourdin, 18 months ago
Passage des thermiques a la nouvelle norme.
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	MODULE lmdz_thermcell_flux2
2	!
3	! $Id: lmdz_thermcell_flux2.F90 4590 2023-06-29 01:03:15Z abarral $
4	!
5	CONTAINS
6
7	SUBROUTINE thermcell_flux2(ngrid,nlay,ptimestep,masse, &
8	& lalim,lmax,alim_star, &
9	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
10	& detr,zqla,lev_out,lunout1,igout)
11	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
12
13
14	!---------------------------------------------------------------------------
15	!thermcell_flux: deduction des flux
16	!---------------------------------------------------------------------------
17
18	USE lmdz_thermcell_ini, ONLY : prt_level,iflag_thermals_optflux
19	IMPLICIT NONE
20
21	! arguments
22	INTEGER, intent(in) :: ngrid,nlay
23	REAL, intent(in) :: ptimestep
24	REAL, intent(in), dimension(ngrid,nlay) :: masse
25	INTEGER, intent(in), dimension(ngrid) :: lalim,lmax
26	REAL, intent(in), dimension(ngrid,nlay) :: alim_star,entr_star,detr_star
27	REAL, intent(in), dimension(ngrid) :: f
28	REAL, intent(in), dimension(ngrid,nlay) :: rhobarz
29	REAL, intent(in), dimension(ngrid,nlay+1) :: zw2,zlev
30	! FH : laisser ca le temps de verifier qu'on a bien fait de commenter les
31	! lignes faisant apparaitre zqla, zmax ...
32	! REAL, intent(in), dimension(ngrid) :: zmax(ngrid)
33	! enlever aussi zqla
34	REAL, intent(in), dimension(ngrid,nlay) :: zqla ! not used
35	integer, intent(in) :: lev_out, lunout1
36
37	REAL,intent(out), dimension(ngrid,nlay) :: entr,detr
38	REAL,intent(out), dimension(ngrid,nlay+1) :: fm
39
40	! local
41	INTEGER ig,l
42	integer igout,lout
43	REAL zfm
44	integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
45	integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
46
47
48	REAL f_old,ddd0,eee0,ddd,eee,zzz
49
50	REAL,SAVE :: fomass_max=0.5
51	REAL,SAVE :: alphamax=0.7
52	!$OMP THREADPRIVATE(fomass_max,alphamax)
53
54	logical check_debug,labort_physic
55
56	character (len=20) :: modname='thermcell_flux2'
57	character (len=80) :: abort_message
58
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,nlay
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,nlay
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,nlay
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,nlay
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,nlay,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,nlay
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,nlay
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,nlay
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,nlay
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,nlay
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,nlay
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,nlay
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	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
398	! lmax(ig)=l+1
399	! zmax(ig)=zlev(ig,lmax(ig))
400	! print*,'alpha>1',l+1,lmax(ig)
401	ncorecalpha=ncorecalpha+1
402	endif
403	endif
404	enddo
405	! enddo
406	!
407
408
409	if (prt_level.ge.10) &
410	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
411	& entr(igout,l),detr(igout,l),fm(igout,l+1)
412
413	! Fin de la grande boucle sur les niveaux verticaux
414	enddo
415
416	! if (prt_level.ge.10) &
417	! & call printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
418	! & ptimestep,masse,entr,detr,fm,'8 ')
419
420
421	!-----------------------------------------------------------------------
422	! On fait en sorte que la quantite totale d'air entraine dans le
423	! panache ne soit pas trop grande comparee a la masse de la maille
424	!-----------------------------------------------------------------------
425
426	if (1.eq.1) then
427	labort_physic=.false.
428	do l=1,nlay-1
429	do ig=1,ngrid
430	eee0=entr(ig,l)
431	ddd0=detr(ig,l)
432	eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
433	ddd=detr(ig,l)-eee
434	if (eee.gt.0.) then
435	ncorecfm3=ncorecfm3+1
436	entr(ig,l)=entr(ig,l)-eee
437	if ( ddd.gt.0.) then
438	! l'entrainement est trop fort mais l'exces peut etre compense par une
439	! diminution du detrainement)
440	detr(ig,l)=ddd
441	else
442	! l'entrainement est trop fort mais l'exces doit etre compense en partie
443	! par un entrainement plus fort dans la couche superieure
444	if(l.eq.lmax(ig)) then
445	detr(ig,l)=fm(ig,l)+entr(ig,l)
446	else
447	if(l.ge.lmax(ig).and.0.eq.1) then
448	igout=ig
449	lout=l
450	labort_physic=.true.
451	endif
452	entr(ig,l+1)=entr(ig,l+1)-ddd
453	detr(ig,l)=0.
454	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
455	detr(ig,l)=0.
456	endif
457	endif
458	endif
459	enddo
460	enddo
461	if (labort_physic) then
462	ig=igout
463	l=lout
464	print*,'ig,l',ig,l
465	print*,'eee0',eee0
466	print*,'ddd0',ddd0
467	print*,'eee',eee
468	print*,'ddd',ddd
469	print*,'entr',entr(ig,l)
470	print*,'detr',detr(ig,l)
471	print*,'masse',masse(ig,l)
472	print*,'fomass_max',fomass_max
473	print,'masse(ig,l)fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
474	print*,'ptimestep',ptimestep
475	print*,'lmax(ig)',lmax(ig)
476	print*,'fm(ig,l+1)',fm(ig,l+1)
477	print*,'fm(ig,l)',fm(ig,l)
478	abort_message = 'probleme dans thermcell_flux'
479	CALL abort_physic (modname,abort_message,1)
480	endif
481	endif
482	!
483	! ddd=detr(ig)-entre
484	!on s assure que tout s annule bien en zmax
485	do ig=1,ngrid
486	fm(ig,lmax(ig)+1)=0.
487	entr(ig,lmax(ig))=0.
488	detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
489	enddo
490
491	!-----------------------------------------------------------------------
492	! Impression du nombre de bidouilles qui ont ete necessaires
493	!-----------------------------------------------------------------------
494
495	!IM 090508 beg
496	! if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) then
497	!
498	! print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
499	! & ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
500	! & ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
501	! & ncorecfm6,'x fm6', &
502	! & ncorecfm7,'x fm7', &
503	! & ncorecfm8,'x fm8', &
504	! & ncorecalpha,'x alpha'
505	! endif
506	!IM 090508 end
507
508	! if (prt_level.ge.10) &
509	! & call printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
510	! & ptimestep,masse,entr,detr,fm,'fin')
511
512
513	RETURN
514	end
515	END MODULE lmdz_thermcell_flux2

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

Download in other formats:

Original Format