source: LMDZ4/trunk/libf/phylmd/thermcell_flux2.F90 @ 972

Last change on this file since 972 was 972, checked in by lmdzadmin, 16 years ago

Version thermique FH/CRio
Ajout tests cas physiques non pris en comptes et ajout/enleve prints
Nouvelle routine thermcell_flux2.F90
IM

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