source: LMDZ4/trunk/libf/phytherm/thermcell_flux.F90 @ 852

Last change on this file since 852 was 852, checked in by Laurent Fairhead, 17 years ago

Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.9 KB
Line 
1      SUBROUTINE thermcell_flux(ngrid,klev,ptimestep,masse, &
2     &       lalim,lmax,alim_star,  &
3     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
4     &       detr,zqla,zmax,lev_out,lunout,igout)
5
6
7!---------------------------------------------------------------------------
8!thermcell_flux: deduction des flux
9!---------------------------------------------------------------------------
10
11      IMPLICIT NONE
12     
13      INTEGER ig,l
14      INTEGER ngrid,klev
15     
16      REAL alim_star(ngrid,klev),entr_star(ngrid,klev)
17      REAL detr_star(ngrid,klev)
18      REAL zw2(ngrid,klev+1)
19      REAL zlev(ngrid,klev+1)
20      REAL masse(ngrid,klev)
21      REAL ptimestep
22      REAL rhobarz(ngrid,klev)
23      REAL f(ngrid)
24      INTEGER lmax(ngrid)
25      INTEGER lalim(ngrid)
26      REAL zqla(ngrid,klev)
27      REAL zmax(ngrid)
28
29      integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
30      integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
31     
32
33      REAL entr(ngrid,klev),detr(ngrid,klev)
34      REAL fm(ngrid,klev+1)
35      REAL zfm
36
37      integer igout
38      integer lev_out
39      integer lunout
40
41      REAL f_old,ddd0,eee0,ddd,eee,zzz
42
43      REAL fomass_max,alphamax
44      save fomass_max,alphamax
45
46      fomass_max=0.5
47      alphamax=0.7
48
49      ncorecfm1=0
50      ncorecfm2=0
51      ncorecfm3=0
52      ncorecfm4=0
53      ncorecfm5=0
54      ncorecfm6=0
55      ncorecfm7=0
56      ncorecfm8=0
57      ncorecalpha=0
58
59!initialisation
60      fm(:,:)=0.
61     
62      if (lev_out.ge.10) then
63         write(lunout,*) 'Dans thermcell_flux 0'
64         write(lunout,*) 'flux base ',f(igout)
65         write(lunout,*) 'lmax ',lmax(igout)
66         write(lunout,*) 'lalim ',lalim(igout)
67         write(lunout,*) 'ig= ',igout
68         write(lunout,*) ' l E*    A*     D*  '
69         write(lunout,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) &
70     &    ,l=1,lmax(igout))
71      endif
72
73
74!-------------------------------------------------------------------------
75! Verification de la nullite des entrainement et detrainement au dessus
76! de lmax(ig)
77!-------------------------------------------------------------------------
78
79      do l=1,klev
80         do ig=1,ngrid
81            if (l.le.lmax(ig)) then
82               if (entr_star(ig,l).gt.1.) then
83                    print*,'WARNING thermcell_flux 1 ig,l,lmax(ig)',ig,l,lmax(ig)
84                    print*,'entr_star(ig,l)',entr_star(ig,l)
85                    print*,'alim_star(ig,l)',alim_star(ig,l)
86                    print*,'detr_star(ig,l)',detr_star(ig,l)
87!                   stop
88               endif
89            else
90               if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0.) then
91                    print*,'cas 1 : ig,l,lmax(ig)',ig,l,lmax(ig)
92                    print*,'entr_star(ig,l)',entr_star(ig,l)
93                    print*,'alim_star(ig,l)',alim_star(ig,l)
94                    print*,'detr_star(ig,l)',detr_star(ig,l)
95                    stop
96               endif
97            endif
98         enddo
99      enddo
100
101!-------------------------------------------------------------------------
102! Multiplication par le flux de masse issu de la femreture
103!-------------------------------------------------------------------------
104
105      do l=1,klev
106         entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
107         detr(:,l)=f(:)*detr_star(:,l)
108      enddo
109
110      if (lev_out.ge.10) then
111         write(lunout,*) 'Dans thermcell_flux 1'
112         write(lunout,*) 'flux base ',f(igout)
113         write(lunout,*) 'lmax ',lmax(igout)
114         write(lunout,*) 'lalim ',lalim(igout)
115         write(lunout,*) 'ig= ',igout
116         write(lunout,*) ' l   E    D     W2'
117         write(lunout,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) &
118     &    ,zw2(igout,l+1),l=1,lmax(igout))
119      endif
120
121      fm(:,1)=0.
122      do l=1,klev
123         do ig=1,ngrid
124            if (l.lt.lmax(ig)) then
125               fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
126            elseif(l.eq.lmax(ig)) then
127               fm(ig,l+1)=0.
128               detr(ig,l)=fm(ig,l)+entr(ig,l)
129            else
130               fm(ig,l+1)=0.
131            endif
132         enddo
133      enddo
134
135
136
137! Test provisoire : pour comprendre pourquoi on corrige plein de fois
138! le cas fm6, on commence par regarder une premiere fois avant les
139! autres corrections.
140
141      do l=1,klev
142         do ig=1,ngrid
143            if (detr(ig,l).gt.fm(ig,l)) then
144               ncorecfm8=ncorecfm8+1
145!              igout=ig
146            endif
147         enddo
148      enddo
149
150      if (lev_out.ge.10) &
151    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
152    &    ptimestep,masse,entr,detr,fm,'2  ')
153
154!-------------------------------------------------------------------------
155! Verification de la positivite des flux de masse
156!-------------------------------------------------------------------------
157
158      do l=1,klev
159         do ig=1,ngrid
160            if (fm(ig,l+1).lt.0.) then
161!              print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
162                ncorecfm1=ncorecfm1+1
163               fm(ig,l+1)=fm(ig,l)
164               detr(ig,l)=entr(ig,l)
165            endif
166         enddo
167      enddo
168
169      if (lev_out.ge.10) &
170    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
171    &    ptimestep,masse,entr,detr,fm,'3  ')
172
173!-------------------------------------------------------------------------
174!Test sur fraca croissant
175!-------------------------------------------------------------------------
176
177      do l=1,klev
178         do ig=1,ngrid
179          if (l.ge.lalim(ig).and.l.le.lmax(ig) &
180     &    .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
181!  zzz est le flux en l+1 a frac constant
182             zzz=fm(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1)  &
183     &                          /(rhobarz(ig,l)*zw2(ig,l))
184             if (fm(ig,l+1).gt.zzz) then
185                detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz
186                fm(ig,l+1)=zzz
187                ncorecfm4=ncorecfm4+1
188             endif
189          endif
190        enddo
191      enddo
192
193      if (lev_out.ge.10) &
194    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
195    &    ptimestep,masse,entr,detr,fm,'4  ')
196
197
198!-------------------------------------------------------------------------
199!test sur flux de masse croissant
200!-------------------------------------------------------------------------
201
202      do l=1,klev
203         do ig=1,ngrid
204            if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
205              f_old=fm(ig,l+1)
206              fm(ig,l+1)=fm(ig,l)
207              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
208               ncorecfm5=ncorecfm5+1
209            endif
210         enddo
211      enddo
212
213      if (lev_out.ge.10) &
214    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
215    &    ptimestep,masse,entr,detr,fm,'5  ')
216
217!-------------------------------------------------------------------------
218!detr ne peut pas etre superieur a fm
219!-------------------------------------------------------------------------
220
221      if(1.eq.1) then
222
223      do l=1,klev
224         do ig=1,ngrid
225            if (entr(ig,l)<0.) then
226               print*,'N1 ig,l,entr',ig,l,entr(ig,l)
227               stop'entr negatif'
228            endif
229            if (detr(ig,l).gt.fm(ig,l)) then
230               ncorecfm6=ncorecfm6+1
231               detr(ig,l)=fm(ig,l)
232!              entr(ig,l)=fm(ig,l+1)
233
234! Dans le cas ou on est au dessus de la couche d'alimentation et que le
235! detrainement est plus fort que le flux de masse, on stope le thermique.
236               if (l.gt.lalim(ig)) then
237                  lmax(ig)=l
238                  fm(ig,l+1)=0.
239                  entr(ig,l)=0.
240               else
241                  ncorecfm7=ncorecfm7+1
242               endif
243            endif
244
245            if(l.gt.lmax(ig)) then
246               detr(ig,l)=0.
247               fm(ig,l+1)=0.
248               entr(ig,l)=0.
249            endif
250
251            if (entr(ig,l).lt.0.) then
252               print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
253               print*,'entr(ig,l)',entr(ig,l)
254               print*,'fm(ig,l)',fm(ig,l)
255               stop'probleme dans thermcell flux'
256            endif
257         enddo
258      enddo
259      endif
260
261
262      if (lev_out.ge.10) &
263    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
264    &    ptimestep,masse,entr,detr,fm,'6  ')
265
266!-------------------------------------------------------------------------
267!fm ne peut pas etre negatif
268!-------------------------------------------------------------------------
269
270      do l=1,klev
271         do ig=1,ngrid
272            if (fm(ig,l+1).lt.0.) then
273               detr(ig,l)=detr(ig,l)+fm(ig,l+1)
274               fm(ig,l+1)=0.
275!              print*,'fm2<0',l+1,lmax(ig)
276               ncorecfm2=ncorecfm2+1
277            endif
278            if (detr(ig,l).lt.0.) then
279               print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig)
280               print*,'detr(ig,l)',detr(ig,l)
281               print*,'fm(ig,l)',fm(ig,l)
282               stop'probleme dans thermcell flux'
283            endif
284        enddo
285     enddo
286
287      if (lev_out.ge.10) &
288    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
289    &    ptimestep,masse,entr,detr,fm,'7  ')
290
291!-----------------------------------------------------------------------
292!la fraction couverte ne peut pas etre superieure a 1           
293!-----------------------------------------------------------------------
294
295
296!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
297! FH Partie a revisiter.
298! Il semble qu'etaient codees ici deux optiques dans le cas
299! F/ (rho *w) > 1
300! soit limiter la hauteur du thermique en considerant que c'est
301! la derniere chouche, soit limiter F a rho w.
302! Dans le second cas, il faut en fait limiter a un peu moins
303! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
304! dans thermcell_main et qu'il semble de toutes facons deraisonable
305! d'avoir des fractions de 1..
306! Ci dessous, et dans l'etat actuel, le premier des  deux if est
307! sans doute inutile.
308!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
309
310     if(1.eq.0) then
311
312     do l=1,klev
313        do ig=1,ngrid
314           if (zw2(ig,l+1).gt.1.e-10) then
315           zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
316           if ( fm(ig,l+1) .gt. zfm) then
317              f_old=fm(ig,l+1)
318              fm(ig,l+1)=zfm
319!             zw2(ig,l+1)=0.
320!             zqla(ig,l+1)=0.
321              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
322!             lmax(ig)=l+1
323!             zmax(ig)=zlev(ig,lmax(ig))
324!             print*,'alpha>1',l+1,lmax(ig)
325              ncorecalpha=ncorecalpha+1
326           endif
327           endif
328        enddo
329     enddo
330!
331      if (lev_out.ge.10) &
332    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
333    &    ptimestep,masse,entr,detr,fm,'8  ')
334
335     endif
336
337!-----------------------------------------------------------------------
338! On fait en sorte que la quantite totale d'air entraine dans le
339! panache ne soit pas trop grande comparee a la masse de la maille
340!-----------------------------------------------------------------------
341
342      if (1.eq.1) then
343      do l=1,klev-1
344         do ig=1,ngrid
345            eee0=entr(ig,l)
346            ddd0=detr(ig,l)
347            eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
348            ddd=detr(ig,l)-eee
349            if (eee.gt.0.) then
350                ncorecfm3=ncorecfm3+1
351                entr(ig,l)=entr(ig,l)-eee
352                if ( ddd.gt.0.) then
353!   l'entrainement est trop fort mais l'exces peut etre compense par une
354!   diminution du detrainement)
355                   detr(ig,l)=ddd
356                else
357!   l'entrainement est trop fort mais l'exces doit etre compense en partie
358!   par un entrainement plus fort dans la couche superieure
359                   if(l.eq.lmax(ig)) then
360                      detr(ig,l)=fm(ig,l)+entr(ig,l)
361                   else
362                      if(l.ge.lmax(ig).and.0.eq.1) then
363                         print*,'ig,l',ig,l
364                         print*,'eee0',eee0
365                         print*,'ddd0',ddd0
366                         print*,'eee',eee
367                         print*,'ddd',ddd
368                         print*,'entr',entr(ig,l)
369                         print*,'detr',detr(ig,l)
370                         print*,'masse',masse(ig,l)
371                         print*,'fomass_max',fomass_max
372                         print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
373                         print*,'ptimestep',ptimestep
374                         print*,'lmax(ig)',lmax(ig)
375                         print*,'fm(ig,l+1)',fm(ig,l+1)
376                         print*,'fm(ig,l)',fm(ig,l)
377                         stop'probleme dans thermcell_flux'
378                      endif
379                      entr(ig,l+1)=entr(ig,l+1)-ddd
380                      detr(ig,l)=0.
381                      fm(ig,l+1)=fm(ig,l)+entr(ig,l)
382                      detr(ig,l)=0.
383                   endif
384                endif
385            endif
386         enddo
387      enddo
388      endif
389!                 
390!              ddd=detr(ig)-entre
391!on s assure que tout s annule bien en zmax
392      do ig=1,ngrid
393         fm(ig,lmax(ig)+1)=0.
394         entr(ig,lmax(ig))=0.
395         detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
396      enddo
397
398!-----------------------------------------------------------------------
399! Impression du nombre de bidouilles qui ont ete necessaires
400!-----------------------------------------------------------------------
401
402      if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) then
403
404          print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
405    &     ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
406    &     ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
407    &     ncorecfm6,'x fm6', &
408    &     ncorecfm7,'x fm7', &
409    &     ncorecfm8,'x fm8', &
410    &     ncorecalpha,'x alpha'
411      endif
412
413      if (lev_out.ge.10) &
414    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
415    &    ptimestep,masse,entr,detr,fm,'fin')
416
417
418      return
419      end
420
421      subroutine printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
422    &    ptimestep,masse,entr,detr,fm,descr)
423
424     implicit none
425
426      integer ngrid,klev,lunout,igout,l,lm
427
428      integer lmax(klev),lalim(klev)
429      real ptimestep,masse(ngrid,klev),entr(ngrid,klev),detr(ngrid,klev)
430      real fm(ngrid,klev+1),f(ngrid)
431
432      character*3 descr
433
434      lm=lmax(igout)+5
435      if(lm.gt.klev) lm=klev
436
437      print*,'Impression jusque lm=',lm
438
439         write(lunout,*) 'Dans thermcell_flux '//descr
440         write(lunout,*) 'flux base ',f(igout)
441         write(lunout,*) 'lmax ',lmax(igout)
442         write(lunout,*) 'lalim ',lalim(igout)
443         write(lunout,*) 'ig= ',igout
444         write(lunout,'(a3,4a14)') 'l','M','E','D','F'
445         write(lunout,'(i4,4e14.4)') (l,masse(igout,l)/ptimestep, &
446     &     entr(igout,l),detr(igout,l) &
447     &    ,fm(igout,l+1),l=1,lm)
448
449
450      do l=lmax(igout)+1,klev
451          if (abs(entr(igout,l))+abs(detr(igout,l))+abs(fm(igout,l)).gt.0.) then
452          print*,'cas 1 : igout,l,lmax(igout)',igout,l,lmax(igout)
453          print*,'entr(igout,l)',entr(igout,l)
454          print*,'detr(igout,l)',detr(igout,l)
455          print*,'fm(igout,l)',fm(igout,l)
456          stop
457          endif
458      enddo
459
460      return
461      end
462
Note: See TracBrowser for help on using the repository browser.