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