SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out &
#ifdef ISO
     &            ,xtpo,xtzo,xtzl &
#endif      
     &   )

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

      USE print_control_mod, ONLY: prt_level
#ifdef ISO
  USE infotrac_phy, ONLY : ntraciso
  USE isotopes_mod, ONLY: pxtice,pxtmelt,iso_eau,iso_HDO
  USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall
#ifdef ISOVERIF
  USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
        iso_verif_aberrant_encadre,iso_verif_noNaN
#endif
#ifdef ISOTRAC
  USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall_trac
#ifdef ISOVERIF
  USE isotopes_verif_mod, ONLY: iso_verif_traceur
#endif
#endif
#endif
      IMPLICIT NONE

#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.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 pqsat(ngrid,nlay)

      INTEGER ig,ll

      real dqsat_dT
      real RLvCp

logical mask(ngrid,nlay)

#ifdef ISO
      REAL xtpo(ntraciso,ngrid,nlay)
      REAL xtzo(ntraciso,ngrid,nlay)
      REAL xtzl(ntraciso,ngrid,nlay)
      integer ixt
      real zxtliq(ntraciso,ngrid),zxtice(ntraciso,ngrid)
      real zfice(ngrid)
      REAL zxtl(ntraciso,ngrid,nlay)
#endif 


!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
! Initialisations :
!------------------

mask(:,:)=.true.
RLvCp = RLVTT/RCPD

!
! 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)
#ifdef ISO
            do ixt=1,ntraciso
              xtzo(ixt,ig,ll)=xtpo(ixt,ig,ll)
              xtzl(ixt,ig,ll)=0. 
            enddo   !do ixt=1,ntraciso 
#endif
         EndDO
       EndDO
!
!
! Condensation :
!---------------
! Calcul de l'humidite a saturation et de la condensation

call thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat)
DO ll=1,nlay
   DO ig=1,ngrid
      zl(ig,ll) = max(0.,po(ig,ll)-pqsat(ig,ll))
      zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)         !   T = Tl + Lv/Cp ql
      zo(ig,ll) = po(ig,ll)-zl(ig,ll)
   ENDDO
ENDDO

#ifdef ISO
#ifdef ISOVERIF
   DO ll=1,nlay
    do ig=1,ngrid 
       call iso_verif_noNaN(po(ig,ll), &
       & 'thermcell_env 99')
       do ixt=1,ntraciso         
         call iso_verif_noNaN(xtpo(ixt,ig,ll), &
       &        'thermcell_env 103')
       enddo !do ixt=1,ntraciso
    enddo !do ig=1,ngrid 
   enddo !DO ll=1,nlay
#endif
#ifdef ISOVERIF 
   write(*,*) 'thermcell_env 94: avant condiso_liq_ice_vectall'
   DO ll=1,nlay
    do ig=1,ngrid   
     if (iso_eau.gt.0) then
       call iso_verif_egalite(xtpo(iso_eau,ig,ll),po(ig,ll), &
       & 'thermcell_env 100')
     endif
    enddo !do ig=1,ngrid 
   enddo !DO ll=1,nlay
#endif      
 DO ll=1,nlay
  do ig=1,ngrid
        zfice(ig)=1.0-(pt(ig,ll)-pxtice)/(pxtmelt-pxtice)
        zfice(ig) = MIN(MAX(zfice(ig),0.0),1.0)
  enddo  !do ig=1,ngrid 

  call condiso_liq_ice_vectall(xtpo(1,1,ll),po(1,ll),zl(1,ll),pt(1,ll),zfice, &
     &            zxtice,zxtliq,ngrid)
#ifdef ISOTRAC
  call condiso_liq_ice_vectall_trac(xtpo(1,1,ll),po(1,ll),zl(1,ll), &
     &            pt(1,ll),zfice,zxtice,zxtliq,ngrid)
#ifdef ISOVERIF
  do ig=1,ngrid            
      call iso_verif_traceur(xtpo(1,ig,ll),'thermcell_env 102')
      call iso_verif_traceur(zxtliq(1,ig),'thermcell_env 103')
      call iso_verif_traceur(zxtice(1,ig),'thermcell_env 104')
  enddo !do ig=1,ngrid 
#endif
#endif
  do ig=1,ngrid
     do ixt=1,ntraciso                 
       zxtl(ixt,ig,ll)=zxtice(ixt,ig)+zxtliq(ixt,ig)
       xtzo(ixt,ig,ll) = xtpo(ixt,ig,ll)-zxtl(ixt,ig,ll)
     enddo ! do ixt=1,ntraciso  
#ifdef ISOVERIF
     if (iso_eau.gt.0) then
       call iso_verif_egalite(zxtl(iso_eau,ig,ll),zl(ig,ll), &
       & 'thermcell_env 120')
       call iso_verif_egalite(xtzo(iso_eau,ig,ll),zo(ig,ll), &
       & 'thermcell_env 136')
     endif
      if (iso_HDO.gt.0) then
      call iso_verif_aberrant_encadre( &
       &          xtzo(iso_HDO,ig,ll)/zo(ig,ll), &
       &          'thermcell_env 141')
      endif     
#endif     
  enddo ! do ig=1,ngrid
enddo !DO ll=1,nlay
#endif
!
!
!-----------------------------------------------------------------------

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

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