source: LMDZ4/trunk/libf/phylmd/thermcell_flux.F90 @ 987

Last change on this file since 987 was 987, checked in by Laurent Fairhead, 16 years ago

Du nettoyage sur le parallelisme, inclusion de nouvelles interfaces pour OPA9
et ORCHIDEE YM
LF

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