Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 2101)
+++ /trunk/LMDZ.GENERIC/README	(revision 2102)
@@ -1455,2 +1455,5 @@
 - Remove useless variable zmax0 in thermcell_main, thermcell_height and physiq_mod.
 - Some minor changes in thermcell_plume and thermcell_main.
+
+==18/02/2018 == AB
+- use detr as thermcell_dq argument (called in thermcell_main) and clean up thermcell_dq
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90	(revision 2101)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90	(revision 2102)
@@ -2,11 +2,7 @@
 !
 !
-      SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,              &
-      &                       masse,q,dq,qa,lmin,lev_out)
-      
-      
-      USE print_control_mod, ONLY: prt_level
-      
-      IMPLICIT NONE
+SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,detr,masse,         &
+                        q,dq,qa,lmin,lev_out)
+      
       
 !==============================================================================
@@ -25,4 +21,9 @@
 !==============================================================================
       
+      USE print_control_mod, ONLY: prt_level
+      
+      IMPLICIT NONE
+      
+      
 !==============================================================================
 ! Declaration
@@ -32,14 +33,16 @@
 !      -------
       
-      INTEGER ngrid,nlay,impl
+      INTEGER ngrid, nlay
+      INTEGER impl
       INTEGER lmin(ngrid)
       INTEGER lev_out                           ! niveau pour les print
       
       REAL ptimestep
-      REAL masse(ngrid,nlay),fm(ngrid,nlay+1)
+      REAL masse(ngrid,nlay)
+      REAL fm(ngrid,nlay+1)
       REAL entr(ngrid,nlay)
+      REAL detr(ngrid,nlay)
       REAL q(ngrid,nlay)
       REAL qa(ngrid,nlay)
-      REAL detr(ngrid,nlay)
       
 !      outputs:
@@ -51,9 +54,10 @@
 !      ------
       
-      INTEGER ig,k
-      INTEGER niter,iter
+      INTEGER ig, l
+      INTEGER niter, iter
       
       REAL cfl
-      REAL qold(ngrid,nlay),fqa(ngrid,nlay+1)
+      REAL qold(ngrid,nlay)
+      REAL fqa(ngrid,nlay+1)
       REAL wqd(ngrid,nlay+1)
       REAL zzm
@@ -78,23 +82,25 @@
       
 !------------------------------------------------------------------------------
-! Calcul du critere CFL pour l'advection dans la subsidence
+! CFL criterion computation for advection in downdraft
 !------------------------------------------------------------------------------
       
       cfl = 0.
       
-      DO k=1,nlay
+      DO l=1,nlay
          DO ig=1,ngrid
-            zzm = masse(ig,k) / ptimestep
-            cfl = max(cfl, fm(ig,k) / zzm)
+            zzm = masse(ig,l) / ptimestep
+            cfl = max(cfl, fm(ig,l) / zzm)
             
-            IF (entr(ig,k).gt.zzm) THEN
+            IF (entr(ig,l).gt.zzm) THEN
                print *, 'ERROR: entrainment is greater than the layer mass!'
-               print *, 'entr*dt,mass', entr(ig,k)*ptimestep, masse(ig,k)
-               
-               print *, 'ig,l', ig, k
-               print *, 'fm-', fm(ig,k)
-               print *, 'entr,detr', entr(ig,k), detr(ig,k)
-               print *, 'fm+', fm(ig,k+1)
-               
+               print *, 'ig,l,entr', ig, l, entr(ig,l)
+               print *, '-------------------------------'
+               print *, 'entr*dt,mass', entr(ig,l)*ptimestep, masse(ig,l)
+               print *, '-------------------------------'
+               print *, 'fm+', fm(ig,l+1)
+               print *, 'entr,detr', entr(ig,l), detr(ig,l)
+               print *, 'fm ', fm(ig,l)
+               print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
+               print *, 'fm-', fm(ig,l-1)
                abort_message = 'entr dt > m, 1st'
                CALL abort_physic(modname,abort_message,1)
@@ -104,26 +110,30 @@
       
 !------------------------------------------------------------------------------
-! Detrainement computation
-!------------------------------------------------------------------------------
-      
-      DO k=1,nlay
-         DO ig=1,ngrid
-            detr(ig,k) = fm(ig,k) - fm(ig,k+1) + entr(ig,k)
+! Computation of tracer concentrations in the ascending plume
+!------------------------------------------------------------------------------
+      
+      DO ig=1,ngrid
+         DO l=1,lmin(ig)
+            qa(ig,l) = q(ig,l)
+         ENDDO
+      ENDDO
+      
+      DO ig=1,ngrid
+         DO l=lmin(ig)+1,nlay
+            if ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-5*masse(ig,l)) then
+               qa(ig,l) = (fm(ig,l) * qa(ig,l-1) + entr(ig,l) * q(ig,l))      &
+               &        / (fm(ig,l+1) + detr(ig,l))
+            else
+               qa(ig,l) = q(ig,l)
+            endif
             
-            IF (detr(ig,k).lt.0.) THEN
-               entr(ig,k) = entr(ig,k) - detr(ig,k)
-               detr(ig,k) = 0.
-!               print *, 'WARNING: detr2 is negative!'
-!               print *, 'detr', detr(ig,k)
-            ENDIF
-            
-!            IF (qa(ig,k).lt.0.) THEN
-!               print *, 'WARNING: fm2 is negative!'
-!               print *, 'fm', fm(ig,k)
+!            IF (qa(ig,l).lt.0.) THEN
+!               print *, 'WARNING: qa is negative!'
+!               print *, 'qa', qa(ig,l)
 !            ENDIF
             
-!            IF (q(ig,k).lt.0.) THEN
-!               print *, 'WARNING: entr2 is negative!'
-!               print *, 'entr', entr(ig,k)
+!            IF (q(ig,l).lt.0.) THEN
+!               print *, 'WARNING: q is negative!'
+!               print *, 'q', q(ig,l)
 !            ENDIF
          ENDDO
@@ -131,40 +141,9 @@
       
 !------------------------------------------------------------------------------
-! Computation of tracer concentrations in the ascending plume
-!------------------------------------------------------------------------------
-      
-      DO ig=1,ngrid
-         DO k=1,lmin(ig)
-            qa(ig,k) = q(ig,k)
-         ENDDO
-      ENDDO
-      
-      DO ig=1,ngrid
-         DO k=lmin(ig)+1,nlay
-            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then
-               qa(ig,k) = (fm(ig,k) * qa(ig,k-1) + entr(ig,k) * q(ig,k))      &
-               &        / (fm(ig,k+1) + detr(ig,k))
-            else
-               qa(ig,k) = q(ig,k)
-            endif
-            
-!            IF (qa(ig,k).lt.0.) THEN
-!               print *, 'WARNING: qa is negative!'
-!               print *, 'qa', qa(ig,k)
-!            ENDIF
-            
-!            IF (q(ig,k).lt.0.) THEN
-!               print *, 'WARNING: q is negative!'
-!               print *, 'q', q(ig,k)
-!            ENDIF
-         ENDDO
-      ENDDO
-      
-!------------------------------------------------------------------------------
 ! Plume vertical flux of q
 !------------------------------------------------------------------------------
       
-      DO k=2,nlay-1
-         fqa(:,k) = fm(:,k) * qa(:,k-1)
+      DO l=2,nlay-1
+         fqa(:,l) = fm(:,l) * qa(:,l-1)
       ENDDO
       
@@ -177,13 +156,13 @@
       
       IF (impl==0) THEN
-         do k=1,nlay-1
-            q(:,k) = q(:,k) + (fqa(:,k) - fqa(:,k+1) - fm(:,k) * q(:,k)       &
-            &      + fm(:,k+1) * q(:,k+1)) * ptimestep / masse(:,k)
+         do l=1,nlay-1
+            q(:,l) = q(:,l) + (fqa(:,l) - fqa(:,l+1) - fm(:,l) * q(:,l)       &
+            &      + fm(:,l+1) * q(:,l+1)) * ptimestep / masse(:,l)
          enddo
       ELSE
-         do k=nlay-1,1,-1
-            q(:,k) = ( q(:,k) + ptimestep / masse(:,k)                        &
-            &      * ( fqa(:,k) - fqa(:,k+1) + fm(:,k+1) * q(:,k+1) ) )       &
-            &      / ( 1. + fm(:,k) * ptimestep / masse(:,k) )
+         do l=nlay-1,1,-1
+            q(:,l) = ( q(:,l) + ptimestep / masse(:,l)                        &
+            &      * ( fqa(:,l) - fqa(:,l+1) + fm(:,l+1) * q(:,l+1) ) )       &
+            &      / ( 1. + fm(:,l) * ptimestep / masse(:,l) )
          ENDDO
       ENDIF
@@ -193,8 +172,8 @@
 !==============================================================================
       
-      DO k=1,nlay
+      DO l=1,nlay
          DO ig=1,ngrid
-            dq(ig,k) = (q(ig,k) - qold(ig,k)) / ptimestep
-            q(ig,k) = qold(ig,k)
+            dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep
+            q(ig,l) = qold(ig,l)
          ENDDO
       ENDDO
@@ -205,24 +184,17 @@
 
 
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-
-! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
-
-
-      subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,  &
-      &           masse,q,dq,qa,lev_out)
-      
-      USE print_control_mod, ONLY: prt_level
-      
-      implicit none
-
-!=======================================================================
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+Subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,masse,            &
+                          q,dq,qa,lev_out)
+      
+      
+!==============================================================================
 !
 !   Calcul du transport verticale dans la couche limite en presence
@@ -230,6 +202,15 @@
 !   calcul du dq/dt une fois qu'on connait les ascendances
 !
-!=======================================================================
-
+!==============================================================================
+      
+      USE print_control_mod, ONLY: prt_level
+      
+      implicit none
+      
+      
+!==============================================================================
+! Declaration
+!==============================================================================
+      
       integer ngrid,nlay,impl
 
@@ -245,5 +226,5 @@
       real zzm
 
-      integer ig,k
+      integer ig,l
       real cfl
 
@@ -258,10 +239,10 @@
 ! Calcul du critere CFL pour l'advection dans la subsidence
       cfl = 0.
-      do k=1,nlay
+      do l=1,nlay
          do ig=1,ngrid
-            zzm=masse(ig,k)/ptimestep
-            cfl=max(cfl,fm(ig,k)/zzm)
-            if (entr(ig,k).gt.zzm) then
-               print*,'entr*dt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
+            zzm=masse(ig,l)/ptimestep
+            cfl=max(cfl,fm(ig,l)/zzm)
+            if (entr(ig,l).gt.zzm) then
+               print*,'entr*dt>m,2',l,entr(ig,l)*ptimestep,masse(ig,l)
                abort_message = 'entr dt > m, 2nd'
                CALL abort_physic (modname,abort_message,1)
@@ -269,7 +250,7 @@
          enddo
       enddo
-
+      
 !IM 090508     print*,'CFL CFL CFL CFL ',cfl
-
+      
 #undef CFL
 #ifdef CFL
@@ -279,63 +260,65 @@
       niter=1
 #endif
-
+      
       ztimestep=ptimestep/niter
       qold=q
-
-
-do iter=1,niter
-      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
-
+      
+      do iter=1,niter
+         if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
+         
 !   calcul du detrainement
-      do k=1,nlay
-         do ig=1,ngrid
-            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
-!           print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
+         do l=1,nlay
+            do ig=1,ngrid
+               detr(ig,l)=fm(ig,l)-fm(ig,l+1)+entr(ig,l)
+!               print*,'Q2 DQ ',detr(ig,l),fm(ig,l),entr(ig,l)
 !test
-            if (detr(ig,k).lt.0.) then
-               entr(ig,k)=entr(ig,k)-detr(ig,k)
-               detr(ig,k)=0.
-!               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
-!     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
-            endif
-            if (fm(ig,k+1).lt.0.) then
-!               print*,'fm2<0!!!'
-            endif
-            if (entr(ig,k).lt.0.) then
-!               print*,'entr2<0!!!'
-            endif
-         enddo
-      enddo
-
-!   calcul de la valeur dans les ascendances
-      do ig=1,ngrid
-         qa(ig,1)=q(ig,1)
-      enddo
-
-      do k=2,nlay
-         do ig=1,ngrid
-            if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt.  &
-     &         1.e-5*masse(ig,k)) then
-         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
-     &         /(fm(ig,k+1)+detr(ig,k))
-            else
-               qa(ig,k)=q(ig,k)
-            endif
-            if (qa(ig,k).lt.0.) then
-!               print*,'qa<0!!!'
-            endif
-            if (q(ig,k).lt.0.) then
-!               print*,'q<0!!!'
-            endif
-         enddo
-      enddo
-
+               if (detr(ig,l).lt.0.) then
+                  entr(ig,l)=entr(ig,l)-detr(ig,l)
+                  detr(ig,l)=0.
+!                  print*,'detr2<0!!!','ig=',ig,'l=',l,'f=',fm(ig,l),
+!     s            'f+1=',fm(ig,l+1),'e=',entr(ig,l),'d=',detr(ig,l)
+               endif
+               
+!               if (fm(ig,l+1).lt.0.) then
+!                  print*,'fm2<0!!!'
+!               endif
+               
+!               if (entr(ig,l).lt.0.) then
+!                  print*,'entr2<0!!!'
+!               endif
+            enddo
+         enddo
+         
+! calcul de la valeur dans les ascendances
+        do ig=1,ngrid
+            qa(ig,1)=q(ig,1)
+         enddo
+         
+         do l=2,nlay
+            do ig=1,ngrid
+               if ((fm(ig,l+1)+detr(ig,l))*ztimestep.gt.1.e-5*masse(ig,l)) then
+                  qa(ig,l) = (fm(ig,l)*qa(ig,l-1)+entr(ig,l)*q(ig,l))         &
+                  &        / (fm(ig,l+1)+detr(ig,l))
+               else
+                  qa(ig,l)=q(ig,l)
+               endif
+               
+               if (qa(ig,l).lt.0.) then
+!                  print*,'qa<0!!!'
+               endif
+               
+               if (q(ig,l).lt.0.) then
+!                  print*,'q<0!!!'
+               endif
+            enddo
+         enddo
+         
 ! Calcul du flux subsident
-
-      do k=2,nlay
-         do ig=1,ngrid
+         
+         do l=2,nlay
+            do ig=1,ngrid
 #undef centre
 #ifdef centre
-             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
+                wqd(ig,l)=fm(ig,l)*0.5*(q(ig,l-1)+q(ig,l))
 #else
 
@@ -343,50 +326,48 @@
 #ifdef plusqueun
 ! Schema avec advection sur plus qu'une maille.
-            zzm=masse(ig,k)/ztimestep
-            if (fm(ig,k)>zzm) then
-               wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
-            else
-               wqd(ig,k)=fm(ig,k)*q(ig,k)
-            endif
+               zzm=masse(ig,l)/ztimestep
+               if (fm(ig,l)>zzm) then
+                 wqd(ig,l)=zzm*q(ig,l)+(fm(ig,l)-zzm)*q(ig,l+1)
+               else
+                  wqd(ig,l)=fm(ig,l)*q(ig,l)
+               endif
 #else
-            wqd(ig,k)=fm(ig,k)*q(ig,k)
+               wqd(ig,l)=fm(ig,l)*q(ig,l)
 #endif
 #endif
 
-            if (wqd(ig,k).lt.0.) then
-!               print*,'wqd<0!!!'
-            endif
+               if (wqd(ig,l).lt.0.) then
+!                  print*,'wqd<0!!!'
+               endif
+            enddo
+         enddo
+         
+         do ig=1,ngrid
+            wqd(ig,1)=0.
+            wqd(ig,nlay+1)=0.
+         enddo
+         
+! Calcul des tendances
+         do l=1,nlay
+            do ig=1,ngrid
+               q(ig,l)=q(ig,l)+(detr(ig,l)*qa(ig,l)-entr(ig,l)*q(ig,l)           &
+               &        -wqd(ig,l)+wqd(ig,l+1))                                  &
+               &        *ztimestep/masse(ig,l)
+!               if (dq(ig,l).lt.0.) then
+!                  print*,'dq<0!!!'
+!               endif
+           enddo
          enddo
       enddo
-      do ig=1,ngrid
-         wqd(ig,1)=0.
-         wqd(ig,nlay+1)=0.
+      
+! Calcul des tendances
+      do l=1,nlay
+         do ig=1,ngrid
+            dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep
+            q(ig,l) = qold(ig,l)
+         enddo
       enddo
-     
-
-! Calcul des tendances
-      do k=1,nlay
-         do ig=1,ngrid
-            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
-     &               -wqd(ig,k)+wqd(ig,k+1))  &
-     &               *ztimestep/masse(ig,k)
-!            if (dq(ig,k).lt.0.) then
-!               print*,'dq<0!!!'
-!            endif
-         enddo
-      enddo
-
-
-enddo
-
-
-! Calcul des tendances
-      do k=1,nlay
-         do ig=1,ngrid
-            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
-            q(ig,k)=qold(ig,k)
-         enddo
-      enddo
-
-      return
-      end
+      
+      
+return
+end
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90	(revision 2101)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90	(revision 2102)
@@ -2,21 +2,21 @@
 !
 !
-SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep                           &
-                         ,pplay,pplev,pphi,firstcall                          &
-                         ,pu,pv,pt,po                                         &
-                         ,pduadj,pdvadj,pdtadj,pdoadj                         &
-                         ,f0,fm0,entr0,detr0                                  &
-                         ,zqta,zqla,ztv,ztva,ztla,zthl,zqsatth                &
-                         ,zw2,fraca                                           &
-                         ,lmin,lmix,lalim,lmax                                &
-                         ,zpspsk,ratqscth,ratqsdiff                           &
-                         ,Ale_bl,Alp_bl,lalim_conv,wght_th                    &
+SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep,                          &
+                          pplay,pplev,pphi,firstcall,                         &
+                          pu,pv,pt,po,                                        &
+                          pduadj,pdvadj,pdtadj,pdoadj,                        &
+                          f0,fm0,entr0,detr0,                                 &
+                          zqta,zqla,ztv,ztva,ztla,zthl,zqsatth,               &
+                          zw2,fraca,                                          &
+                          lmin,lmix,lalim,lmax,                               &
+                          zpspsk,ratqscth,ratqsdiff,                          &
+                          Ale_bl,Alp_bl,lalim_conv,wght_th,                   &
 !!! nrlmd le 10/04/2012
-                         ,pbl_tke,pctsrf,omega,airephy                        &
-                         ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0   &
-                         ,n2,s2,ale_bl_stat                                   &
-                         ,therm_tke_max,env_tke_max                           &
-                         ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke          &
-                         ,alp_bl_conv,alp_bl_stat)
+                          pbl_tke,pctsrf,omega,airephy,                       &
+                          zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,  &
+                          n2,s2,ale_bl_stat,                                  &
+                          therm_tke_max,env_tke_max,                          &
+                          alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,         &
+                          alp_bl_conv,alp_bl_stat)
 !!! fin nrlmd le 10/04/2012
       
@@ -213,9 +213,15 @@
 !------Sorties
       real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid)
-      real therm_tke_max0(ngrid),env_tke_max0(ngrid)
+      real therm_tke_max0(ngrid)
+      real env_tke_max0(ngrid)
       real n2(ngrid),s2(ngrid)
       real ale_bl_stat(ngrid)
-      real therm_tke_max(ngrid,nlay),env_tke_max(ngrid,nlay)
-      real alp_bl_det(ngrid),alp_bl_fluct_m(ngrid),alp_bl_fluct_tke(ngrid),alp_bl_conv(ngrid),alp_bl_stat(ngrid)
+      real therm_tke_max(ngrid,nlay)
+      real env_tke_max(ngrid,nlay)
+      real alp_bl_det(ngrid)
+      real alp_bl_fluct_m(ngrid)
+      real alp_bl_fluct_tke(ngrid)
+      real alp_bl_conv(ngrid)
+      real alp_bl_stat(ngrid)
 !------Local
       integer nsrf
@@ -225,16 +231,11 @@
       real interp(ngrid)                        ! Coef d'interpolation pour le LCL
 !------Triggering
-      real Su                                   ! Surface unite: celle d'un updraft elementaire
-      parameter(Su=4e4)
-      real hcoef                                ! Coefficient directeur pour le calcul de s2
-      parameter(hcoef=1)
-      real hmincoef                             ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2 
-      parameter(hmincoef=0.3)
-      real eps1                                 ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd)
-      parameter(eps1=0.3)
+      real,parameter :: Su=4e4                  ! Surface unite: celle d'un updraft elementaire
+      real,parameter :: hcoef                   ! Coefficient directeur pour le calcul de s2
+      real,parameter :: hmincoef=0.3            ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2 (hmincoef=0.3)
+      real,parameter :: eps1=0.3                ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd)
       real hmin(ngrid)                          ! Ordonnee a l'origine pour le calcul de s2
       real zmax_moy(ngrid)                      ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
-      real zmax_moy_coef
-      parameter(zmax_moy_coef=0.33)
+      real,parameter :: zmax_moy_coef=0.33
       real depth(ngrid)                         ! Epaisseur moyenne du cumulus
       real w_max(ngrid)                         ! Vitesse max statistique 
@@ -244,20 +245,19 @@
       real pbl_tke_max0(ngrid)                  ! TKE moyenne au LCL
       real w_ls(ngrid,nlay)                     ! Vitesse verticale grande echelle (m/s)
-      real coef_m                               ! On considere un rendement pour alp_bl_fluct_m
-      parameter(coef_m=1.)
-      real coef_tke                             ! On considere un rendement pour alp_bl_fluct_tke
-      parameter(coef_tke=1.)
+      real,parameter :: coef_m=1.               ! On considere un rendement pour alp_bl_fluct_m
+      real,parameter :: coef_tke=1.             ! On considere un rendement pour alp_bl_fluct_tke
       
 !!! fin nrlmd le 10/04/2012
       
 ! Nouvelles variables pour la convection
+      integer lalim_conv(ngrid)
+      integer n_int(ngrid)
       real Ale_bl(ngrid)
       real Alp_bl(ngrid)
-      real alp_int(ngrid),dp_int(ngrid),zdp
+      real alp_int(ngrid)
+      real dp_int(ngrid),zdp
       real ale_int(ngrid)
-      integer n_int(ngrid)
       real fm_tot(ngrid)
       real wght_th(ngrid,nlay)
-      integer lalim_conv(ngrid)
       
       CHARACTER*2 str2
@@ -515,5 +515,5 @@
       IF (tau_thermals>1.) THEN
          lambda = exp(-ptimestep/tau_thermals)
-         f0 = (1.-lambda) * f + lambda * f0
+         f0(:) = (1.-lambda) * f(:) + lambda * f0(:)
       ELSE
          f0(:) = f(:)
@@ -527,6 +527,4 @@
          print *, 'f0 =', f0(1)
          CALL abort_physic(modname,abort_message,1)
-      ELSE
-         print *, 'f0 =', f0(1)
       ENDIF
       
@@ -587,8 +585,8 @@
 !------------------------------------------------------------------------------
       
-      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
+      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse,    &
       &                 zthl,zdthladj,zta,lmin,lev_out)
       
-      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
+      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse,    &
       &                 po,pdoadj,zoa,lmin,lev_out)
       
@@ -614,8 +612,8 @@
 ! Calcul purement conservatif pour le transport de V=(u,v)
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,       &
+         CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
          &                 zu,pduadj,zua,lmin,lev_out)
          
-         CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,       &
+         CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
          &                 zv,pdvadj,zva,lmin,lev_out)
       ENDIF
@@ -631,51 +629,51 @@
 !------------------------------------------------------------------------------
          
-         if (prt_level.ge.1) print*,'14a OK convect8'
-         
-         do ig=1,ngrid
+         DO ig=1,ngrid
             nivcon(ig) = 0
             zcon(ig) = 0.
-         enddo 
+         ENDDO 
 !nouveau calcul
          do ig=1,ngrid
             CHI = zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
             pcon(ig) = pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
-         enddo
+         ENDDO
          
 !IM       do k=1,nlay
          do k=1,nlay-1
             do ig=1,ngrid
-               if ((pcon(ig).le.pplay(ig,k)).and.(pcon(ig).gt.pplay(ig,k+1))) then
-                  zcon2(ig) = zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
-               endif
-            enddo
-         enddo
+               IF ((pcon(ig).le.pplay(ig,k)).and.(pcon(ig).gt.pplay(ig,k+1))) then
+                  zcon2(ig) = zlay(ig,k) - (pcon(ig) - pplay(ig,k))           &
+                            / (RG * rho(ig,k)) / 100.
+               ENDIF
+            ENDDO
+         ENDDO
          
          ierr = 0
          
          do ig=1,ngrid
-           if (pcon(ig).le.pplay(ig,nlay)) then 
-              zcon2(ig) = zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
+           IF (pcon(ig).le.pplay(ig,nlay)) then 
+              zcon2(ig) = zlay(ig,nlay) - (pcon(ig) - pplay(ig,nlay))         &
+                        / (RG * rho(ig,nlay)) / 100.
               ierr = 1
-            endif
-         enddo
-         
-         if (ierr==1) then
+            ENDIF
+         ENDDO
+         
+         IF (ierr==1) then
               write(*,*) 'ERROR: thermal plumes rise the model top'
               CALL abort
-         endif
-         
-         if (prt_level.ge.1) print*,'14b OK convect8'
+         ENDIF
+         
+         IF (prt_level.ge.1) print*,'14b OK convect8'
          
          do k=nlay,1,-1
             do ig=1,ngrid
-               if (zqla(ig,k).gt.1e-10) then
+               IF (zqla(ig,k).gt.1e-10) then
                   nivcon(ig) = k
                   zcon(ig) = zlev(ig,k)
-               endif
-            enddo
-         enddo
-         
-         if (prt_level.ge.1) print*,'14c OK convect8'
+               ENDIF
+            ENDDO
+         ENDDO
+         
+         IF (prt_level.ge.1) print*,'14c OK convect8'
          
 !------------------------------------------------------------------------------
@@ -690,8 +688,8 @@
                ratqscth(ig,l) = 0.
                ratqsdiff(ig,l) = 0.
-            enddo
-         enddo
-         
-         if (prt_level.ge.1) print*,'14d OK convect8'
+            ENDDO
+         ENDDO
+         
+         IF (prt_level.ge.1) print*,'14d OK convect8'
          
          do l=1,nlay
@@ -701,17 +699,17 @@
                thetath2(ig,l) = zf2*(ztla(ig,l)-zthl(ig,l))**2
                
-               if (zw2(ig,l).gt.1.e-10) then
+               IF (zw2(ig,l).gt.1.e-10) then
                   wth2(ig,l) = zf2*(zw2(ig,l))**2
                else
                   wth2(ig,l) = 0.
-               endif
+               ENDIF
                
-               wth3(ig,l) = zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))               &
+               wth3(ig,l) = zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))            &
                &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
                q2(ig,l) = zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
 !test: on calcul q2/po=ratqsc
                ratqscth(ig,l) = sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
-            enddo
-         enddo
+            ENDDO
+         ENDDO
          
 !------------------------------------------------------------------------------
@@ -724,18 +722,18 @@
                wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
                wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
-            enddo
-         enddo
-         
-         CALL thermcell_alp(ngrid,nlay,ptimestep,                                &
-         &                  pplay,pplev,                                         &
-         &                  fm0,entr0,lmax,                                      &
-         &                  Ale_bl,Alp_bl,lalim_conv,wght_th,                    &
-         &                  zw2,fraca,                                           &
-         &                  pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,  &
-         &                  pbl_tke,pctsrf,omega,airephy,                        &
-         &                  zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,   &
-         &                  n2,s2,ale_bl_stat,                                   &
-         &                  therm_tke_max,env_tke_max,                           &
-         &                  alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,          &
+            ENDDO
+         ENDDO
+         
+         CALL thermcell_alp(ngrid,nlay,ptimestep,                             &
+         &                  pplay,pplev,                                      &
+         &                  fm0,entr0,lmax,                                   &
+         &                  Ale_bl,Alp_bl,lalim_conv,wght_th,                 &
+         &                  zw2,fraca,pcon,                                   &
+         &                  rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,    &
+         &                  pbl_tke,pctsrf,omega,airephy,                     &
+         &                  zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,&
+         &                  n2,s2,ale_bl_stat,                                &
+         &                  therm_tke_max,env_tke_max,                        &
+         &                  alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,       &
          &                  alp_bl_conv,alp_bl_stat)
          
@@ -756,5 +754,5 @@
          ENDDO
          
-         if (prt_level.ge.1) print*,'14f OK convect8'
+         IF (prt_level.ge.1) print*,'14f OK convect8'
       
          DO l=1,nlay
@@ -763,5 +761,6 @@
                   zf  = fraca(ig,l)
                   zf2 = zf / (1.-zf)
-                  vardiff = vardiff + alim_star(ig,l) * (zqta(ig,l) * 1000. - var)**2
+                  vardiff = vardiff + alim_star(ig,l)                         &
+                          * (zqta(ig,l) * 1000. - var)**2
                ENDIF
             ENDDO
@@ -840,13 +839,13 @@
 !  test sur la hauteur des thermiques ...
       do i=1,ngrid
-!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
-        if (prt_level.ge.10) then
+!IMtemp IF (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
+        IF (prt_level.ge.10) then
             print *, 'WARNING ',comment,' au point ',i,' K= ',long(i)
             print *, '  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
             do k=1,nlay
                write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
-            enddo
-        endif
-      enddo
+            ENDDO
+        ENDIF
+      ENDDO
       
       
@@ -911,5 +910,5 @@
          detr0(:,k) = fm0(:,k) - fm0(:,k+1) + entr0(:,k)
          masse0(:,k) = (pplev(:,k) - pplev(:,k+1)) / RG
-      enddo
+      ENDDO
       
 !------------------------------------------------------------------------------
@@ -927,5 +926,5 @@
          detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
          fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
-      enddo
+      ENDDO
       
       fm(:,nlay+1)=0.
@@ -935,5 +934,5 @@
 !      do ig=1,ngrid
 !         qa(ig,1) = q(ig,1)
-!      enddo
+!      ENDDO
       
       q(:,:)=therm_tke_max(:,:)
@@ -941,20 +940,20 @@
       do ig=1,ngrid
          qa(ig,1)=q(ig,1)
-      enddo
-      
-      if (1==1) then
+      ENDDO
+      
+      IF (1==1) then
          do k=2,nlay
             do ig=1,ngrid
-               if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then
+               IF ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then
                   qa(ig,k) = (fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))         &
                   &        / (fm(ig,k+1)+detr(ig,k))
                else
                   qa(ig,k)=q(ig,k)
-               endif
+               ENDIF
                
-!               if (qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!'
-!               if (q(ig,k).lt.0.) print *, 'WARNING: q is negative!'
-            enddo
-         enddo
+!               IF (qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!'
+!               IF (q(ig,k).lt.0.) print *, 'WARNING: q is negative!'
+            ENDDO
+         ENDDO
          
 !------------------------------------------------------------------------------
@@ -965,12 +964,12 @@
             do ig=1,ngrid
                wqd(ig,k) = fm(ig,k) * q(ig,k)
-               if (wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!'
-            enddo
-         enddo
+               IF (wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!'
+            ENDDO
+         ENDDO
          
          do ig=1,ngrid
             wqd(ig,1) = 0.
             wqd(ig,nlay+1) = 0.
-         enddo
+         ENDDO
          
 !------------------------------------------------------------------------------
@@ -983,7 +982,7 @@
                &       * (detr(ig,k) * qa(ig,k) - entr(ig,k) * q(ig,k)        &
                &       - wqd(ig,k) + wqd(ig,k+1))
-            enddo
-         enddo
-      endif
+            ENDDO
+         ENDDO
+      ENDIF
       
       therm_tke_max(:,:) = q(:,:)
