MODULE lmdz_thermcell_qsat

  REAL, PARAMETER :: DDT0=.01

CONTAINS

subroutine thermcell_qsat(klon, nlev, active,pplev,ztemp,zqta,zqsat)
USE yoethf_mod_h
  USE yomcst_mod_h


implicit none


  INCLUDE "FCTTRE.h"


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

! Arguments
INTEGER, INTENT(IN) :: klon
INTEGER, INTENT(IN) :: nlev  ! number of vertical to apply qsat
REAL zpspsk(klon, nlev),pplev(klon, nlev)
REAL ztemp(klon, nlev),zqta(klon,nlev),zqsat(klon,nlev)
LOGICAL active(klon, nlev)

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

LOGICAL afaire(klon, nlev),tout_converge
INTEGER :: l
!====================================================================
! INITIALISATIONS
!====================================================================

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


!====================================================================
! 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.
!====================================================================
do l=1,nlev
  do ig=1,klon
     if (active(ig,l)) then
               Tbef(ig,l)=ztemp(ig,l)
               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig,l)))
               qsatbef= R2ES * FOEEW(Tbef(ig,l),zdelta)/pplev(ig,l)
               qsatbef=MIN(0.5,qsatbef)
               zcor=1./(1.-retv*qsatbef)
               qsatbef=qsatbef*zcor
               qlbef=max(0.,zqta(ig,l)-qsatbef)
               DT(ig) = 0.5*RLvCp*qlbef
               zqsat(ig,l)=qsatbef
     endif
  enddo
enddo
! Traitement du cas ou il y a condensation mais faible
! On ne condense pas mais on dit que le qsat est le qta
do l=1,nlev
  do ig=1,klon
     if (active(ig,l)) then
       if (0.<abs(DT(ig,l)).and.abs(DT(ig,l))<=DDT0) then
           zqsat(ig,l)=zqta(ig,l)
        endif
     endif
  enddo
enddo

do iter=1,10
    do l=1,nlev
      afaire(:,l)=abs(DT(:,l)).gt.DDT0
      do ig=1,klon
               if (afaire(ig,l)) then
                 Tbef(ig,l)=Tbef(ig,l)+DT(ig,l)
                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig,l)))
                 qsatbef= R2ES * FOEEW(Tbef(ig,l),zdelta)/pplev(ig,l)
                 qsatbef=MIN(0.5,qsatbef)
                 zcor=1./(1.-retv*qsatbef)
                 qsatbef=qsatbef*zcor
                 qlbef=zqta(ig,l)-qsatbef
                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig,l)))
                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
                 zcor=1./(1.-retv*qsatbef)
                 dqsat_dT=FOEDE(Tbef(ig,l),zdelta,zcvm5,qsatbef,zcor)
                 num=-Tbef(ig,l)+ztemp(ig,l)+RLvCp*qlbef
                 denom=1.+RLvCp*dqsat_dT
                 zqsat(ig,l) = qsatbef
                 DT(ig,l)=num/denom
               endif
      enddo
    enddo
enddo

return
end subroutine thermcell_qsat
END MODULE lmdz_thermcell_qsat