SUBROUTINE borne_2m(klon,   &
         t1, q1, pplay,            &
         ftsol, qsurf, paprs, pctsrf, &
         t2m, q2m,                    &
         zt2m_cor, zq2m_cor,           &
         zrh2m_cor, zqsat2m_cor)
USE yomcst_mod_h
USE yoethf_mod_h
USE indice_sol_mod
IMPLICIT NONE

!==================================================================
! Declarations
!==================================================================

! arguments
INTEGER, INTENT(IN) :: klon
REAL,DIMENSION(klon),INTENT(IN) :: t1, q1
REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: t2m, q2m, pctsrf
REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol
REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: qsurf
REAL,DIMENSION(klon),INTENT(IN) :: pplay
REAL,DIMENSION(klon),INTENT(IN) :: paprs
REAL,DIMENSION (klon),INTENT(OUT) :: zt2m_cor, zq2m_cor
REAL,DIMENSION (klon),INTENT(OUT)  :: zrh2m_cor, zqsat2m_cor

! local
INTEGER                    :: i,nsrf
REAL,DIMENSION(klon)       :: tpot, t2mpot
REAL,DIMENSION(klon,nbsrf) :: t2mpot_cor, t2m_cor, q2m_cor
REAL,DIMENSION(klon,nbsrf) :: rh2m_cor, qsat2m_cor
REAL,DIMENSION(klon,nbsrf) :: tmin, tmax, qmin, qmax, rhmin, rhmax, qsatmin,qsatmax
REAL,DIMENSION(klon)       :: pres2m, zgeo1
REAL                       :: zx_qsa, zx_qss, zx_qs2m, zcor1, zdelta1
REAL                       :: zx_rhs, zx_rh2m, zx_rh1, zref

include "FCTTRE.h"
!==================================================================
! Correction of 2m variables
!==================================================================
zt2m_cor=0.
zq2m_cor=0.
zrh2m_cor=0.
zqsat2m_cor=0.

zref=2
DO i=1,klon
 zgeo1(i) = RD * t1(i) / (0.5*(paprs(i)+pplay(i))) &
               * (paprs(i)-pplay(i))
 pres2m(i) = paprs(i) * (1. - (zref*RG)/zgeo1(i)) + &
             pplay(i) * (zref*RG)/zgeo1(i)
ENDDO 

DO nsrf=1,nbsrf
   DO i=1,klon
!
!!! potential temperature at first level and 2m
!
      tpot(i) = t1(i)* (paprs(i)/pplay(i))**RKAPPA
      t2mpot(i) = t2m(i,nsrf)* (paprs(i)/pres2m(i))**RKAPPA
!
!!! threshold on potential 2m temperature
!
      tmin(i, nsrf) = MIN(tpot(i),ftsol(i,nsrf))
      tmax(i, nsrf) = MAX(tpot(i),ftsol(i,nsrf))
      t2mpot_cor(i,nsrf)=MAX(MIN(t2mpot(i),tmax(i,nsrf)),tmin(i,nsrf))
!
!!! compute t2m temperature from potential 2m temperature
!
      t2m_cor(i,nsrf) = t2mpot_cor(i,nsrf)* (paprs(i)/pres2m(i))**(-RKAPPA)
      zt2m_cor(i)   = zt2m_cor(i)   + t2m_cor(i,nsrf) * pctsrf(i,nsrf)
!
!!! q2m should be between surface and 1st level values
!
      qmin(i, nsrf) = MIN(q1(i),qsurf(i,nsrf))
      qmax(i, nsrf) = MAX(q1(i),qsurf(i,nsrf))
      q2m_cor(i,nsrf)=MAX(MIN(q2m(i,nsrf),qmax(i,nsrf)),qmin(i,nsrf))
      zq2m_cor(i)   = zq2m_cor(i)   + q2m_cor(i,nsrf) * pctsrf(i,nsrf)
!
!!! saturation humidity at first level
!
      zdelta1 = MAX(0.,SIGN(1., rtt-t1(i)))
      zx_qsa  = r2es * FOEEW(t1(i),zdelta1)/pplay(i)
      zx_qsa  = MIN(0.5,zx_qsa)
      zcor1   = 1./(1.-RETV*zx_qsa)
      zx_qsa  = zx_qsa*zcor1
      zx_rh1  = q1(i)/zx_qsa
!
!!! saturation humidity at 2m
!
      zdelta1 = MAX(0.,SIGN(1., rtt-t2m_cor(i,nsrf) ))
      zx_qs2m  = r2es * FOEEW(t2m_cor(i,nsrf),zdelta1)/pres2m(i)
      zx_qs2m  = MIN(0.5,zx_qs2m)
      zcor1   = 1./(1.-RETV*zx_qs2m)
      zx_qs2m  = zx_qs2m*zcor1
      zx_rh2m  = q2m_cor(i,nsrf)/zx_qs2m
!
!!! saturation humidity at the surface
!
      if( nsrf == is_ter ) THEN
       zdelta1 = MAX(0.,SIGN(1., rtt-ftsol(i,nsrf)))
       zx_qss  = r2es * FOEEW(ftsol(i,nsrf),zdelta1)/paprs(i)
       zx_qss  = MIN(0.5,zx_qss)
       zcor1   = 1./(1.-RETV*zx_qss)
       zx_qss  = zx_qss*zcor1
       zx_rhs  = qsurf(i,nsrf)/zx_qss
!
!!!IM cf. EV : do not allow over saturation at the surface over the land
!
       zx_rhs  = MIN(qsurf(i,nsrf)/zx_qss,1.)
      else
       zx_rhs=1.
      endif
!
!!! rh2m should be between surface and 1st level values
!
      rhmin(i,nsrf)=MIN(zx_rhs,zx_rh1)
      rhmax(i,nsrf)=MAX(zx_rhs,zx_rh1)
      rh2m_cor(i,nsrf)=MAX(MIN(zx_rh2m,rhmax(i,nsrf)),rhmin(i,nsrf))
      zrh2m_cor(i)   = zrh2m_cor(i)   + rh2m_cor(i,nsrf) * pctsrf(i,nsrf)
!
!!! qsat2m must be comprised between surface and 1st level values
!
      qsatmin(i, nsrf)=MIN(zx_qss,zx_qsa)
      qsatmax(i, nsrf)=MAX(zx_qss,zx_qsa)
      qsat2m_cor(i,nsrf)=MAX(MIN(zx_qs2m,qsatmax(i,nsrf)),qsatmin(i,nsrf))
      zqsat2m_cor(i) = zqsat2m_cor(i) + qsat2m_cor(i,nsrf) * pctsrf(i,nsrf)
!
   ENDDO
ENDDO
RETURN
END SUBROUTINE borne_2m