MODULE lmdz_thermcell_qsat
CONTAINS

  SUBROUTINE thermcell_qsat(klon, active, pplev, ztemp, zqta, zqsat)
    USE lmdz_yoethf
    USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
    USE lmdz_yomcst

    IMPLICIT NONE

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

    ! Arguments
    INTEGER klon
    REAL zpspsk(klon), pplev(klon)
    REAL ztemp(klon), zqta(klon), zqsat(klon)
    LOGICAL 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)

    LOGICAL afaire(klon), tout_converge

    !====================================================================
    ! 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 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
    END DO

    ! Traitement du cas ou il y a condensation mais faible
    ! On ne condense pas mais on dit que le qsat est le qta
    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
    END DO

    DO iter = 1, 10
      afaire(:) = abs(DT(:))>DDT0
      do ig = 1, klon
        IF (afaire(ig)) 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
    END DO

    RETURN
  END
END MODULE lmdz_thermcell_qsat