Context Navigation

lmdz_thermcell_flux2.f90 @ 5449

Last change on this file since 5449 was 5390, checked in by yann meurdesoif, 3 weeks ago

Remove UTF8 character that inihibit fortran parsing with GPU morphosis
Add missing END SUBROUTINE instead of simple END, that inhibit correct parsing with regulat expression parser (quick and dirty parsing)

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.8 KB

Line
1	MODULE lmdz_thermcell_flux2
2	!
3	! $Id: lmdz_thermcell_flux2.f90 5390 2024-12-05 16:09:25Z fhourdin $
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 SUBROUTINE thermcell_flux2
515	END MODULE lmdz_thermcell_flux2

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

Download in other formats:

Original Format