SUBROUTINE thermcell_updown_dq(ngrid,nlay,lmax,eup,dup,edn,ddn,masse,theta,dtheta,thetau,thetad)

!--------------------------------------------------------------
! thermcell_updown_dq: calcul du transport d'un traceur en présence 
! d'up/down drafts
!--------------------------------------------------------------

   ! Suite du travail :
   ! Calculer les tendances d'un traceur (ici theta) en tenant compte
   ! des up et down drafts et de la subsidence compensatoire.


   IMPLICIT NONE

! arguments

   integer,intent(in) :: ngrid,nlay
   real,intent(in), dimension(ngrid,nlay) :: eup,dup,edn,ddn,masse
   real,intent(in), dimension(ngrid,nlay) :: theta
   real,intent(out), dimension(ngrid,nlay) :: thetau,thetad,dtheta
   integer, intent(out), dimension(ngrid) :: lmax

   
! Local

   real, dimension(ngrid,nlay+1) :: fup,fdn

   integer ig,ilay

   fdn(:,:)=0.
   thetad(:,:)=0.

   ! lmax : indice tel que fu(kmax+1)=0
   
   ! Dans ce cas, pas besoin d'initialiser thetad(lmax) ( =theta(lmax) )

   do ilay=nlay,1,-1
      do ig=1,ngrid
         if (ilay.le.lmax(ig)) then
            fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay)
            thetad(ig,ilay)=( fdn(ig,ilay+1)*thetad(ig,ilay+1) + edn(ig,ilay)*theta(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay))
         endif
      enddo 
   enddo



!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
! Initialisations :
!------------------


!
 RETURN
   END

!=========================================================================

   SUBROUTINE thermcell_down(ngrid,nlay,po,pt,pu,pv,pplay,pplev,  &
     &           lmax,fup,eup,dup,theta)

!--------------------------------------------------------------
!thermcell_down: calcul des propriétés du panache descendant.
!--------------------------------------------------------------


   USE thermcell_ini_mod, ONLY : prt_level,RLvCp,RKAPPA,RETV
   IMPLICIT NONE

! arguments

   integer,intent(in) :: ngrid,nlay
   real,intent(in), dimension(ngrid,nlay) :: po,pt,pu,pv,pplay,eup,dup
   real,intent(in), dimension(ngrid,nlay) :: theta
   real,intent(in), dimension(ngrid,nlay+1) :: pplev,fup
   integer, intent(in), dimension(ngrid) :: lmax


   
! Local

   real, dimension(ngrid,nlay) :: edn,ddn,thetad
   real, dimension(ngrid,nlay+1) :: fdn

   integer ig,ilay
   real dqsat_dT
   logical mask(ngrid,nlay)

   edn(:,:)=0.
   ddn(:,:)=0.
   fdn(:,:)=0.
   thetad(:,:)=0.

   ! lmax : indice tel que fu(kmax+1)=0
   
   ! Dans ce cas, pas besoin d'initialiser thetad(lmax) ( =theta(lmax) )

! FH MODIFS APRES REUNIONS POUR COMMISSIONS
! quelques erreurs de declaration
! probleme si lmax=1 ce qui a l'air d'être le cas en début de simu. Devrait être 0 ?
! Remarques :
! on pourrait écrire la formule de thetad
!    www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
!    thetad(ig,ilay)= www * thetad(ig,ilay+1) + (1.-www) * theta(ig,ilay) 
! Elle a l'avantage de bien montré la conservation, l'idée fondamentale dans le 
!   transport qu'on ne fait que sommer des "sources" au travers d'un "propagateur"
!   (Green)
! Elle montre aussi beaucoup plus clairement pourquoi on n'a pas à se souccier (trop)
!   de la possible nulité du dénominateur


   do ilay=nlay,1,-1
      do ig=1,ngrid
         if (ilay.le.lmax(ig).and.lmax(ig)>1) then
            edn(ig,ilay)=0.5*dup(ig,ilay)
            ddn(ig,ilay)=0.5*eup(ig,ilay)
            fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay)
            thetad(ig,ilay)=( fdn(ig,ilay+1)*thetad(ig,ilay+1) + edn(ig,ilay)*theta(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay))
         endif
      enddo 
   enddo

   ! Suite du travail :
   ! Ecrire la conservervation de theta_l dans le panache descendant
   ! Eventuellement faire la transformation theta_l -> theta_v
   ! Si l'air est sec (et qu'on oublie le côté theta_v) on peut
   ! se contenter de conserver theta.
   !
   ! Connaissant thetadn, on peut calculer la flotabilité.
   ! Connaissant la flotabilité, on peut calculer w de proche en proche
   ! On peut calculer le detrainement de facon à garder alpha*rho = cste
   ! On en déduit l'entrainement latéral
   ! C'est le modèle des mini-projets.

!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
! Initialisations :
!------------------


!
 RETURN
   END