Context Navigation

lmdz_thermcell_flux2.f90 @ 5512

Last change on this file since 5512 was 5512, checked in by yann meurdesoif, 42 hours ago

Implement GPU automatic port for :

Thermics
acama_gwd_rando
flott_gwd_rando

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

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

Download in other formats:

Original Format