Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 2102)
+++ trunk/LMDZ.GENERIC/README	(revision 2103)
@@ -1457,3 +1457,6 @@
 
 ==18/02/2018 == AB
-- use detr as thermcell_dq argument (called in thermcell_main) and clean up thermcell_dq
+- use detr as thermcell_dq argument (called in thermcell_main) and clean up thermcell_dq.
+- fix a bug in thermcell_flux with lmax.
+- fix a bug in thermcell_flux when updraft fraction is too high.
+- add new correction in thermcell_flux when entrainement is too high.
Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90	(revision 2102)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90	(revision 2103)
@@ -56,5 +56,5 @@
 !      ------
       
-      INTEGER ig,l
+      INTEGER ig,l,k
       INTEGER lout
       
@@ -66,5 +66,5 @@
       REAL ddd                            ! ddd0 - eee
       REAL zzz                            ! temporary variable set to mass flux
-      
+      REAL fact
       REAL test
       
@@ -75,6 +75,4 @@
       INTEGER ncorecfm5
       INTEGER ncorecfm6
-      INTEGER ncorecfm7
-      INTEGER ncorecfm8
       INTEGER ncorecalpha
       
@@ -95,6 +93,4 @@
       ncorecfm5 = 0
       ncorecfm6 = 0
-      ncorecfm7 = 0
-      ncorecfm8 = 0
       
       ncorecalpha = 0
@@ -178,6 +174,74 @@
          IF (labort_physic) THEN
             print *, 'ERROR: mass flux has negative value(s)!'
-            print *, 'ig,l,fm',ig,l,entr(igout,lout)
-            abort_message = 'negative incoming fm in thermcell_flux'
+            print *, 'ig,l,fm',igout, lout, fm(igout,lout)
+            print *, 'lmin,lmax', lmin(igout), lmax(igout)
+            print *, '-------------------------------'
+            print *, 'fm+', fm(igout,lout+1)
+            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
+            print *, 'fm ', fm(igout,lout)
+            print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
+            print *, 'fm-', fm(igout,lout-1)
+            abort_message = 'Negative incoming fm in thermcell_flux'
+            CALL abort_physic(modname,abort_message,1)
+         ENDIF
+         
+!------------------------------------------------------------------------------
+! Test entrainment positif
+!------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB : According to the previous derivations, we MUST have positive entrainment
+!      and detrainment everywhere!
+!      Indeed, they are set to max between zero and a computed value.
+!      Then tests are uncompromising.
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         
+         labort_physic=.false.
+         
+         DO ig=1,ngrid
+            IF (entr(ig,l).lt.0.) THEN
+               labort_physic = .true.
+               igout = ig
+               lout = l
+            ENDIF
+         ENDDO
+         
+         IF (labort_physic) THEN
+            print *, 'ERROR: entrainment has negative value(s)!'
+            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
+            print *, 'lmin,lmax', lmin(igout), lmax(igout)
+            print *, '-------------------------------'
+            print *, 'fm+', fm(igout,lout+1)
+            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
+            print *, 'fm ', fm(igout,lout)
+            print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
+            print *, 'fm-', fm(igout,lout-1)
+            abort_message = 'Negative entrainment in thermcell_flux'
+            CALL abort_physic(modname,abort_message,1)
+         ENDIF
+         
+!------------------------------------------------------------------------------
+! Test detrainment positif
+!------------------------------------------------------------------------------
+         
+         DO ig=1,ngrid
+            IF (detr(ig,l).lt.0.) THEN
+               labort_physic = .true.
+               igout = ig
+               lout = l
+            ENDIF
+         ENDDO
+         
+         IF (labort_physic) THEN
+            print *, 'ERROR: detrainment has negative value(s)!'
+            print *, 'ig,l,detr', igout, lout, detr(igout,lout)
+            print *, 'lmin,lmax', lmin(igout), lmax(igout)
+            print *, '-------------------------------'
+            print *, 'fm+', fm(igout,lout+1)
+            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
+            print *, 'fm ', fm(igout,lout)
+            print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
+            print *, 'fm-', fm(igout,lout-1)
+            abort_message = 'Negative detrainment in thermcell_flux'
             CALL abort_physic(modname,abort_message,1)
          ENDIF
@@ -220,47 +284,4 @@
          
 !------------------------------------------------------------------------------
-! Test entrainment et detrainment positif
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : According to the previous derivations, we MUST have positive entrainment
-!      and detrainment everywhere!
-!      Indeed, they are set to max between zero and a computed value.
-!      Then tests are uncompromising.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         
-         labort_physic=.false.
-         
-         DO ig=1,ngrid
-            IF (entr(ig,l).lt.0.) THEN
-               labort_physic = .true.
-               igout = ig
-               lout = l
-            ENDIF
-         ENDDO
-         
-         IF (labort_physic) THEN
-            print *, 'ERROR: entrainment has negative value(s)!'
-            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
-            abort_message = 'negative entrainment in thermcell_flux'
-            CALL abort_physic(modname,abort_message,1)
-         ENDIF
-         
-         DO ig=1,ngrid
-            IF (detr(ig,l).lt.0.) THEN
-               labort_physic = .true.
-               igout = ig
-               lout = l
-            ENDIF
-         ENDDO
-         
-         IF (labort_physic) THEN
-            print *, 'ERROR: detrainment has negative value(s)!'
-            print *, 'ig,l,detr', igout, lout, detr(igout,lout)
-            abort_message = 'negative detrainment in thermcell_flux'
-            CALL abort_physic(modname,abort_message,1)
-         ENDIF
-         
-!------------------------------------------------------------------------------
 ! Test detrainment inferieur au flux de masse
 !------------------------------------------------------------------------------
@@ -281,5 +302,5 @@
                
                IF (prt_level.ge.10) THEN
-                  print *, 'Detrainment is modified in thermcell_flux'
+                  print *, 'WARNING: Detrainment is modified in thermcell_flux!'
                   print *, 'ig,l,lmax', ig, l, lmax(ig)
                ENDIF
@@ -299,6 +320,5 @@
 !
 ! AB : Do we have to stop the plume here?
-!
-! AB : WARNING: if lmax is modified, be sure that fm, zw2, entr and detr are
+! AB : Attention, if lmax is modified, be sure that fm, zw2, entr and detr are
 !      set to zero above this level!
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -308,5 +328,5 @@
 !                  entr(ig,l)=0.
 !               ELSE
-!                  ncorecfm7=ncorecfm7+1
+!                  ncorecfm2=ncorecfm2+1
 !               ENDIF
             ENDIF
@@ -332,6 +352,8 @@
 !            print *, '       entrainment cannot compensate it!'
 !            print *, 'ig,l,lmax(ig)', igout, l, lmax(igout)
-!            print *, 'entr(ig,l)', entr(igout,l)
-!            print *, 'fm(ig,l)', fm(igout,l)
+!            print *, '--------------------'
+!            print *, 'fm+', fm(igout,lout+1)
+!            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
+!            print *, 'fm-', fm(igout,lout)
 !            abort_message = 'problem in thermcell_flux'
 !            CALL abort_physic (modname,abort_message,1)
@@ -348,5 +370,8 @@
 !      If it's not enough, we can increase entr in the layer above and decrease
 !      the outgoing mass flux in the current layer.
-!      If it's still unsufficient, code will abort.
+!      If it's still unsufficient, code will abort (now commented).
+!      Else we reset entr to its intial value and we renormalize entrainment,
+!      detrainment and mass flux profiles such as we do not exceed the maximal
+!      authorized entrained mass.
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          
@@ -363,20 +388,32 @@
                ncorecfm3  = ncorecfm3 + 1
                
-               print *, 'CORRECTION: entr is to large!'
-               print *, 'ig,l,lmax', ig, l, lmax(ig)
-               
-               IF (ddd.gt.0.) THEN
+               IF (prt_level.ge.10) THEN
+                  print *, 'WARNING: Entrainment is modified in thermcell_flux!'
+                  print *, 'ig,l,lmax', ig, l, lmax(ig)
+               ENDIF
+               
+               IF (ddd.gt.0.) THEN                          ! on reduit le detr dans la couche
                   detr(ig,l) = ddd
+               ELSEIF (l.eq.lmax(ig)) THEN                  ! on ajuste le detr dans la couche
+                  detr(ig,l) = fm(ig,l) + entr(ig,l)
+               ELSEIF (entr(ig,l+1).gt.ABS(ddd)) THEN       ! on reduit le detr dans la couche et l'entr au-dessus
+                  detr(ig,l) = 0.
+                  fm(ig,l+1) = fm(ig,l) + entr(ig,l)
+                  entr(ig,l+1) = entr(ig,l+1) + ddd
+!               ELSE
+!                  entr(ig,l) = entr(ig,l) + eee
+!                  igout=ig
+!                  lout=l
+!                  labort_physic=.true.
                ELSE
-                  IF (l.eq.lmax(ig)) THEN
-                     detr(ig,l) = fm(ig,l) + entr(ig,l)
-                  ELSEIF (l.lt.lmax(ig)) THEN
-                     entr(ig,l+1) = entr(ig,l+1) - ddd
-                     detr(ig,l) = 0.
-                     fm(ig,l+1) = fm(ig,l) + entr(ig,l)
-                  ELSE
-                     igout=ig
-                     lout=l
-                     labort_physic=.true.
+                  fact = (eee0 - eee) / eee0
+                  entr(ig,l) = eee0
+                  detr(ig,:) = detr(ig,:) * fact
+                  entr(ig,:) = entr(ig,:) * fact
+                  fm(ig,:)   = fm(ig,:) * fact
+                  
+                  IF (prt_level.ge.1) THEN
+                     print *, 'WARNING: Normalisation is modified in thermcell_flux!'
+                     print *, 'old f, new f :', f(ig), fact * f(ig)
                   ENDIF
                ENDIF
@@ -385,18 +422,25 @@
          
          IF (labort_physic) THEN
-            print *, ' ERROR: Entrainment is greater than its maximal authorized value!'
-            print *, '        Nor detrainment neither entrainment can compensate it!'
-            print *, 'ig,l,lmax', igout, lout, lmax(igout)
-            print *, 'entr  :', eee0
-            print *, 'detr  :', ddd0
+            print *, 'ERROR: Entrainment is greater than its maximal authorized value!'
+            print *, '       Nor detrainment neither entrainment can compensate it!'
+            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
+            print *, 'lmin,lmax', lmin(igout), lmax(igout)
+            print *, '--------------------'
             print *, 'e_max :', masse(igout,lout)*fomass_max/ptimestep
-            print *, ' --- fomass_max :', fomass_max
-            print *, ' --- masse      :', masse(igout,lout)
-            print *, ' --- ptimestep  :', ptimestep
-            print *, 'eee :', eee
-            print *, 'ddd :', ddd
-            print *, 'fm-', fm(ig,l)
-            print *, 'entr,detr', entr(ig,l), detr(ig,l)
-            print *, 'fm+', fm(ig,l+1)
+            print *, '   fomass_max :', fomass_max
+            print *, '   masse      :', masse(igout,lout)
+            print *, '   ptimestep  :', ptimestep
+            print *, 'norm  :', f(igout)
+            print *, 'alim* :', alim_star(igout,lout)
+            print *, 'entr* :', entr_star(igout,lout)
+            print *, '--------------------'
+            print *, 'fm l+1   ', fm(igout,lout+2)
+            print *, 'entr <-- ', entr(igout,lout+1)
+            print *, 'detr  -->', detr(igout,lout+1)
+            print *, 'fm l     ', fm(igout,lout+1)
+            print *, 'entr <-- ', entr(igout,lout)
+            print *, 'detr  -->', detr(igout,lout)
+            print *, 'fm l-1   ', fm(igout,lout)
+            call flush
             abort_message = 'problem in thermcell_flux'
             CALL abort_physic (modname,abort_message,1)
@@ -429,7 +473,10 @@
 ! AB : The previous change doesn't observe the equation df/dz = e - d.
 !      To avoid this issue, what is better to do? I choose to increase the
-!      entrainment in the layer above.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-               entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1)
+!      entrainment in the layer above when l<lmax. If l=lmax then fm(l+1)=0 and
+!      there are never any problems.
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+               IF (l.lt.lmax(ig)) THEN
+                  entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1)
+               ENDIF
                
                IF (prt_level.ge.10) THEN
@@ -437,5 +484,5 @@
                   print *, 'ig,l,lmax', ig, l, lmax(ig)
                ENDIF
-                  
+               
                IF (prt_level.ge.100) THEN
                   print *, 'fm-', fm(ig,l)
@@ -444,42 +491,8 @@
                ENDIF
                
-               ncorecalpha = ncorecalpha+1
-               
-            ENDIF
-         ENDDO
-         
-!------------------------------------------------------------------------------
-! Test sur fm sortant positif
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : Is it usefull to do that?
-!      We will check fm in the first test with the next value of l
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!         DO ig=1,ngrid
-!            IF (fm(ig,l+1).lt.0.) THEN
-!               detr(ig,l) = detr(ig,l) + fm(ig,l+1)
-!               fm(ig,l+1) = 0.
-!               ncorecfm2  = ncorecfm2 + 1
-!            ENDIF
-!         ENDDO
-!         
-!         labort_physic=.false.
-!         
-!         DO ig=1,ngrid
-!            IF (detr(ig,l).lt.0.) THEN
-!               labort_physic = .true.
-!               igout = ig
-!               lout = l
-!            ENDIF
-!         ENDDO
-!         
-!         IF (labort_physic) THEN
-!            print *, ' ERROR: mass flux has negative value(s) and detrainement'
-!            print *, '        cannot compensate it!'
-!            print*,'ig,l,entr', igout, lout, entr(igout,lout)
-!            abort_message = 'negative outgoing fm in thermcell flux'
-!            CALL abort_physic (modname,abort_message,1)
-!         ENDIF
+               ncorecalpha = ncorecalpha + 1
+               
+            ENDIF
+         ENDDO
          
       ENDDO
@@ -491,12 +504,13 @@
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : test added to avoid problem when lmax = 0, which is not a fortran index.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      IF (lmax(ig).gt.0) THEN
-         DO ig=1,ngrid
+!      Theoretically, we always get lmax >= lmin >= linf > 0
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      DO ig=1,ngrid
+         IF (lmax(ig).gt.0) THEN
             detr(ig,lmax(ig)) = fm(ig,lmax(ig)) + entr(ig,lmax(ig))
             fm(ig,lmax(ig)+1) = 0.
             entr(ig,lmax(ig)) = 0.
-         ENDDO
-      ENDIF
+         ENDIF
+      ENDDO
       
 !------------------------------------------------------------------------------
@@ -505,10 +519,4 @@
       
       IF (prt_level.ge.1) THEN
-         
-! AB : those outputs are commented for their uselessness.
-!         IF (ncorecfm2.ge.1) THEN
-!            print *, 'WARNING: Outgoing mass flux has negative value(s)!'
-!            print *, 'in', ncorecfm2, 'point(s)'
-!         ENDIF
          
          IF (ncorecfm4.ge.1) THEN
@@ -528,16 +536,16 @@
          
 ! AB : those outputs are commented because corresponding test is also commented.
-!         IF (ncorecfm7.ge.1) THEN
+!         IF (ncorecfm2.ge.1) THEN
 !            print *, 'WARNING: Detrainment is greater than mass flux!'
-!            print *, 'in', ncorecfm7, 'point(s)'
+!            print *, 'in', ncorecfm2, 'point(s)'
 !         ENDIF
          
          IF (ncorecfm3.ge.1) THEN
-            print *, 'WARNING: Entrained mass is greater than its maximal authorized value!'
+            print *, 'WARNING: Entrained mass is greater than maximal authorized value!'
             print *, 'in', ncorecfm3, 'point(s)'
          ENDIF
          
          IF (ncorecalpha.ge.1) THEN
-            print *, 'WARNING: Updraft fraction is greater than its maximal authorized value!'
+            print *, 'WARNING: Updraft fraction is greater than maximal authorized value!'
             print *, 'in', ncorecalpha, 'point(s)'
          ENDIF
@@ -554,7 +562,9 @@
 !               print *, 'WARNING: df/dz != entr - detr'
 !               print *, 'ig,l', ig, l
-!               print *, 'fm-  :', fm(ig,l)
+!               print *, 'fm+  :', fm(ig,l+1)
 !               print *, 'entr,detr', entr(ig,l), detr(ig,l)
-!               print *, 'fm+  :', fm(ig,l+1)
+!               print *, 'fm   :', fm(ig,l)
+!               print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
+!               print *, 'fm-  :', fm(ig,l-1)
 !               print *, 'err. :', test
 !            ENDIF
