SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)

!--------------------------------------------------------------
!thermcell_env: calcule les caracteristiques de l environnement
!necessaires au calcul des proprietes dans le thermique
!--------------------------------------------------------------

      IMPLICIT NONE

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

      INTEGER ngrid,nlay
      REAL po(ngrid,nlay)
      REAL pt(ngrid,nlay)
      REAL pu(ngrid,nlay)
      REAL pv(ngrid,nlay)
      REAL pplay(ngrid,nlay)
      REAL pplev(ngrid,nlay+1)
      integer lev_out                           ! niveau pour les print

      REAL zo(ngrid,nlay)
      REAL zl(ngrid,nlay)
      REAL zh(ngrid,nlay)
      REAL ztv(ngrid,nlay)
      REAL zthl(ngrid,nlay)
      REAL zpspsk(ngrid,nlay)
      REAL zu(ngrid,nlay)
      REAL zv(ngrid,nlay)
      REAL zqsat(ngrid,nlay)

      INTEGER ig,l,ll

      real zcor,zdelta,zcvm5,qlbef
      real Tbef,qsatbef
      real dqsat_dT,DT,num,denom
      REAL RLvCp,DDT0
      PARAMETER (DDT0=.01)
      LOGICAL Zsat

      Zsat=.false.
!
! 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
!
! calcul des caracteristiques de l environnement
       DO  ll=1,nlay
         DO ig=1,ngrid
            zo(ig,ll)=po(ig,ll)
            zl(ig,ll)=0.
            zh(ig,ll)=pt(ig,ll)
            zqsat(ig,ll)=0.
         EndDO
       EndDO
!
!
!recherche de saturation dans l environnement
       DO ll=1,nlay
! les points insatures sont definitifs
         DO ig=1,ngrid
            Tbef=pt(ig,ll)
            zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
            qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll)
            qsatbef=MIN(0.5,qsatbef)
            zcor=1./(1.-retv*qsatbef)
            qsatbef=qsatbef*zcor
            Zsat = (max(0.,po(ig,ll)-qsatbef) .gt. 1.e-10)
            if (Zsat) then
            qlbef=max(0.,po(ig,ll)-qsatbef)
! si sature: ql est surestime, d'ou la sous-relax
            DT = 0.5*RLvCp*qlbef
! on pourra enchainer 2 ou 3 calculs sans Do while
            do while (abs(DT).gt.DDT0)
! il faut verifier si c,a conserve quand on repasse en insature ...
              Tbef=Tbef+DT
              zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
              qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll)
              qsatbef=MIN(0.5,qsatbef)
              zcor=1./(1.-retv*qsatbef)
              qsatbef=qsatbef*zcor
! on veut le signe de qlbef
              qlbef=po(ig,ll)-qsatbef
!          dqsat_dT
              zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
              zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
              zcor=1./(1.-retv*qsatbef)
              dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
              num=-Tbef+pt(ig,ll)+RLvCp*qlbef
              denom=1.+RLvCp*dqsat_dT
              if (denom.lt.1.e-10) then
                  print*,'pb denom'
              endif
              DT=num/denom
            enddo
! on ecrit de maniere conservative (sat ou non)
            zl(ig,ll) = max(0.,qlbef)
!          T = Tl +Lv/Cp ql
            zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)
            zo(ig,ll) = po(ig,ll)-zl(ig,ll)
           endif
!on ecrit zqsat 
            zqsat(ig,ll)=qsatbef     
         EndDO
       EndDO
!
!
!-----------------------------------------------------------------------
!   incrementation eventuelle de tendances precedentes:
!   ---------------------------------------------------

      if (prt_level.ge.1) print*,'0 OK convect8'

      DO 1010 l=1,nlay
         DO 1015 ig=1,ngrid
             zpspsk(ig,l)=(pplay(ig,l)/100000.)**RKAPPA
             zu(ig,l)=pu(ig,l)
             zv(ig,l)=pv(ig,l)
!attention zh est maintenant le profil de T et plus le profil de theta !
!
!   T-> Theta
            ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
!Theta_v
            ztv(ig,l)=ztv(ig,l)*(1.+RETV*(zo(ig,l))  &
     &           -zl(ig,l))
!Thetal
            zthl(ig,l)=pt(ig,l)/zpspsk(ig,l)
!            
1015     CONTINUE
1010  CONTINUE
 
      RETURN
      END