subroutine thermcell_qsat(klon,mask,pplev,ptemp,pqt,pqsat)
implicit none

#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"

!====================================================================
! DECLARATIONS
!====================================================================

! Arguments
INTEGER klon
REAL pplev(klon)
REAL ptemp(klon),pqt(klon),pqsat(klon)
LOGICAL mask(klon)

! Variables locales
INTEGER ig,iter
REAL Tbef(klon),DT(klon)
REAL tdelta,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT
logical Zsat
REAL RLvCp
REAL, SAVE :: DDT0=.01
LOGICAL afaire(klon),tout_converge

!====================================================================
! INITIALISATIONS
!====================================================================

RLvCp = RLVTT/RCPD
tout_converge=.false.
afaire(:)=mask(:)
DT(:)=0.


!====================================================================
! Routine a vectoriser en copiant mask dans converge et en mettant
! la boucle sur les iterations a l'exterieur est en mettant
! converge= false des que la convergence est atteinte.
! Pr Tprec=Tl calcul de qsat
! Si qsat>qT T=Tl, q=qT
! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
! On cherche DDT < DDT0
!====================================================================

do ig=1,klon
   if (mask(ig)) then
               Tbef(ig)=ptemp(ig)
               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
               pqsat(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
               pqsat(ig)=MIN(0.5,pqsat(ig))
               zcor=1./(1.-retv*pqsat(ig))
               pqsat(ig)=pqsat(ig)*zcor
               qlbef=max(0.,pqt(ig)-pqsat(ig))
               afaire(ig)=qlbef>1.e-10
               DT(ig) = 0.5*RLvCp*qlbef
     endif
enddo


do iter=1,10
    afaire(:)=(abs(DT(:))>DDT0).and.afaire(:)
    do ig=1,klon
               if (afaire(ig)) then
                 Tbef(ig)=Tbef(ig)+DT(ig)
                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
                 pqsat(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
                 pqsat(ig)=MIN(0.5,pqsat(ig))
                 zcor=1./(1.-retv*pqsat(ig))
                 pqsat(ig)=pqsat(ig)*zcor
                 qlbef=pqt(ig)-pqsat(ig)
                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
                 zcor=1./(1.-retv*pqsat(ig))
                 dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,pqsat(ig),zcor)
                 num=-Tbef(ig)+ptemp(ig)+RLvCp*qlbef
                 denom=1.+RLvCp*dqsat_dT
                 DT(ig)=num/denom
               endif
    enddo
enddo

return
end