subroutine thermcell_qsat(klon,active,pplev,ztemp,zqta,zqsat)
  USE thermcell_ini_mod, ONLY : RLvCp,RETV, RCPD, RLVTT, RTT
  implicit none

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


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

! Arguments
INTEGER, INTENT(IN) :: klon
!REAL zpspsk(klon),pplev(klon)
!REAL ztemp(klon),zqta(klon),zqsat(klon)
REAL, INTENT(IN) :: pplev(klon)
REAL, INTENT(IN) :: ztemp(klon),zqta(klon)
REAL, INTENT(OUT) :: zqsat(klon)
LOGICAL, INTENT(IN) :: active(klon)

! Variables locales
INTEGER ig,iter
REAL Tbef(klon),DT(klon)
REAL tdelta,qsatbef,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT
logical Zsat
!REAL RLvCp

REAL, SAVE :: DDT0=.01
!$OMP THREADPRIVATE(DDT0)

     !$acc data create (dt, tbef) &
     !$acc present(pplev, ztemp, zqta, active) &
     !$acc present(zqsat) &
     !$acc &

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

!RLvCp = RLVTT/RCPD
!$acc kernels default(none)
DT(:)=0.
!!LF !$acc end kernels

!====================================================================
! Routine a vectoriser en copiant active dans converge et en mettant
! la boucle sur les iterations a l'exterieur est en mettant
! converge= false des que la convergence est atteinte.
!====================================================================

!!LF   !$acc kernels default(none)
do ig=1,klon
   if (active(ig)) then
               Tbef(ig)=ztemp(ig)
               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
               qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
               qsatbef=MIN(0.5,qsatbef)
               zcor=1./(1.-retv*qsatbef)
               qsatbef=qsatbef*zcor
               qlbef=max(0.,zqta(ig)-qsatbef)
               DT(ig) = 0.5*RLvCp*qlbef
               zqsat(ig)=qsatbef
     endif
enddo
!!LF   !$acc end kernels

! Traitement du cas ou il y a condensation mais faible
! On ne condense pas mais on dit que le qsat est le qta
!!LF    !$acc kernels default(none)
do ig=1,klon
   if (active(ig)) then
      if (0.<abs(DT(ig)).and.abs(DT(ig))<=DDT0) then
         zqsat(ig)=zqta(ig)
      endif
   endif
enddo
!!LF   !$acc end kernels

!!LF   !$acc kernels default(none)
!!   !$acc parallel loop private(tbef,dt,zqsat)
do iter=1,10
    do ig=1,klon
       if (abs(DT(ig)).gt.DDT0) then
          Tbef(ig)=Tbef(ig)+DT(ig)
          zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
          qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
          qsatbef=MIN(0.5,qsatbef)
          zcor=1./(1.-retv*qsatbef)
          qsatbef=qsatbef*zcor
          qlbef=zqta(ig)-qsatbef
          zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
          zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
          zcor=1./(1.-retv*qsatbef)
          dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef,zcor)
          num=-Tbef(ig)+ztemp(ig)+RLvCp*qlbef
          denom=1.+RLvCp*dqsat_dT
          zqsat(ig) = qsatbef
          DT(ig)=num/denom
       endif
    enddo
enddo
   !$acc end kernels


   !$acc end data 

return
end