My Project
 All Classes Files Functions Variables Macros
thermcell_flux.F90
Go to the documentation of this file.
1 !
2 ! $Id: thermcell_flux.F90 1403 2010-07-01 09:02:53Z fairhead $
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  IMPLICIT NONE
17 #include "iniprint.h"
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_gcm(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_gcm(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_gcm(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_gcm(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_gcm(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_gcm(modname,abort_message,1)
514  endif
515  enddo
516 
517  return
518  end
519