source: LMDZ5/trunk/libf/phylmd/thermcell_flux.F90 @ 2271

Last change on this file since 2271 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.7 KB
Line 
1!
2! $Id: thermcell_flux.F90 1907 2013-11-26 13:10:46Z musat $
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
Note: See TracBrowser for help on using the repository browser.