subroutine thermcell_dqupdown(ngrid,nlayer,ptimestep,fm0,entr0,  &
     &    detr0,masse0,q_therm,dq_therm,ztvd,fm_down,ztv,charvar,lmax)
      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
!   Version modifiee pour prendre les downdrafts a la place de la 
!   subsidence compensatoire
!=======================================================================

#include "dimensions.h"
#include "dimphys.h"

! ============================ INPUTS ============================

      INTEGER, INTENT(IN) :: ngrid,nlayer
      REAL, INTENT(IN) :: ptimestep
      REAL, INTENT(IN) :: fm0(ngridmx,nlayermx+1)
      REAL, INTENT(IN) :: entr0(ngridmx,nlayermx),detr0(ngridmx,nlayermx)
      REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx)
      REAL, INTENT(IN) :: fm_down(ngridmx,nlayermx+1)
      REAL, INTENT(IN) :: ztvd(ngridmx,nlayermx)
      REAL, INTENT(IN) :: ztv(ngridmx,nlayermx)
      CHARACTER (LEN=20), INTENT(IN) :: charvar
      REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
      INTEGER, INTENT(IN) :: lmax(ngridmx)

! ============================ OUTPUTS ===========================

      REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx)

! ============================ LOCAL =============================

!      REAL detr0(ngridmx,nlayermx)
      REAL detrd(ngridmx,nlayermx)
      REAL entrd(ngridmx,nlayermx)      
      REAL fmd(ngridmx,nlayermx+1)
      REAL q(ngridmx,nlayermx)
      REAL qa(ngridmx,nlayermx)
      REAL qd(ngridmx,nlayermx)
      INTEGER ig,k
      LOGICAL active(ngridmx,nlayermx)
      INTEGER lmax_down(ngridmx),lmin_down(ngridmx)
      INTEGER ncorec

! =========== Init ==============================================

      entrd(:,:)=0.
      detrd(:,:)=0.
      qa(:,:)=q_therm(:,:)
      q(:,:)=q_therm(:,:)
      qd(:,:)=q_therm(:,:)
      active(:,:)=.false.

! previous calculation of zdthl_down uses the divergence of fmd
! so it can be negative without problem. Here we include the sign
! of fmd in the equations, so it has to be positive
!
      fmd(:,:)=-fm_down(:,:)
!
!! ========== Entrainment, Detrainement and Mass =================
!
!! ========== DOWNDRAFT TRANSPORT DISABLED FOR NOW ===============
!
!      do ig=1,ngridmx
!          if (ztv(ig,nlayermx)-ztvd(ig,nlayermx) .gt. 0.5) then
!            print*,"downdraft non nul derniere couche !!! (dqupdown)"
!          endif
!          detrd(ig,nlayermx)=0.
!          entrd(ig,nlayermx)=0.
!      enddo
!
!      do k=nlayermx-1,1,-1
!          do ig=1,ngridmx
!
!          if (ztv(ig,k)-ztvd(ig,k) .gt. 0.0001) then
!          detrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztv(ig,k)-ztvd(ig,k+1)))     &
!     &    /(ztv(ig,k)-ztvd(ig,k)) - fmd(ig,k))
!          entrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztvd(ig,k)-ztvd(ig,k+1)))    &
!     &    /(ztv(ig,k)-ztvd(ig,k)))     
!
!          endif
!        enddo
!      enddo
!
!! ======= We have computed entrainment and detrainment from a prescribed
!! mass flux and potential temp profile. Due to the way downdraft are parametrized,
!! this can yield negative entr and detr. We force it to be positive, but in 
!! order to conserve tracers, we need to recompute an adequate mass flux 
!! and modify interface rates, to preserve consistency.
!      lmax_down(:)=1
!      lmin_down(:)=1
!
!      do k=1,nlayermx
!        do ig=1,ngridmx
!        if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
!!         if (entrd(ig,k).gt.detrd(ig,k)) then
!        lmax_down(ig)=min(k,lmax(ig))
!        endif
!        enddo
!      enddo
!      do k=nlayermx,1,-1
!        do ig=1,ngridmx
!        if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
!!         if (detrd(ig,k).gt.entrd(ig,k)) then
!        lmin_down(ig)=k
!        endif
!        enddo
!      enddo
!
!
!      fmd(:,:)=0.
!       
!      do ig=1,ngridmx
!          if ((lmax_down(ig).gt.1) .and. ((lmax_down(ig)-lmin_down(ig)).gt.1)) then
!!          fmd(ig,lmax_down(ig))=0.
!!          entrd(ig,lmax_down(ig))=detrd(ig,lmax_down(ig))
!!          detrd(ig,lmax_down(ig))=0.
!!           print*,lmin_down(ig),lmax_down(ig),lmax(ig)
!
!            fmd(ig,lmax_down(ig)+1)=0.
!
!      do k=lmax_down(ig),lmin_down(ig)+1,-1
!            fmd(ig,k)=fmd(ig,k+1)+entrd(ig,k)-detrd(ig,k)
!      enddo
!            
!            fmd(ig,lmin_down(ig))=0.
!            detrd(ig,lmin_down(ig))=fmd(ig,lmin_down(ig)+1)+entrd(ig,lmin_down(ig))
!
!          else
!           entrd(ig,:)=0.
!           detrd(ig,:)=0.
!           active(ig,:)=.false.
!          endif
!
!      enddo
!         ncorec=0
!         do k=nlayermx,2,-1
!         do ig=1,ngridmx
!            if (fmd(ig,k).lt.0.) then
!!               detrd(ig,k)=max(0.,detrd(ig,k)+fmd(ig,k-1))
!!               fmd(ig,k-1)=0.
!!               entrd(ig,k-1)=0.
!!               detrd(ig,k-1)=0.
!!               lmin_down(ig)=k-1
!                fmd(ig,k)=fmd(ig,k+1)
!                detrd(ig,k)=entrd(ig,k)
!                ncorec=ncorec+1
!!                fmd(ig,k)=0.
!!                detrd(ig,k)=entrd(ig,k)+fmd(ig,k+1)
!
!            endif
!         enddo
!         enddo
!
!         if (ncorec .ne. 0) then
!         print*, 'corrections for negative downward mass flux :',ncorec
!         endif
!         print*, lmin_down(:),lmax_down(:)
!
!      do k=2,nlayermx
!      do ig=1,ngridmx
!          active(ig,k)=(k.ge.lmin_down(ig)).and.(k.le.lmax_down(ig)) &
!     & .and.(((fmd(ig,k)+detrd(ig,k))*ptimestep).gt.1.e-6*masse0(ig,k))
!      enddo
!      enddo
!
!
!      do ig=1,ngridmx
!      do k=lmin_down(ig),lmax_down(ig)
!          if(.not.active(ig,k)) then
!             active(ig,:)=.false.
!          endif
!      enddo
!      enddo
!
!      if(charvar .eq. 'tke') then
!      active(:,:)=.false.
!      endif
!
!!      do ig=1,ngridmx
!!         active(ig,lmax_down(ig))=(((fmd(ig,lmax_down(ig))+detrd(ig,lmax_down(ig)))* &
!!     &       ptimestep).gt.1.e-8*masse0(ig,lmax_down(ig)))
!!      enddo
!!
!!      do ig=1,ngridmx
!!          if (lmax_down(ig).gt.1) then
!!            do k=lmax_down(ig)-1,lmin_down(ig),-1
!!            active(ig,k)=(((fmd(ig,k)+detrd(ig,k))* &
!!     &       ptimestep).gt.1.e-8*masse0(ig,k)) &
!!     &     .and. active(ig,k+1)
!!            enddo
!!          else 
!!            active(ig,:)=.false.
!!          endif
!!      enddo
!!
!! ========== qa : q in updraft ==================================
!
      do k=2,nlayermx
         do ig=1,ngridmx
            if ((fm0(ig,k+1)+detr0(ig,k))*ptimestep.gt.  &
     &         1.e-5*masse0(ig,k)) then
         qa(ig,k)=(fm0(ig,k)*qa(ig,k-1)+entr0(ig,k)*q(ig,k))  &
     &         /(fm0(ig,k+1)+detr0(ig,k))

            if ((qa(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then
                 print*,'qa<0!!!',charvar,ig,k,fm0(ig,k),qa(ig,k-1),entr0(ig,k),q(ig,k),fm0(ig,k+1),detr0(ig,k)
                 print*,'---------> Cancelling qa'
                 qa(ig,k)=q(ig,k)
            endif
            else
               qa(ig,k)=q(ig,k)
            endif
            enddo
      enddo

! ========== qd : q in downdraft =================================
!
!      do k=nlayermx-1,1,-1
!         do ig=1,ngridmx
!            if (active(ig,k)) then
!         qd(ig,k)=(fmd(ig,k+1)*qd(ig,k+1)+entrd(ig,k)*q(ig,k))  &
!     &         /(fmd(ig,k)+detrd(ig,k))
!
!            if ((qd(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then
!     print*,'qd<0!!!',charvar,ig,k,fmd(ig,k),qd(ig,k),entrd(ig,k),q(ig,k),fmd(ig,k+1),detrd(ig,k),lmin_down(ig),lmax_down(ig)
!     print*, '---------> cancelling qd, no downdraft for this gridpoint'
!                     qd(ig,k)=q(ig,k)
!                     active(ig,:)=.false.
!            endif
!            else
!               qd(ig,k)=q(ig,k)
!            endif
!!          print*,'active,k,entr,detr,q,qd (down) :',active(ig,k),k,entrd(ig,k),detrd(ig,k),q(ig,k),qd(ig,k)
!         enddo
!      enddo
!
!
! ====== dq ======================================================

      do ig=1,ngridmx
         if(active(ig,1)) then

         dq_therm(ig,1)=(detr0(ig,1)*qa(ig,1)+detrd(ig,1)*qd(ig,1) &
      &               +fm0(ig,2)*q(ig,2)   &
      &               -entr0(ig,1)*q(ig,1)-entrd(ig,1)*q(ig,1)   &
      &               -fmd(ig,2)*q(ig,1)) &
      &               *ptimestep/masse0(ig,1)

         else
         dq_therm(ig,1)=(detr0(ig,1)*qa(ig,1)+fm0(ig,2)*q(ig,2) &
      &               -entr0(ig,1)*q(ig,1)) &
      &               *ptimestep/masse0(ig,1)

         endif
       enddo
      
       do k=2,nlayermx-1
         do ig=1, ngridmx

         if(active(ig,k)) then

         dq_therm(ig,k)=(detr0(ig,k)*qa(ig,k)+detrd(ig,k)*qd(ig,k) &
      &               +fm0(ig,k+1)*q(ig,k+1)+fmd(ig,k)*q(ig,k-1)   &
      &               -entr0(ig,k)*q(ig,k)-entrd(ig,k)*q(ig,k)   &
      &               -fm0(ig,k)*q(ig,k)-fmd(ig,k+1)*q(ig,k))      &
      &               *ptimestep/masse0(ig,k)


         else
         dq_therm(ig,k)=(detr0(ig,k)*qa(ig,k)+fm0(ig,k+1)*q(ig,k+1) &
      &               -entr0(ig,k)*q(ig,k)-fm0(ig,k)*q(ig,k))  &
      &               *ptimestep/masse0(ig,k)


         endif

         enddo
      enddo

         do ig=1, ngridmx

         if(active(ig,nlayermx)) then

         dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx)+detrd(ig,nlayermx)*qd(ig,nlayermx) &
      &               +fmd(ig,nlayermx)*q(ig,nlayermx-1)   &
      &         -entr0(ig,nlayermx)*q(ig,nlayermx)-entrd(ig,nlayermx)*q(ig,nlayermx)   &
      &               -fm0(ig,nlayermx)*q(ig,nlayermx)) &
      &               *ptimestep/masse0(ig,nlayermx)

         else
         dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx) &
      &             -entr0(ig,nlayermx)*q(ig,nlayermx)-fm0(ig,nlayermx)*q(ig,nlayermx)) &
      &               *ptimestep/masse0(ig,nlayermx)
         endif
         
         enddo
      return
      end