subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm0,entr0,  &
     &    alpha,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)
      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)
      REAL, INTENT(IN) :: alpha(ngridmx,nlayermx+1)

! ============================ 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),fmid(ngridmx,nlayermx)
      REAL gammae(ngridmx,nlayermx)
      REAL alphamid(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))
              gammae(ig,k)=fm0(ig,k+1)/(fm0(ig,k)+entr0(ig,k))
            else
              gammaf(ig,k)=0.
              gammae(ig,k)=1.
            endif
         enddo
      enddo


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

!      qa(:,1)=q_therm(:,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

! =========== Environnement =======================================

!     q(:,nlayermx)=q_therm(:,nlayermx)

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

         enddo
      enddo

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

      dq_therm(:,:)=0.
      fmid(:,:)=0.
      alphamid(:,:)=0.

        do ig=1, ngridmx
           do k=1,nlayermx
              fmid(ig,k) = 0.5*(fm0(ig,k)+fm0(ig,k+1))
              alphamid(ig,k) = 0.5*(alpha(ig,k)+alpha(ig,k+1))
           enddo
        enddo

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

           enddo
        enddo
      return
      end