My Project
 All Classes Files Functions Variables Macros
thermcell_flux2.F90
Go to the documentation of this file.
1 !
2 ! $Id: thermcell_flux2.F90 1403 2010-07-01 09:02:53Z fairhead $
3 !
4  SUBROUTINE thermcell_flux2(ngrid,klev,ptimestep,masse, &
5  & lalim,lmax,alim_star, &
6  & entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
7  & detr,zqla,lev_out,lunout1,igout)
8 !IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
9 
10 
11 !---------------------------------------------------------------------------
12 !thermcell_flux: deduction des flux
13 !---------------------------------------------------------------------------
14 
15  IMPLICIT NONE
16 #include "iniprint.h"
17 #include "thermcell.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,lout
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 
52  logical check_debug,labort_gcm
53 
54  character (len=20) :: modname='thermcell_flux2'
55  character (len=80) :: abort_message
56 
57  fomass_max=0.5
58  alphamax=0.7
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,klev
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_gcm=.true.
111  CALL abort_gcm(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,klev
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,klev
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,klev
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,klev,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,klev
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,klev
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,klev
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,klev
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,klev
269 
270 
271 
272  labort_gcm=.false.
273  do ig=1,ngrid
274  if (entr(ig,l)<0.) then
275  labort_gcm=.true.
276  igout=ig
277  lout=l
278  endif
279  enddo
280 
281  if (labort_gcm) then
282  print*,'N1 ig,l,entr',igout,lout,entr(igout,lout)
283  abort_message = 'entr negatif'
284  CALL abort_gcm(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_gcm=.false.
313  do ig=1,ngrid
314  if (entr(ig,l).lt.0.) then
315  labort_gcm=.true.
316  igout=ig
317  endif
318  enddo
319  if (labort_gcm) 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_gcm(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,klev
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_gcm=.false.
351  do ig=1,ngrid
352  if (detr(ig,l).lt.0.) then
353  labort_gcm=.true.
354  igout=ig
355  endif
356  enddo
357  if (labort_gcm) 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_gcm(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,klev
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 ! zw2(ig,l+1)=0.
398 ! zqla(ig,l+1)=0.
399  detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
400 ! lmax(ig)=l+1
401 ! zmax(ig)=zlev(ig,lmax(ig))
402 ! print*,'alpha>1',l+1,lmax(ig)
403  ncorecalpha=ncorecalpha+1
404  endif
405  endif
406  enddo
407 ! enddo
408 !
409 
410 
411  if (prt_level.ge.10) &
412  & write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
413  & entr(igout,l),detr(igout,l),fm(igout,l+1)
414 
415 ! Fin de la grande boucle sur les niveaux verticaux
416  enddo
417 
418 ! if (prt_level.ge.10) &
419 ! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
420 ! & ptimestep,masse,entr,detr,fm,'8 ')
421 
422 
423 !-----------------------------------------------------------------------
424 ! On fait en sorte que la quantite totale d'air entraine dans le
425 ! panache ne soit pas trop grande comparee a la masse de la maille
426 !-----------------------------------------------------------------------
427 
428  if (1.eq.1) then
429  labort_gcm=.false.
430  do l=1,klev-1
431  do ig=1,ngrid
432  eee0=entr(ig,l)
433  ddd0=detr(ig,l)
434  eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
435  ddd=detr(ig,l)-eee
436  if (eee.gt.0.) then
437  ncorecfm3=ncorecfm3+1
438  entr(ig,l)=entr(ig,l)-eee
439  if ( ddd.gt.0.) then
440 ! l'entrainement est trop fort mais l'exces peut etre compense par une
441 ! diminution du detrainement)
442  detr(ig,l)=ddd
443  else
444 ! l'entrainement est trop fort mais l'exces doit etre compense en partie
445 ! par un entrainement plus fort dans la couche superieure
446  if(l.eq.lmax(ig)) then
447  detr(ig,l)=fm(ig,l)+entr(ig,l)
448  else
449  if(l.ge.lmax(ig).and.0.eq.1) then
450  igout=ig
451  lout=l
452  labort_gcm=.true.
453  endif
454  entr(ig,l+1)=entr(ig,l+1)-ddd
455  detr(ig,l)=0.
456  fm(ig,l+1)=fm(ig,l)+entr(ig,l)
457  detr(ig,l)=0.
458  endif
459  endif
460  endif
461  enddo
462  enddo
463  if (labort_gcm) then
464  ig=igout
465  l=lout
466  print*,'ig,l',ig,l
467  print*,'eee0',eee0
468  print*,'ddd0',ddd0
469  print*,'eee',eee
470  print*,'ddd',ddd
471  print*,'entr',entr(ig,l)
472  print*,'detr',detr(ig,l)
473  print*,'masse',masse(ig,l)
474  print*,'fomass_max',fomass_max
475  print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
476  print*,'ptimestep',ptimestep
477  print*,'lmax(ig)',lmax(ig)
478  print*,'fm(ig,l+1)',fm(ig,l+1)
479  print*,'fm(ig,l)',fm(ig,l)
480  abort_message = 'probleme dans thermcell_flux'
481  CALL abort_gcm(modname,abort_message,1)
482  endif
483  endif
484 !
485 ! ddd=detr(ig)-entre
486 !on s assure que tout s annule bien en zmax
487  do ig=1,ngrid
488  fm(ig,lmax(ig)+1)=0.
489  entr(ig,lmax(ig))=0.
490  detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
491  enddo
492 
493 !-----------------------------------------------------------------------
494 ! Impression du nombre de bidouilles qui ont ete necessaires
495 !-----------------------------------------------------------------------
496 
497 !IM 090508 beg
498 ! if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) then
499 !
500 ! print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
501 ! & ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
502 ! & ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
503 ! & ncorecfm6,'x fm6', &
504 ! & ncorecfm7,'x fm7', &
505 ! & ncorecfm8,'x fm8', &
506 ! & ncorecalpha,'x alpha'
507 ! endif
508 !IM 090508 end
509 
510 ! if (prt_level.ge.10) &
511 ! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
512 ! & ptimestep,masse,entr,detr,fm,'fin')
513 
514 
515  return
516  end