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

Last change on this file since 2537 was 2480, checked in by aboissinot, 4 years ago

Commit the last changes in the thermal plume model (more english comments).

File size: 15.5 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+1) * zw2(ig,l+1) * alpha_max
281!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
282! AB: The plume mass flux can be reduced.
283!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
284!            IF (fff0 > fmax) THEN
285!               fm(ig,l+1) = fmax
286!               detr(ig,l) = detr(ig,l) + fff0 - fmax
287!               ncorecalpha = ncorecalpha + 1
288!               entr(ig,l+1) = entr(ig,l+1) + fff0 - fmax
289!            ENDIF
290!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
291! AB: The plume can be stopped here.
292!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
293            IF (fff0 > fmax) THEN
294               ncorecalpha = ncorecalpha + 1
295               DO k=l+1,lmax(ig)
296                  entr(ig,k) = 0.
297                  detr(ig,k) = 0.
298                  fm(ig,k) = 0.
299               ENDDO
300               lmax(ig) = l
301               entr(ig,l) = 0.
302               detr(ig,l) = fm(ig,l)
303            ENDIF
304         ENDDO
305         
306!-------------------------------------------------------------------------------
307! Is detrainment lesser than incoming mass flux ?
308!-------------------------------------------------------------------------------
309         
310!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
311! AB : Even if fm has no negative value, it can be lesser than detr.
312!      That is not suitable because when we will mix the plume with the
313!      environment, it will detrain more mass than it is physically able to do.
314!      When it occures, that imply that entr + fm is greater than detr,
315!      otherwise we get a negative outgoing mass flux (cf. above).
316!      That is why we decrease entrainment and detrainment as follows.
317!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
318         
319         DO ig=1,ngrid
320            IF (detr(ig,l) > fm(ig,l)) THEN
321               detr(ig,l) = fm(ig,l)
322               entr(ig,l) = fm(ig,l+1)
323               ncorecdetr  = ncorecdetr + 1
324            ENDIF
325         ENDDO
326         
327      ENDDO
328     
329!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
330! AB: The renormalization can be applied everywhere.
331!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
332!      IF (fact > 0.) THEN
333!         entr(:,:) = entr(:,:) / fact
334!         detr(:,:) = detr(:,:) / fact
335!         fm(:,:) = fm(:,:) / fact
336!      ENDIF
337     
338!-------------------------------------------------------------------------------
339! Is equation df/dz = e - d still verified ?
340!-------------------------------------------------------------------------------
341     
342!      DO l=1,nlay
343!         DO ig=1,ngrid
344!            test = abs(fm(ig,l) + entr(ig,l) - detr(ig,l) - fm(ig,l+1))
345!            IF (test > 1.e-10) THEN
346!               nerrorequa = nerrorequa + 1
347!            ENDIF
348!         ENDDO
349!      ENDDO
350     
351!-------------------------------------------------------------------------------
352! Reset top boundary conditions
353!-------------------------------------------------------------------------------
354     
355      DO ig=1,ngrid
356         IF (lmax(ig) > 0) THEN
357            detr(ig,lmax(ig)) = fm(ig,lmax(ig))
358            fm(ig,lmax(ig)+1) = 0.
359            entr(ig,lmax(ig)) = 0.
360         ENDIF
361      ENDDO
362     
363!===============================================================================
364! Outputs
365!===============================================================================
366     
367      IF (prt_level > 0) THEN
368         
369         IF (ncorecdetr > 0) THEN
370            print *, 'WARNING: Detrainment is greater than mass flux!'
371            print *, 'In', ncorecdetr, 'grid point(s) over', nlay, 'x', ngrid
372         ENDIF
373         
374         IF (ncorecentr > 0) THEN
375            print *, 'WARNING: Entrained mass is greater than maximal authorized value!'
376            print *, 'In', ncorecentr, 'grid point(s) over', nlay, 'x', ngrid
377         ENDIF
378         
379         IF (ncorecfact > 0) THEN
380            print *, 'WARNING: Entrained mass needs renormalization to be ok!'
381            print *, 'In', ncorecfact, 'grid point(s) over', nlay, 'x', ngrid
382!            print *, 'WARNING: Entr fact:', fact
383         ENDIF
384         
385!         IF (nerrorequa > 0) THEN
386!            print *, 'WARNING: !'
387!            print *, 'in', nerrorequa, 'grid point(s) over', nlay, 'x', ngrid
388!         ENDIF
389         
390         IF (ncorecalpha > 0) THEN
391            print *, 'WARNING: Updraft fraction is greater than maximal authorized value!'
392            print *, 'In', ncorecalpha, 'grid point(s) over', nlay, 'x', ngrid
393         ENDIF
394         
395      ENDIF
396     
397     
398RETURN
399END
Note: See TracBrowser for help on using the repository browser.