source: LMDZ5/trunk/libf/phylmd/thermcell_flux2.F90 @ 2102

Last change on this file since 2102 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:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 16.3 KB
Line 
1!
2! $Id: thermcell_flux2.F90 1907 2013-11-26 13:10:46Z musat $
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
Note: See TracBrowser for help on using the repository browser.