Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90	(revision 2103)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90	(revision 2127)
@@ -3,12 +3,14 @@
 !
 SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
-                          lalim,lmin,lmax,alim_star,entr_star,detr_star,      &
-                          f,rhobarz,zlev,zw2,fm,entr,                         &
-                          detr,zqla,lev_out,lunout1,igout)
-      
-      
-!==============================================================================
-! thermcell_flux: deduction des flux
-!==============================================================================
+                          lalim,lmin,lmax,entr_star,detr_star,                &
+                          f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
+      
+      
+!===============================================================================
+!  Purpose: deduction des flux
+!  
+!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
+!  
+!===============================================================================
       
       USE print_control_mod, ONLY: prt_level
@@ -18,21 +20,18 @@
       
       
-!==============================================================================
+!===============================================================================
 ! Declaration
-!==============================================================================
-      
-!      inputs:
-!      ------- 
-      
-      INTEGER ngrid
-      INTEGER nlay
-      INTEGER igout
-      INTEGER lev_out
-      INTEGER lunout1
+!===============================================================================
+      
+!     Inputs:
+!     ------- 
+      
+      INTEGER ngrid, nlay
       INTEGER lmin(ngrid)
       INTEGER lmax(ngrid)
       INTEGER lalim(ngrid)
       
-      REAL alim_star(ngrid,nlay)
+! AB : I remove alimentation, which is in fact included in entr_star. EN COURS
+!      REAL alim_star(ngrid,nlay)
       REAL entr_star(ngrid,nlay)
       REAL detr_star(ngrid,nlay)
@@ -46,6 +45,6 @@
       REAL zmax(ngrid)
       
-!      outputs:
-!      --------      
+!     Outputs:
+!     --------      
       
       REAL entr(ngrid,nlay)
@@ -53,17 +52,17 @@
       REAL fm(ngrid,nlay+1)
       
-!      local:
-!      ------
-      
-      INTEGER ig,l,k
-      INTEGER lout
-      
-      REAL zfm                            ! mass flux such as updraft fraction is equal to its maximal authorized value alphamax
-      REAL f_old                          ! save initial value of mass flux
-      REAL eee0                           ! save initial value of entrainment
-      REAL ddd0                           ! save initial value of detrainment
+!     Local:
+!     ------
+      
+      INTEGER ig, l, k
+      INTEGER igout, lout                 ! Error grid point number and level
+      
+      REAL zfm                            ! Mass flux such as updraft fraction is equal to its maximal authorized value alphamax
+      REAL f_old                          ! Save initial value of mass flux
+      REAL eee0                           ! Save initial value of entrainment
+      REAL ddd0                           ! Save initial value of detrainment
       REAL eee                            ! eee0 - layer mass * maximal authorized mass fraction / dt
       REAL ddd                            ! ddd0 - eee
-      REAL zzz                            ! temporary variable set to mass flux
+      REAL zzz                            ! Temporary variable set to mass flux
       REAL fact
       REAL test
@@ -83,7 +82,7 @@
       CHARACTER (len=80) :: abort_message
       
-!==============================================================================
+!===============================================================================
 ! Initialization
-!==============================================================================
+!===============================================================================
       
       ncorecfm1 = 0
@@ -100,32 +99,34 @@
       fm(:,:)   = 0.
       
-      IF (prt_level.ge.10) THEN
-         write(lunout1,*) 'Dans thermcell_flux 0'
-         write(lunout1,*) 'flux base ', f(igout)
-         write(lunout1,*) 'lmax ', lmax(igout)
-         write(lunout1,*) 'lalim ', lalim(igout)
-         write(lunout1,*) 'ig= ', igout
-         write(lunout1,*) ' l E*    A*     D*  '
-         write(lunout1,'(i4,3e15.5)') (l,entr_star(igout,l),                  &
-         &                             alim_star(igout,l),detr_star(igout,l), &
-         &                             l=1,lmax(igout))
-      ENDIF
-      
-!==============================================================================
+!===============================================================================
 ! Calcul de l'entrainement, du detrainement et du flux de masse
-!==============================================================================
-      
-!------------------------------------------------------------------------------
+!===============================================================================
+      
+!-------------------------------------------------------------------------------
 ! Multiplication par la norme issue de la relation de fermeture
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       DO l=1,nlay
-         entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l))
+! AB : I remove alimentation, which is in fact included in entr_star. EN COURS
+!         entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l))
+         entr(:,l) = f(:) * entr_star(:,l)
          detr(:,l) = f(:) * detr_star(:,l)
       ENDDO
       
-!------------------------------------------------------------------------------
+! AB : temporary test
+      DO l=1,nlay
+         DO ig=1,ngrid
+            IF (detr(ig,l)<0.) THEN
+               print *, 'ERROR: detr is negative in thermcell_flux!'
+               print *, 'ig,l', ig, l
+               print *, 'detr_star', detr_star(ig,l)
+               print *, 'detr', detr(ig,l)
+            ENDIF
+         ENDDO
+      ENDDO
+      
+!-------------------------------------------------------------------------------
 ! Calcul du flux de masse
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       DO l=1,nlay
@@ -144,15 +145,15 @@
       ENDDO
       
-!==============================================================================
+!===============================================================================
 ! Tests de validite
-!==============================================================================
+!===============================================================================
       
       DO l=1,nlay
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Verification de la positivite du flux de masse entrant
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : I remove the correction and replace it by an uncompromising test.
 !      According to the previous derivations, we MUST have positive mass flux
@@ -160,5 +161,5 @@
 !      Then the only value which can be negative is the mass flux at the top of
 !      the plume. However, this value was set to zero a few lines above.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          
          labort_physic=.false.
@@ -177,25 +178,24 @@
             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)
+            DO k=nlay,1,-1
+               print *, 'fm ', fm(igout,k+1)
+               print *, 'entr,detr', entr(igout,k), detr(igout,k)
+            ENDDO
+            print *, 'fm-', fm(igout,1)
+            print *, '-------------------------------'
             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
@@ -212,16 +212,17 @@
             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)
+            DO k=nlay,1,-1
+               print *, 'fm ', fm(igout,k+1)
+               print *, 'entr,detr', entr(igout,k), detr(igout,k)
+            ENDDO
+            print *, 'fm-', fm(igout,1)
+            print *, '-------------------------------'
             abort_message = 'Negative entrainment in thermcell_flux'
             CALL abort_physic(modname,abort_message,1)
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test detrainment positif
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
          
          DO ig=1,ngrid
@@ -238,16 +239,17 @@
             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)
+            DO k=nlay,1,-1
+               print *, 'fm ', fm(igout,k+1)
+               print *, 'entr,detr', entr(igout,k), detr(igout,k)
+            ENDDO
+            print *, 'fm-', fm(igout,1)
+            print *, '-------------------------------'
             abort_message = 'Negative detrainment in thermcell_flux'
             CALL abort_physic(modname,abort_message,1)
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test sur fraca constante ou croissante au-dessus de lalim
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
          
 ! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0
@@ -267,7 +269,7 @@
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test sur flux de masse constant ou decroissant au-dessus de lalim
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
          
 ! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0
@@ -283,9 +285,9 @@
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test detrainment inferieur au flux de masse
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : Even if fm has no negative value, it can be lesser than detr.
 !      That's not suitable because when we will mix the plume with the
@@ -294,5 +296,5 @@
 !      otherwise we get a negative outgoing mass flux (cf. above).
 !      That is why we correct the detrainment as follows.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          
          DO ig=1,ngrid
@@ -314,5 +316,5 @@
                ncorecfm6  = ncorecfm6 + 1
                
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! Dans le cas ou on est au dessus de la couche d'alimentation et que le
 ! detrainement est plus fort que le flux de masse, on stope le thermique.
@@ -320,7 +322,7 @@
 !
 ! AB : Do we have to stop the plume here?
-! AB : Attention, if lmax is modified, be sure that fm, zw2, entr and detr are
-!      set to zero above this level!
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB : Attention, if lmax is modified, be sure that fm, zw2, entr, detr are set
+!      to zero above this level!
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !               IF (l.gt.lalim(ig)) THEN
 !                  lmax(ig)=l
@@ -331,5 +333,5 @@
 !               ENDIF
             ENDIF
-!            
+            
 !            IF (l.gt.lmax(ig)) THEN
 !               detr(ig,l) = 0.
@@ -360,9 +362,9 @@
 !         ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test fraction masse entrainee inferieure a fomass_max
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : Entrainment is bigger than the maximal authorized value.
 !      If we consider that the excess entrainement is in fact plume air which
@@ -374,5 +376,5 @@
 !      detrainment and mass flux profiles such as we do not exceed the maximal
 !      authorized entrained mass.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          
          labort_physic=.false.
@@ -432,5 +434,5 @@
             print *, '   ptimestep  :', ptimestep
             print *, 'norm  :', f(igout)
-            print *, 'alim* :', alim_star(igout,lout)
+!            print *, 'alim* :', alim_star(igout,lout)
             print *, 'entr* :', entr_star(igout,lout)
             print *, '--------------------'
@@ -447,9 +449,9 @@
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test fraction couverte inferieure a alphamax
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! FH : Partie a revisiter.
 !      Il semble qu'etaient codees ici deux optiques dans le cas F/(rho*w) > 1
@@ -461,8 +463,9 @@
 !      semble de toutes facons deraisonable d'avoir des fractions de 1...
 !      Ci-dessous, et dans l'etat actuel, le premier des deux if est sans doute inutile.  
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                 
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                 
          
          DO ig=1,ngrid
             zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alphamax
+            
             IF (fm(ig,l+1) .gt. zfm) THEN
                f_old = fm(ig,l+1)
@@ -470,10 +473,10 @@
                detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1)
                
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! 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 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)
@@ -498,12 +501,12 @@
       ENDDO
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! We check if fm, entr and detr vanish correctly in layer lmax
-!------------------------------------------------------------------------------
-      
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+      
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : test added to avoid problem when lmax = 0, which is not a fortran index.
 !      Theoretically, we always get lmax >= lmin >= linf > 0
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       DO ig=1,ngrid
          IF (lmax(ig).gt.0) THEN
@@ -514,7 +517,7 @@
       ENDDO
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Impression des bidouilles utilisees pour corriger les problemes
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       IF (prt_level.ge.1) THEN
@@ -553,7 +556,7 @@
       ENDIF
       
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : temporary test added to check the validity of eq. df/dz = e - d
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !      DO l=1,nlay
 !         DO ig=1,ngrid
