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

Last change on this file since 2647 was 2647, checked in by dbardet, 3 years ago

Thermal plume model: Fix an index error for calculating the maximal authorized interlayer mass flux of the plume.

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