source: trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90 @ 2176

Last change on this file since 2176 was 2146, checked in by aboissinot, 5 years ago

fix a bug in a thermcell_flux "bidouille".
File size: 13.4 KB
RevLine 
[2060]1!
2!
3!
[2101]4SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
[2132]5                          lmin,lmax,entr_star,detr_star,                      &
[2143]6                          f,rhobarz,zlev,zw2,fm,entr,detr)
[2060]7     
8     
[2127]9!===============================================================================
[2143]10!  Purpose: deduction des flux
[2127]11! 
12!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
13! 
14!===============================================================================
[2060]15     
16      USE print_control_mod, ONLY: prt_level
[2132]17      USE thermcell_mod, ONLY: fomass_max, alpha_max
[2060]18     
19      IMPLICIT NONE
20     
21     
[2127]22!===============================================================================
[2060]23! Declaration
[2127]24!===============================================================================
[2060]25     
[2127]26!     Inputs:
27!     -------
[2060]28     
[2127]29      INTEGER ngrid, nlay
[2060]30      INTEGER lmin(ngrid)
31      INTEGER lmax(ngrid)
32     
[2101]33      REAL entr_star(ngrid,nlay)
34      REAL detr_star(ngrid,nlay)
35      REAL zw2(ngrid,nlay+1)
36      REAL zlev(ngrid,nlay+1)
37      REAL masse(ngrid,nlay)
[2060]38      REAL ptimestep
[2101]39      REAL rhobarz(ngrid,nlay)
[2060]40      REAL f(ngrid)
41     
[2127]42!     Outputs:
43!     --------     
[2060]44     
[2101]45      REAL entr(ngrid,nlay)
46      REAL detr(ngrid,nlay)
47      REAL fm(ngrid,nlay+1)
[2060]48     
[2127]49!     Local:
50!     ------
[2060]51     
[2127]52      INTEGER ig, l, k
[2143]53      INTEGER igout, lout                 ! Error grid point and level
[2060]54     
[2132]55      REAL zfm                            ! Mass flux such as updraft fraction is equal to its maximal authorized value alpha_max
[2127]56      REAL f_old                          ! Save initial value of mass flux
57      REAL eee0                           ! Save initial value of entrainment
58      REAL ddd0                           ! Save initial value of detrainment
[2060]59      REAL eee                            ! eee0 - layer mass * maximal authorized mass fraction / dt
60      REAL ddd                            ! ddd0 - eee
[2127]61      REAL zzz                            ! Temporary variable set to mass flux
[2143]62      REAL fact
63      REAL test
[2060]64     
[2143]65      INTEGER ncorecentr
66      INTEGER ncorecdetr
67      INTEGER nerrorequa
68      INTEGER ncorecfact
69      INTEGER ncorecalpha
[2060]70     
71      LOGICAL labort_physic
72     
[2127]73!===============================================================================
[2060]74! Initialization
[2127]75!===============================================================================
[2060]76     
[2143]77      nerrorequa = 0
78      ncorecentr = 0
79      ncorecdetr = 0
80      ncorecfact = 0
[2060]81      ncorecalpha = 0
82     
83      entr(:,:) = 0.
84      detr(:,:) = 0.
85      fm(:,:)   = 0.
86     
[2143]87      labort_physic = .false.
88     
89      fact = 0.
90     
[2127]91!===============================================================================
[2143]92! Calcul de l'entrainement, du detrainement et du flux de masse
[2127]93!===============================================================================
[2060]94     
[2127]95!-------------------------------------------------------------------------------
[2060]96! Multiplication par la norme issue de la relation de fermeture
[2127]97!-------------------------------------------------------------------------------
[2060]98     
[2101]99      DO l=1,nlay
[2127]100         entr(:,l) = f(:) * entr_star(:,l)
[2060]101         detr(:,l) = f(:) * detr_star(:,l)
102      ENDDO
103     
[2127]104!-------------------------------------------------------------------------------
[2143]105! Calcul du flux de masse
[2127]106!-------------------------------------------------------------------------------
[2060]107     
[2101]108      DO l=1,nlay
[2060]109         DO ig=1,ngrid
[2143]110            IF (l < lmax(ig) .and. l >= lmin(ig)) THEN
[2060]111               fm(ig,l+1) = fm(ig,l) + entr(ig,l) - detr(ig,l)
[2143]112            ELSEIF (l == lmax(ig)) THEN
[2060]113               fm(ig,l+1) = 0.
114               detr(ig,l) = fm(ig,l) + entr(ig,l)
115            ELSE
116               fm(ig,l+1) = 0.
117               entr(ig,l) = 0.
118               detr(ig,l) = 0.
119            ENDIF
120         ENDDO
121      ENDDO
122     
[2127]123!===============================================================================
[2143]124! Checking
[2127]125!===============================================================================
[2060]126     
[2101]127      DO l=1,nlay
[2060]128         
[2127]129!-------------------------------------------------------------------------------
[2132]130! Is incoming mass flux positive ?
[2127]131!-------------------------------------------------------------------------------
[2060]132         
133         DO ig=1,ngrid
[2143]134            IF (fm(ig,l) < 0.) THEN
[2060]135               labort_physic = .true.
136               igout = ig
137               lout = l
138            ENDIF
139         ENDDO
140         
[2127]141!-------------------------------------------------------------------------------
[2132]142! Is entrainment positive ?
[2127]143!-------------------------------------------------------------------------------
[2101]144         
[2060]145         DO ig=1,ngrid
[2143]146            IF (entr(ig,l) < 0.) THEN
[2060]147               labort_physic = .true.
148               igout = ig
149               lout = l
150            ENDIF
151         ENDDO
152         
[2127]153!-------------------------------------------------------------------------------
[2132]154! Is detrainment positive ?
[2127]155!-------------------------------------------------------------------------------
[2103]156         
[2060]157         DO ig=1,ngrid
[2143]158            IF (detr(ig,l) < 0.) THEN
[2060]159               labort_physic = .true.
160               igout = ig
161               lout = l
162            ENDIF
163         ENDDO
164         
[2143]165!-------------------------------------------------------------------------------
166! Abort
167!-------------------------------------------------------------------------------
168         
[2060]169         IF (labort_physic) THEN
[2143]170            print *, '---------------------------------------------------------'
171            print *, 'ERROR: mass flux has negative value(s)!'
172            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
[2103]173            print *, 'lmin,lmax', lmin(igout), lmax(igout)
[2143]174            print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
[2127]175            DO k=nlay,1,-1
176               print *, 'fm ', fm(igout,k+1)
177               print *, 'entr,detr', entr(igout,k), detr(igout,k)
178            ENDDO
[2132]179            print *, 'fm ', fm(igout,1)
[2143]180            print *, '---------------------------------------------------------'
181            CALL abort
[2060]182         ENDIF
183         
[2127]184!-------------------------------------------------------------------------------
[2143]185! Is detrainment lessser than incoming mass flux ?
[2127]186!-------------------------------------------------------------------------------
[2103]187         
[2127]188!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]189! AB : Even if fm has no negative value, it can be lesser than detr.
190!      That's not suitable because when we will mix the plume with the
191!      environment, it will detrain more mass than it is physically able to do.
192!      When it occures, that imply that entr + fm is greater than detr,
193!      otherwise we get a negative outgoing mass flux (cf. above).
194!      That is why we correct the detrainment as follows.
[2127]195!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]196         
197         DO ig=1,ngrid
198            IF (detr(ig,l).gt.fm(ig,l)) THEN
199               detr(ig,l) = fm(ig,l)
200               entr(ig,l) = fm(ig,l+1)
[2132]201               ncorecdetr  = ncorecdetr + 1
[2060]202            ENDIF
203         ENDDO
204         
[2127]205!-------------------------------------------------------------------------------
[2143]206! Is entrained mass lesser than fomass_max ?
[2127]207!-------------------------------------------------------------------------------
[2060]208         
[2127]209!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]210! AB : Entrainment is bigger than the maximal authorized value.
211!      If we consider that the excess entrainement is in fact plume air which
212!      is not detrained then we compensate it by decreasing detr.
213!      If it's not enough, we can increase entr in the layer above and decrease
214!      the outgoing mass flux in the current layer.
[2103]215!      If it's still unsufficient, code will abort (now commented).
216!      Else we reset entr to its intial value and we renormalize entrainment,
217!      detrainment and mass flux profiles such as we do not exceed the maximal
218!      authorized entrained mass.
[2127]219!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]220         
221         DO ig=1,ngrid
222            eee0 = entr(ig,l)
223            ddd0 = detr(ig,l)
224            eee  = entr(ig,l) - masse(ig,l) * fomass_max / ptimestep
225            ddd  = detr(ig,l) - eee
[2143]226            IF (eee > 0.) THEN
[2060]227               entr(ig,l) = entr(ig,l) - eee
[2132]228               ncorecentr  = ncorecentr + 1
[2143]229               IF (ddd > 0.) THEN
[2060]230                  detr(ig,l) = ddd
[2143]231               ELSEIF (l == lmax(ig)) THEN
[2103]232                  detr(ig,l) = fm(ig,l) + entr(ig,l)
[2143]233               ELSEIF (entr(ig,l+1) > ABS(ddd)) THEN
[2103]234                  detr(ig,l) = 0.
235                  fm(ig,l+1) = fm(ig,l) + entr(ig,l)
236                  entr(ig,l+1) = entr(ig,l+1) + ddd
237!               ELSE
238!                  entr(ig,l) = entr(ig,l) + eee
[2143]239!                  igout = ig
240!                  lout = l
241!                  labort_physic = .true.
[2060]242               ELSE
[2143]243                  fact = max(fact, eee0 / (eee0 - eee))
[2103]244                  entr(ig,l) = eee0
[2143]245                  ncorecfact = ncorecfact + 1
[2060]246               ENDIF
247            ENDIF
248         ENDDO
249         
[2143]250         IF (labort_physic) THEN
251            print *, '---------------------------------------------------------'
252            print *, 'ERROR: Entrainment is greater than maximal authorized value!'
253            print *, '       Nor detrainment neither entrainment can compensate it!'
254            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
255            print *, 'lmin,lmax', lmin(igout), lmax(igout)
256            print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
257            print *, 'e_max :', masse(igout,lout)*fomass_max/ptimestep
258            print *, '   fomass_max :', fomass_max
259            print *, '   masse      :', masse(igout,lout)
260            print *, '   ptimestep  :', ptimestep
261            print *, 'norm  :', f(igout)
262            print *, 'entr* :', entr_star(igout,lout)
263            print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
264            DO k=nlay,1,-1
265               print *, 'fm ', fm(igout,k+1)
266               print *, 'entr,detr', entr(igout,k), detr(igout,k)
267            ENDDO
268            print *, 'fm ', fm(igout,1)
269            print *, '---------------------------------------------------------'
270            CALL abort
271         ENDIF
[2060]272         
[2127]273!-------------------------------------------------------------------------------
[2143]274! Is Updraft fraction lesser than alpha_max ?
[2127]275!-------------------------------------------------------------------------------
[2060]276         
277         DO ig=1,ngrid
[2132]278            zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha_max
[2127]279           
[2143]280            IF (fm(ig,l+1) > zfm) THEN
[2060]281               f_old = fm(ig,l+1)
282               fm(ig,l+1) = zfm
283               detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1)
[2143]284               ncorecalpha = ncorecalpha + 1
[2127]285!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]286! AB : The previous change doesn't observe the equation df/dz = e - d.
[2127]287!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2143]288               entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1)
[2060]289            ENDIF
290         ENDDO
291         
292      ENDDO
293     
[2143]294      IF (fact > 0.) THEN
295         entr(:,:) = entr(:,:) / fact
[2146]296         detr(:,:) = detr(:,:) / fact
297         fm(:,:) = fm(:,:) / fact
[2143]298      ENDIF
299     
[2127]300!-------------------------------------------------------------------------------
[2143]301! Is equation df/dz = e - d still verified ?
[2127]302!-------------------------------------------------------------------------------
[2060]303     
[2143]304!      DO l=1,nlay
305!         DO ig=1,ngrid
306!            test = abs(fm(ig,l) + entr(ig,l) - detr(ig,l) - fm(ig,l+1))
307!            IF (test > 1.e-10) THEN
308!               nerrorequa = nerrorequa + 1
309!            ENDIF
310!         ENDDO
311!      ENDDO
312     
313!-------------------------------------------------------------------------------
314! Reset top boundary conditions
315!-------------------------------------------------------------------------------
316     
[2103]317      DO ig=1,ngrid
[2143]318         IF (lmax(ig).GT.0) THEN
[2060]319            detr(ig,lmax(ig)) = fm(ig,lmax(ig)) + entr(ig,lmax(ig))
320            fm(ig,lmax(ig)+1) = 0.
321            entr(ig,lmax(ig)) = 0.
[2103]322         ENDIF
323      ENDDO
[2060]324     
[2143]325!===============================================================================
326! Outputs
327!===============================================================================
[2060]328     
[2143]329      IF (prt_level > 0) THEN
330         
331         IF (ncorecdetr.ge.1) THEN
[2060]332            print *, 'WARNING: Detrainment is greater than mass flux!'
[2143]333            print *, 'In', ncorecdetr, 'grid point(s) over', nlay, 'x', ngrid
[2060]334         ENDIF
335         
[2143]336         IF (ncorecentr.ge.1) THEN
[2103]337            print *, 'WARNING: Entrained mass is greater than maximal authorized value!'
[2143]338            print *, 'in', ncorecentr, 'grid point(s) over', nlay, 'x', ngrid
[2060]339         ENDIF
340         
[2143]341         IF (ncorecfact.ge.1) THEN
342            print *, 'WARNING: Entrained mass needs renormalization to be ok!'
343            print *, 'in', ncorecfact, 'grid point(s) over', nlay, 'x', ngrid
344         ENDIF
345         
346!         IF (nerrorequa.ge.1) THEN
347!            print *, 'WARNING: !'
348!            print *, 'in', nerrorequa, 'grid point(s) over', nlay, 'x', ngrid
349!         ENDIF
350         
351         IF (ncorecalpha.ge.1) THEN
[2103]352            print *, 'WARNING: Updraft fraction is greater than maximal authorized value!'
[2143]353            print *, 'in', ncorecalpha, 'grid point(s) over', nlay, 'x', ngrid
[2060]354         ENDIF
355         
356      ENDIF
357     
358     
359RETURN
360END
Note: See TracBrowser for help on using the repository browser.