MODULE lmdz_thermcell_env
CONTAINS

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

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


   USE lmdz_thermcell_ini, ONLY: prt_level,RLvCp,RKAPPA,RETV
   USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
   IMPLICIT NONE

! arguments

   INTEGER,INTENT(IN) :: ngrid,nlay,lev_out
   REAL,INTENT(IN), DIMENSION(ngrid,nlay) :: po,pt,pu,pv,pplay
   REAL,INTENT(IN), DIMENSION(ngrid,nlay+1) :: pplev
   REAL,INTENT(OUT), DIMENSION(ngrid,nlay) :: zo,zl,zh,ztv,zthl
   REAL,INTENT(OUT), DIMENSION(ngrid,nlay) :: zpspsk,zu,zv,pqsat
   REAL, INTENT(OUT), DIMENSION(ngrid) :: lintercong
   INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lcong
! Local

   INTEGER ig,ll
   REAL dqsat_dT
   LOGICAL mask(ngrid,nlay)


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

   mask(:,:)=.TRUE.

! 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)
     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

!-----------------------------------------------------------------------
   IF (prt_level>=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
!CR: Calcul du niveau de congelation
   DO ig=1,ngrid
      lcong(ig)=1
      lintercong(ig)=0.
   enddo
   DO ig=1,ngrid
      DO ll=1,nlay-1
         IF ((zh(ig,ll)>273.15).AND.(zh(ig,ll+1)<=273.15)) THEN
            lcong(ig)=ll+1
            lintercong(ig)=(ll*(zh(ig,ll+1)-zh(ig,ll))  &
                 -zh(ig,ll)+273.15)/(zh(ig,ll+1)-zh(ig,ll))
         endif
      enddo
   enddo 
 
 RETURN
   END
END MODULE lmdz_thermcell_env