subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm0,entr0,  &
     &    detr0,masse0,q_therm,dq_therm,charvar,zdz)
      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)
      CHARACTER (LEN=20), INTENT(IN) :: charvar
      REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
      REAL, INTENT(IN) :: zdz(ngridmx,nlayermx)

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

      REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx)  ! dq/dt -> derivative

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

      REAL q(ngridmx,nlayermx)
      REAL qa(ngridmx,nlayermx)
      REAL qd(ngridmx,nlayermx)
      INTEGER ig,k
      REAL gammaf(ngridmx,nlayermx)

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

      qa(:,:)=q_therm(:,:)
      q(:,:)=q_therm(:,:)

      do ig=1,ngridmx
         do k=1, nlayermx
            if (fm0(ig,k)+entr0(ig,k) .gt. 0.) then
              gammaf(ig,k)=fm0(ig,k)/(fm0(ig,k)+entr0(ig,k))
            else
              gammaf(ig,k)=0.
            endif
         enddo
      enddo


! =========== Updraft ============================================

!      qa(:,1)=q(:,1)

      do ig=1,ngridmx
         do k=2, nlayermx
 
             qa(ig,k)=gammaf(ig,k)*qa(ig,k-1) +(1.-gammaf(ig,k))*q(ig,k)

         enddo
      enddo



! ====== dq ======================================================

!      do ig=1,ngridmx
!         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)
!       enddo
!       do k=2,nlayermx-1
!         do ig=1, ngridmx
!         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)
!         enddo
!      enddo
!
!         do ig=1, ngridmx
!         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)
!         
!         enddo

        do ig=1, ngridmx
           do k=1,nlayermx-1
              dq_therm(ig,k)=-(1./masse0(ig,k))*(  &
     &           fm0(ig,k+1)*(qa(ig,k+1)-q(ig,k+1)) -   &
     &           fm0(ig,k)*(qa(ig,k)-q(ig,k))          ) &
     &       /zdz(ig,k)

           enddo
        enddo

      return
      end