Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

lmdz_thermcell_flux2.F90

Last change on this file was 5158, checked in by abarral, 3 months ago
Add missing klon on strataer_emiss_mod.F90 Correct various missing explicit declarations Replace tabs by spaces (tabs are not part of the fortran charset) Continue cleaning modules Removed unused arguments and variables
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: 17.2 KB

Line
1	MODULE lmdz_thermcell_flux2
2
3	! $Id: lmdz_thermcell_flux2.F90 5158 2024-08-02 12:12:03Z fairhead $
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	USE lmdz_abort_physic, ONLY: abort_physic
20	IMPLICIT NONE
21
22	! arguments
23	INTEGER, INTENT(IN) :: ngrid, nlay
24	REAL, INTENT(IN) :: ptimestep
25	REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: masse
26	INTEGER, INTENT(IN), DIMENSION(ngrid) :: lalim, lmax
27	REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: alim_star, entr_star, detr_star
28	REAL, INTENT(IN), DIMENSION(ngrid) :: f
29	REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: rhobarz
30	REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: zw2, zlev
31	! FH : laisser ca le temps de verifier qu'on a bien fait de commenter les
32	! lignes faisant apparaitre zqla, zmax ...
33	! REAL, INTENT(IN), DIMENSION(ngrid) :: zmax(ngrid)
34	! enlever aussi zqla
35	REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zqla ! not used
36	INTEGER, INTENT(IN) :: lev_out, lunout1
37
38	REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: entr, detr
39	REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: fm
40
41	! local
42	INTEGER ig, l
43	INTEGER igout, lout
44	REAL zfm
45	INTEGER ncorecfm1, ncorecfm2, ncorecfm3, ncorecalpha
46	INTEGER ncorecfm4, ncorecfm5, ncorecfm6, ncorecfm7, ncorecfm8
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	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>=10) THEN
73	WRITE(lunout1, *) 'Dans thermcell_flux 0'
74	WRITE(lunout1, *) 'flux base ', f(igout)
75	WRITE(lunout1, *) 'lmax ', lmax(igout)
76	WRITE(lunout1, *) 'lalim ', lalim(igout)
77	WRITE(lunout1, *) 'ig= ', igout
78	WRITE(lunout1, ) ' l E A* D* '
79	WRITE(lunout1, '(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	! Active uniquement si check_debug=.TRUE. ou prt_level>=10
88	!-------------------------------------------------------------------------
89
90	check_debug = .false..or.prt_level>=10
91
92	IF (check_debug) THEN
93	DO l = 1, nlay
94	DO ig = 1, ngrid
95	IF (l<=lmax(ig)) THEN
96	IF (entr_star(ig, l)>1.) THEN
97	PRINT*, 'WARNING thermcell_flux 1 ig,l,lmax(ig)', ig, l, lmax(ig)
98	PRINT*, 'entr_star(ig,l)', entr_star(ig, l)
99	PRINT*, 'alim_star(ig,l)', alim_star(ig, l)
100	PRINT*, 'detr_star(ig,l)', detr_star(ig, l)
101	endif
102	else
103	IF (abs(entr_star(ig, l)) + abs(alim_star(ig, l)) + abs(detr_star(ig, l))>0.) THEN
104	PRINT*, 'cas 1 : ig,l,lmax(ig)', ig, l, lmax(ig)
105	PRINT*, 'entr_star(ig,l)', entr_star(ig, l)
106	PRINT*, 'alim_star(ig,l)', alim_star(ig, l)
107	PRINT*, 'detr_star(ig,l)', detr_star(ig, l)
108	abort_message = ''
109	labort_physic = .TRUE.
110	CALL abort_physic (modname, abort_message, 1)
111	endif
112	endif
113	enddo
114	enddo
115	endif
116
117	!-------------------------------------------------------------------------
118	! Multiplication par le flux de masse issu de la femreture
119	!-------------------------------------------------------------------------
120
121	DO l = 1, nlay
122	entr(:, l) = f(:) * (entr_star(:, l) + alim_star(:, l))
123	detr(:, l) = f(:) * detr_star(:, l)
124	enddo
125
126	IF (prt_level>=10) THEN
127	WRITE(lunout1, *) 'Dans thermcell_flux 1'
128	WRITE(lunout1, *) 'flux base ', f(igout)
129	WRITE(lunout1, *) 'lmax ', lmax(igout)
130	WRITE(lunout1, *) 'lalim ', lalim(igout)
131	WRITE(lunout1, *) 'ig= ', igout
132	WRITE(lunout1, *) ' l E D W2'
133	WRITE(lunout1, '(i4,3e15.5)') (l, entr(igout, l), detr(igout, l) &
134	, zw2(igout, l + 1), l = 1, lmax(igout))
135	endif
136
137	fm(:, 1) = 0.
138	DO l = 1, nlay
139	DO ig = 1, ngrid
140	IF (l<lmax(ig)) THEN
141	fm(ig, l + 1) = fm(ig, l) + entr(ig, l) - detr(ig, l)
142	elseif(l==lmax(ig)) THEN
143	fm(ig, l + 1) = 0.
144	detr(ig, l) = fm(ig, l) + entr(ig, l)
145	else
146	fm(ig, l + 1) = 0.
147	endif
148	enddo
149	enddo
150
151
152
153	! Test provisoire : pour comprendre pourquoi on corrige plein de fois
154	! le cas fm6, on commence par regarder une premiere fois avant les
155	! autres corrections.
156
157	DO l = 1, nlay
158	DO ig = 1, ngrid
159	IF (detr(ig, l)>fm(ig, l)) THEN
160	ncorecfm8 = ncorecfm8 + 1
161	! igout=ig
162	endif
163	enddo
164	enddo
165
166	! if (prt_level.ge.10) &
167	! & CALL printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
168	! & ptimestep,masse,entr,detr,fm,'2 ')
169
170
171
172	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173	! FH Version en cours de test;
174	! par rapport a thermcell_flux, on fait une grande boucle sur "l"
175	! et on modifie le flux avec tous les contrôles appliques d'affilee
176	! pour la meme couche
177	! Momentanement, on duplique le calcule du flux pour pouvoir comparer
178	! les flux avant et apres modif
179	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
180
181	DO l = 1, nlay
182
183	DO ig = 1, ngrid
184	IF (l<lmax(ig)) THEN
185	fm(ig, l + 1) = fm(ig, l) + entr(ig, l) - detr(ig, l)
186	elseif(l==lmax(ig)) THEN
187	fm(ig, l + 1) = 0.
188	detr(ig, l) = fm(ig, l) + entr(ig, l)
189	else
190	fm(ig, l + 1) = 0.
191	endif
192	enddo
193
194
195	!-------------------------------------------------------------------------
196	! Verification de la positivite des flux de masse
197	!-------------------------------------------------------------------------
198
199	! do l=1,nlay
200	DO ig = 1, ngrid
201	IF (fm(ig, l + 1)<0.) THEN
202	! PRINT*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
203	ncorecfm1 = ncorecfm1 + 1
204	fm(ig, l + 1) = fm(ig, l)
205	detr(ig, l) = entr(ig, l)
206	endif
207	enddo
208	! enddo
209
210	IF (prt_level>=10) &
211	WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
212	entr(igout, l), detr(igout, l), fm(igout, l + 1)
213
214	!-------------------------------------------------------------------------
215	!Test sur fraca croissant
216	!-------------------------------------------------------------------------
217	IF (iflag_thermals_optflux==0) THEN
218	! do l=1,nlay
219	DO ig = 1, ngrid
220	IF (l>=lalim(ig).AND.l<=lmax(ig) &
221	.AND.(zw2(ig, l + 1)>1.e-10).AND.(zw2(ig, l)>1.e-10)) THEN
222	! zzz est le flux en l+1 a frac constant
223	zzz = fm(ig, l) * rhobarz(ig, l + 1) * zw2(ig, l + 1) &
224	/ (rhobarz(ig, l) * zw2(ig, l))
225	IF (fm(ig, l + 1)>zzz) THEN
226	detr(ig, l) = detr(ig, l) + fm(ig, l + 1) - zzz
227	fm(ig, l + 1) = zzz
228	ncorecfm4 = ncorecfm4 + 1
229	endif
230	endif
231	enddo
232	! enddo
233	endif
234
235	IF (prt_level>=10) &
236	WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
237	entr(igout, l), detr(igout, l), fm(igout, l + 1)
238
239
240	!-------------------------------------------------------------------------
241	!test sur flux de masse croissant
242	!-------------------------------------------------------------------------
243	IF (iflag_thermals_optflux==0) THEN
244	! do l=1,nlay
245	DO ig = 1, ngrid
246	IF ((fm(ig, l + 1)>fm(ig, l)).AND.(l>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	endif
255
256	IF (prt_level>=10) &
257	WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
258	entr(igout, l), detr(igout, l), fm(igout, l + 1)
259
260	!fin 1.EQ.0
261	!-------------------------------------------------------------------------
262	!detr ne peut pas etre superieur a fm
263	!-------------------------------------------------------------------------
264
265	IF(1==1) THEN
266	! do l=1,nlay
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)>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>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)<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	IF (prt_level>=10) &
329	WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
330	entr(igout, l), detr(igout, l), fm(igout, l + 1)
331
332	!-------------------------------------------------------------------------
333	!fm ne peut pas etre negatif
334	!-------------------------------------------------------------------------
335
336	! do l=1,nlay
337	DO ig = 1, ngrid
338	IF (fm(ig, l + 1)<0.) THEN
339	detr(ig, l) = detr(ig, l) + fm(ig, l + 1)
340	fm(ig, l + 1) = 0.
341	ncorecfm2 = ncorecfm2 + 1
342	endif
343	enddo
344
345	labort_physic = .FALSE.
346	DO ig = 1, ngrid
347	IF (detr(ig, l)<0.) THEN
348	labort_physic = .TRUE.
349	igout = ig
350	endif
351	enddo
352	IF (labort_physic) THEN
353	ig = igout
354	PRINT*, 'cas 2 : ig,l,lmax(ig)', ig, l, lmax(ig)
355	PRINT*, 'detr(ig,l)', detr(ig, l)
356	PRINT*, 'fm(ig,l)', fm(ig, l)
357	abort_message = 'probleme dans thermcell flux'
358	CALL abort_physic (modname, abort_message, 1)
359	endif
360	! enddo
361
362	IF (prt_level>=10) &
363	WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
364	entr(igout, l), detr(igout, l), fm(igout, l + 1)
365
366	!-----------------------------------------------------------------------
367	!la fraction couverte ne peut pas etre superieure a 1
368	!-----------------------------------------------------------------------
369
370
371	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
372	! FH Partie a revisiter.
373	! Il semble qu'etaient codees ici deux optiques dans le cas
374	! F/ (rho *w) > 1
375	! soit limiter la hauteur du thermique en considerant que c'est
376	! la derniere chouche, soit limiter F a rho w.
377	! Dans le second cas, il faut en fait limiter a un peu moins
378	! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
379	! dans thermcell_main et qu'il semble de toutes facons deraisonable
380	! d'avoir des fractions de 1..
381	! Ci dessous, et dans l'etat actuel, le premier des deux if est
382	! sans doute inutile.
383	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
384
385	! do l=1,nlay
386	DO ig = 1, ngrid
387	IF (zw2(ig, l + 1)>1.e-10) THEN
388	zfm = rhobarz(ig, l + 1) * zw2(ig, l + 1) * alphamax
389	IF (fm(ig, l + 1) > zfm) THEN
390	f_old = fm(ig, l + 1)
391	fm(ig, l + 1) = zfm
392	detr(ig, l) = detr(ig, l) + f_old - fm(ig, l + 1)
393	! lmax(ig)=l+1
394	! zmax(ig)=zlev(ig,lmax(ig))
395	! PRINT*,'alpha>1',l+1,lmax(ig)
396	ncorecalpha = ncorecalpha + 1
397	endif
398	endif
399	enddo
400	! enddo
401
402	IF (prt_level>=10) &
403	WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
404	entr(igout, l), detr(igout, l), fm(igout, l + 1)
405
406	! Fin de la grande boucle sur les niveaux verticaux
407	enddo
408
409	! if (prt_level.ge.10) &
410	! & CALL printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
411	! & ptimestep,masse,entr,detr,fm,'8 ')
412
413
414	!-----------------------------------------------------------------------
415	! On fait en sorte que la quantite totale d'air entraine dans le
416	! panache ne soit pas trop grande comparee a la masse de la maille
417	!-----------------------------------------------------------------------
418
419	IF (1==1) THEN
420	labort_physic = .FALSE.
421	DO l = 1, nlay - 1
422	DO ig = 1, ngrid
423	eee0 = entr(ig, l)
424	ddd0 = detr(ig, l)
425	eee = entr(ig, l) - masse(ig, l) * fomass_max / ptimestep
426	ddd = detr(ig, l) - eee
427	IF (eee>0.) THEN
428	ncorecfm3 = ncorecfm3 + 1
429	entr(ig, l) = entr(ig, l) - eee
430	IF (ddd>0.) THEN
431	! l'entrainement est trop fort mais l'exces peut etre compense par une
432	! diminution du detrainement)
433	detr(ig, l) = ddd
434	else
435	! l'entrainement est trop fort mais l'exces doit etre compense en partie
436	! par un entrainement plus fort dans la couche superieure
437	IF(l==lmax(ig)) THEN
438	detr(ig, l) = fm(ig, l) + entr(ig, l)
439	else
440	IF(l>=lmax(ig).AND.0==1) THEN
441	igout = ig
442	lout = l
443	labort_physic = .TRUE.
444	endif
445	entr(ig, l + 1) = entr(ig, l + 1) - ddd
446	detr(ig, l) = 0.
447	fm(ig, l + 1) = fm(ig, l) + entr(ig, l)
448	detr(ig, l) = 0.
449	endif
450	endif
451	endif
452	enddo
453	enddo
454	IF (labort_physic) THEN
455	ig = igout
456	l = lout
457	PRINT*, 'ig,l', ig, l
458	PRINT*, 'eee0', eee0
459	PRINT*, 'ddd0', ddd0
460	PRINT*, 'eee', eee
461	PRINT*, 'ddd', ddd
462	PRINT*, 'entr', entr(ig, l)
463	PRINT*, 'detr', detr(ig, l)
464	PRINT*, 'masse', masse(ig, l)
465	PRINT*, 'fomass_max', fomass_max
466	PRINT, 'masse(ig,l)fomass_max/ptimestep', masse(ig, l) * fomass_max / ptimestep
467	PRINT*, 'ptimestep', ptimestep
468	PRINT*, 'lmax(ig)', lmax(ig)
469	PRINT*, 'fm(ig,l+1)', fm(ig, l + 1)
470	PRINT*, 'fm(ig,l)', fm(ig, l)
471	abort_message = 'probleme dans thermcell_flux'
472	CALL abort_physic (modname, abort_message, 1)
473	endif
474	endif
475
476	! ddd=detr(ig)-entre
477	!on s assure que tout s annule bien en zmax
478	DO ig = 1, ngrid
479	fm(ig, lmax(ig) + 1) = 0.
480	entr(ig, lmax(ig)) = 0.
481	detr(ig, lmax(ig)) = fm(ig, lmax(ig)) + entr(ig, lmax(ig))
482	enddo
483
484	!-----------------------------------------------------------------------
485	! Impression du nombre de bidouilles qui ont ete necessaires
486	!-----------------------------------------------------------------------
487
488	!IM 090508 beg
489	! if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) THEN
490	! PRINT*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
491	! & ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
492	! & ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
493	! & ncorecfm6,'x fm6', &
494	! & ncorecfm7,'x fm7', &
495	! & ncorecfm8,'x fm8', &
496	! & ncorecalpha,'x alpha'
497	! endif
498	!IM 090508 end
499
500	! if (prt_level.ge.10) &
501	! & CALL printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
502	! & ptimestep,masse,entr,detr,fm,'fin')
503
504	RETURN
505	END
506	END MODULE lmdz_thermcell_flux2

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

Download in other formats:

Original Format