subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,  &
     &           masse,q,dq,qa,lev_out)
      USE print_control_mod, ONLY: prt_level
      implicit none

!=======================================================================
!
!   Calcul du transport verticale dans la couche limite en presence
!   de "thermiques" explicitement representes
!   calcul du dq/dt une fois qu'on connait les ascendances
!
! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
!  Introduction of an implicit computation of vertical advection in
!  the environment of thermal plumes in thermcell_dq
!  impl =     0 : explicit, 1 : implicit, -1 : old version
!
!=======================================================================

      integer ngrid,nlay,impl

      real ptimestep
      real masse(ngrid,nlay),fm(ngrid,nlay+1)
      real entr(ngrid,nlay)
      real q(ngrid,nlay)
      real dq(ngrid,nlay)
      integer lev_out                           ! niveau pour les print

      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)

      real zzm

      integer ig,k
      real cfl

      real qold(ngrid,nlay),fqa(ngrid,nlay+1)
      integer niter,iter
      CHARACTER (LEN=20) :: modname='thermcell_dq'
      CHARACTER (LEN=80) :: abort_message


! Old explicite scheme
      if (impl<=-1) then
         call thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,  &
     &           masse,q,dq,qa,lev_out)
         return
      endif

! Calcul du critere CFL pour l'advection dans la subsidence
      cfl = 0.
      do k=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,1',k,entr(ig,k)*ptimestep,masse(ig,k)
               abort_message = 'entr dt > m, 1st'
               CALL abort_physic (modname,abort_message,1)
            endif
         enddo
      enddo

      qold=q


      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)
!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

! Computation of tracer concentrations in the ascending plume
      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))*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*,'qa<0!!!'
            endif
            if (q(ig,k).lt.0.) then
!               print*,'q<0!!!'
            endif
         enddo
      enddo

! Plume vertical flux
      do k=2,nlay-1
         fqa(:,k)=fm(:,k)*qa(:,k-1)
      enddo
      fqa(:,1)=0. ; fqa(:,nlay)=0.


! Trace species evolution
   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)
      enddo
   else
      do k=nlay-1,1,-1
! FH debut de modif : le calcul ci dessous modifiait numériquement
! la concentration quand le flux de masse etait nul car on divisait
! puis multipliait par masse/ptimestep.
!        q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
!    &               /(fm(:,k)+masse(:,k)/ptimestep)
         q(:,k)=(q(:,k)+ptimestep/masse(:,k)*(fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1))) &
      &               /(1.+fm(:,k)*ptimestep/masse(:,k))
! FH fin de modif.
      enddo
   endif

! Tendencies
      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

#ifdef ISO      
      subroutine thermcell_dq_iso(ngrid,nlay,impl,ptimestep,fm,entr,  &
     &           masse,q,dq,qa,lev_out,xt,dxt,xta)
  USE infotrac_phy, ONLY : ntraciso
#ifdef ISOVERIF
  USE isotopes_mod, ONLY: iso_eau,iso_HDO
  USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
        iso_verif_aberrant_enc_vect2D,iso_verif_egalite_vect2D
#endif
      implicit none


#include "iniprint.h"
!=======================================================================
!
!   Calcul du transport verticale dans la couche limite en presence
!   de "thermiques" explicitement representes
!   calcul du dq/dt une fois qu'on connait les ascendances
!
! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
!  Introduction of an implicit computation of vertical advection in
!  the environment of thermal plumes in thermcell_dq
!  impl =     0 : explicit, 1 : implicit, -1 : old version
!
!=======================================================================

      integer ngrid,nlay,impl

      real ptimestep
      real masse(ngrid,nlay),fm(ngrid,nlay+1)
      real entr(ngrid,nlay)
      real q(ngrid,nlay)
      real dq(ngrid,nlay)
      integer lev_out                           ! niveau pour les print

      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)

      real zzm

      integer ig,k
      real cfl

      real qold(ngrid,nlay),fqa(ngrid,nlay+1)
      integer niter,iter
      CHARACTER (LEN=20) :: modname='thermcell_dq'
      CHARACTER (LEN=80) :: abort_message

! ifdef ISO
      real xt(ntraciso,ngrid,nlay)
      real dxt(ntraciso,ngrid,nlay)
      real xta(ntraciso,ngrid,nlay)
      real q_avant_entr(ngrid,nlay)
      integer ixt
      real wxtd(ntraciso,ngrid,nlay+1)
      real xtold(ntraciso,ngrid,nlay)
      real fxta(ntraciso,ngrid,nlay+1)
! endif

! Old explicite scheme
      if (impl==-1) then
         call thermcell_dq_o_iso(ngrid,nlay,ptimestep,fm,entr,  &
     &           masse,q,dq,qa,lev_out,xt,dxt,xta)
         return
      endif

! Calcul du critere CFL pour l'advection dans la subsidence
      cfl = 0.
      do k=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 ',entr(ig,k)*ptimestep,masse(ig,k)
               abort_message = ''
               CALL abort_gcm (modname,abort_message,1)
            endif
         enddo
      enddo

      qold=q
!#ifdef ISO
      xtold=xt
!#endif

      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)
!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
      ! bla90

! Computation of tracer concentrations in the ascending plume
      do ig=1,ngrid
         qa(ig,1)=q(ig,1)
      enddo
!#ifdef ISO
      do ig=1,ngrid
       do ixt=1,ntraciso
         xta(ixt,ig,1)=xt(ixt,ig,1)
       enddo
      enddo
!#endif      

      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
         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
     &         /(fm(ig,k+1)+detr(ig,k))
!#ifdef ISO
            do ixt=1,ntraciso
               xta(ixt,ig,k)=(fm(ig,k)*xta(ixt,ig,k-1) &
     &                  +entr(ig,k)*xt(ixt,ig,k))  &
     &                  /(fm(ig,k+1)+detr(ig,k))
             enddo
!#endif      
            else
               qa(ig,k)=q(ig,k)
!#ifdef ISO
               do ixt=1,ntraciso
                xta(ixt,ig,k)=xt(ixt,ig,k)
               enddo
!#endif      
            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

!#ifdef ISO
#ifdef ISOVERIF
      if (iso_HDO.gt.0) then
      call iso_verif_aberrant_enc_vect2D( &
     &           xta,qa, &
     &           'thermcell_dq_iso 320, qa',ntraciso,ngrid,nlay)
      endif        
#endif
!#endif
!
! Plume vertical flux
      do k=2,nlay-1
         fqa(:,k)=fm(:,k)*qa(:,k-1)
      enddo
      fqa(:,1)=0. ; fqa(:,nlay)=0.
!#ifdef ISO
      do k=2,nlay-1
        do ixt=1,ntraciso
         fxta(ixt,:,k)=fm(:,k)*xta(ixt,:,k-1)
        enddo
      enddo
      fxta(ixt,:,1)=0. ; fxta(ixt,:,nlay)=0.
!#endif


! Trace species evolution
   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)
!#ifdef ISO  
         do ixt=1,ntraciso
           xt(ixt,:,k)=xt(ixt,:,k)+(fxta(ixt,:,k)-fxta(ixt,:,k+1) &
     &               -fm(:,k)*xt(ixt,:,k)+fm(:,k+1)*xt(ixt,:,k+1)) &
     &               *ptimestep/masse(:,k)
         enddo ! do ixt=1,ntraciso
!#endif   
      enddo !do k=1,nlay-1
   else !if (impl==0) then
      do k=nlay-1,1,-1
         q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
     &               /(fm(:,k)+masse(:,k)/ptimestep)
!#ifdef ISO  
         do ixt=1,ntraciso
           xt(ixt,:,k)=xt(ixt,:,k)+(fxta(ixt,:,k)-fxta(ixt,:,k+1) &
     &               -fm(:,k)*xt(ixt,:,k)+fm(:,k+1)*xt(ixt,:,k+1)) &
     &               *ptimestep/masse(:,k)
         enddo ! do ixt=1,ntraciso
!#endif
      enddo !do k=nlay-1,1,-1
   endif !if (impl==0) then

! Tendencies
      do k=1,nlay
         do ig=1,ngrid
            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
            q(ig,k)=qold(ig,k)
!#ifdef ISO  
         do ixt=1,ntraciso
           dxt(ixt,ig,k)=(xt(ixt,ig,k)-xtold(ixt,ig,k))/ptimestep
           xt(ixt,ig,k)=xtold(ixt,ig,k)
         enddo ! do ixt=1,ntraciso
!#endif
         enddo
      enddo

#ifdef ISOVERIF
      if (iso_HDO.gt.0) then
      call iso_verif_aberrant_enc_vect2D( &
     &           xt,q, &
     &           'thermcell_dq_iso 219, q',ntraciso,ngrid,nlay)
      endif 
      if (iso_eau.gt.0) then  
      call iso_verif_egalite_vect2D(xt,q, &
     &           'thermcell_dq_iso 223, q',ntraciso,ngrid,nlay)
      call iso_verif_egalite_vect2D(dxt,dq, &
     &           'thermcell_dq_iso 224, dq',ntraciso,ngrid,nlay)
      endif 
#endif

return
end

#endif
! ifdef iso
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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

!=======================================================================
!
!   Calcul du transport verticale dans la couche limite en presence
!   de "thermiques" explicitement representes
!   calcul du dq/dt une fois qu'on connait les ascendances
!
!=======================================================================

      integer ngrid,nlay,impl

      real ptimestep
      real masse(ngrid,nlay),fm(ngrid,nlay+1)
      real entr(ngrid,nlay)
      real q(ngrid,nlay)
      real dq(ngrid,nlay)
      integer lev_out                           ! niveau pour les print

      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)

      real zzm

      integer ig,k
      real cfl

      real qold(ngrid,nlay)
      real ztimestep
      integer niter,iter
      CHARACTER (LEN=20) :: modname='thermcell_dq'
      CHARACTER (LEN=80) :: abort_message



! Calcul du critere CFL pour l'advection dans la subsidence
      cfl = 0.
      do k=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)
               abort_message = 'entr dt > m, 2nd'
               CALL abort_physic (modname,abort_message,1)
            endif
         enddo
      enddo

!IM 090508     print*,'CFL CFL CFL CFL ',cfl

#undef CFL
#ifdef CFL
! On subdivise le calcul en niter pas de temps.
      niter=int(cfl)+1
#else
      niter=1
#endif

      ztimestep=ptimestep/niter
      qold=q


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)
!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

! Calcul du flux subsident

      do k=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))
#else

#define plusqueun
#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
#else
            wqd(ig,k)=fm(ig,k)*q(ig,k)
#endif
#endif

            if (wqd(ig,k).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 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


#ifdef ISO      
      subroutine thermcell_dq_o_iso(ngrid,nlay,ptimestep,fm,entr,  &
     &           masse,q,dq,qa,lev_out,xt,dxt,xta)
#ifdef ISO
    use infotrac_phy, ONLY: ntraciso
#ifdef ISOVERIF
  USE isotopes_mod, ONLY: iso_eau,iso_HDO
  USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
        iso_verif_aberrant_enc_vect2D,iso_verif_egalite_vect2D
#endif
#endif
      implicit none

#include "iniprint.h"
!=======================================================================
!
!   Calcul du transport verticale dans la couche limite en presence
!   de "thermiques" explicitement representes
!   calcul du dq/dt une fois qu'on connait les ascendances
!
!=======================================================================

      integer ngrid,nlay

      real ptimestep
      real masse(ngrid,nlay),fm(ngrid,nlay+1)
      real entr(ngrid,nlay)
      real q(ngrid,nlay)
      real dq(ngrid,nlay)
      integer lev_out                           ! niveau pour les print

      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)

      real zzm

      integer ig,k
      real cfl

      real qold(ngrid,nlay)
      real ztimestep
      integer niter,iter
      CHARACTER (LEN=20) :: modname='thermcell_dq'
      CHARACTER (LEN=80) :: abort_message

! ifdef ISO
      real xt(ntraciso,ngrid,nlay)
      real dxt(ntraciso,ngrid,nlay)
      real xta(ntraciso,ngrid,nlay)
      real xtold(ntraciso,ngrid,nlay)
      real q_avant_entr(ngrid,nlay)
      integer ixt
      real wxtd(ntraciso,ngrid,nlay+1)
! endif

! Calcul du critere CFL pour l'advection dans la subsidence
      cfl = 0.
      do k=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 ',entr(ig,k)*ptimestep,masse(ig,k)
               abort_message = ''
               CALL abort_gcm (modname,abort_message,1)
            endif
         enddo
      enddo

!IM 090508     print*,'CFL CFL CFL CFL ',cfl

#undef CFL
#ifdef CFL
! On subdivise le calcul en niter pas de temps.
      niter=int(cfl)+1
#else
      niter=1
#endif

      ztimestep=ptimestep/niter
      qold=q
!#ifdef ISO
      xtold=xt
!#endif

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)
!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
!#ifdef ISO
      do ig=1,ngrid
       do ixt=1,ntraciso
         xta(ixt,ig,1)=xt(ixt,ig,1)
       enddo
      enddo
!#endif      

      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))
!#ifdef ISO
            do ixt=1,ntraciso
               xta(ixt,ig,k)=(fm(ig,k)*xta(ixt,ig,k-1) &
     &                  +entr(ig,k)*xt(ixt,ig,k))  &
     &                  /(fm(ig,k+1)+detr(ig,k))
             enddo
!#endif      
            else
               qa(ig,k)=q(ig,k)
!#ifdef ISO
               do ixt=1,ntraciso
                xta(ixt,ig,k)=xt(ixt,ig,k)
               enddo
!#endif      
            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

!#ifdef ISO
#ifdef ISOVERIF
      if (iso_HDO.gt.0) then
      call iso_verif_aberrant_enc_vect2D( &
     &           xta,qa, &
     &           'thermcell_dq_iso 320, qa',ntraciso,ngrid,nlay)
      endif        
#endif
!#endif

! Calcul du flux subsident

      do k=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))
!#ifdef ISO
! attention à l'advection: on suppose que R varie linéairement entre les
! couches, comme q.
             do ixt=1,ntraciso
               wxtd(ixt,ig,k)=wqd(ig,k)*0.5*( &
      &         xt(ixt,ig,k-1)/q(ig,k-1)+xt(ixt,ig,k)/q(ig,k))
             enddo !do ixt=1,ntraciso
!#endif
#else

#define plusqueun
#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)
!#ifdef ISO            
             do ixt=1,ntraciso
               wxtd(ixt,ig,k)=zzm*xt(ixt,ig,k)+(fm(ig,k)-zzm)*xt(ixt,ig,k+1)
             enddo !do ixt=1,ntraciso
!#endif             
            else
               wqd(ig,k)=fm(ig,k)*q(ig,k)
!#ifdef ISO            
             do ixt=1,ntraciso
               wxtd(ixt,ig,k)=fm(ig,k)*xt(ixt,ig,k)
             enddo !do ixt=1,ntraciso
!#endif             
            endif
#else
            wqd(ig,k)=fm(ig,k)*q(ig,k)
!#ifdef ISO            
             do ixt=1,ntraciso
               wxtd(ixt,ig,k)=fm(ig,k)*xt(ixt,ig,k)
             enddo !do ixt=1,ntraciso
!#endif             
#endif
#endif

            if (wqd(ig,k).lt.0.) then
!               print*,'wqd<0!!!'
            endif
         enddo
      enddo
      do ig=1,ngrid
         wqd(ig,1)=0.
         wqd(ig,nlay+1)=0.
!#ifdef ISO            
         do ixt=1,ntraciso
           wxtd(ixt,ig,1)=0.0
           wxtd(ixt,ig,nlay+1)=0.
         enddo !do ixt=1,ntraciso
!#endif          
      enddo
     

! Calcul des tendances
      do k=1,nlay
         do ig=1,ngrid
!#ifdef ISO         
            q_avant_entr(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)  &
     &               -wqd(ig,k)+wqd(ig,k+1))  &
     &               *ztimestep/masse(ig,k)
!#endif
            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

!#ifdef ISO
! Hypothèse: on ajoute le détrainement, et on fait la subsidence. Ensuite, on
! retranche la perte vers l'entrainement en supposant que c'est un puit non
! fractionnant
      do k=1,nlay
         do ig=1,ngrid
           do ixt=1,ntraciso
             xt(ixt,ig,k)=xt(ixt,ig,k)+(detr(ig,k)*xta(ixt,ig,k)  &
     &               -wxtd(ixt,ig,k)+wxtd(ixt,ig,k+1))  &
     &               *ztimestep/masse(ig,k)
             xt(ixt,ig,k)=xt(ixt,ig,k)/q_avant_entr(ig,k)*q(ig,k)
           enddo !do ixt=1,ntraciso
         enddo !do ig=1,ngrid
      enddo !do k=1,nlay
#ifdef ISOVERIF
      if (iso_eau.gt.0) then
      call iso_verif_egalite_vect2D( &
     &           xt,q, &
     &           'thermcell_dq_iso 421',ntraciso,ngrid,nlay)
      endif
      if (iso_HDO.gt.0) then
      call iso_verif_aberrant_enc_vect2D( &
     &           xt,q, &
     &           'thermcell_dq_iso 426',ntraciso,ngrid,nlay)
      endif      
#endif      
!#endif

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)
!#ifdef ISO
            do ixt=1,ntraciso
              dxt(ixt,ig,k)=(xt(ixt,ig,k)-xtold(ixt,ig,k))/ptimestep
              xt(ixt,ig,k)=xtold(ixt,ig,k)
            enddo !  do ixt=1,ntraciso
!#endif
         enddo
      enddo

      return
      end

#endif
! ifdef iso