GCC Code Coverage Report


Directory: ./
File: phys/thermcell_flux.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 207 0.0%
Branches: 0 148 0.0%

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