source: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/thermcell_flux2.F90 @ 5160

Last change on this file since 5160 was 3102, checked in by oboucher, 7 years ago

Removing x permission from these files

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