subroutine thermcell_dq(ngrid,nlay,ptimestep,fm,entr,  &
     &           masse,q,dq,qa,lev_out)
      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

      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)

      integer ig,k

!   calcul du detrainement

      if (lev_out.ge.1) print*,'Q2 THERMCEL_DQ 0'

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

      do k=2,nlay
         do ig=1,ngrid
!             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
            wqd(ig,k)=fm(ig,k)*q(ig,k)
            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
     
      do k=1,nlay
         do ig=1,ngrid
            dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
     &               -wqd(ig,k)+wqd(ig,k+1))  &
     &               /masse(ig,k)
!            if (dq(ig,k).lt.0.) then
!               print*,'dq<0!!!'
!            endif
         enddo
      enddo

      return
      end