subroutine SISVAT_TS2
  ! #ES.                     (ETSo_0,ETSo_1,ETSo_d)

  ! +------------------------------------------------------------------------+
  ! | MAR          SISVAT_TS2                            Mon 16-08-2009  MAR |
  ! |   SubRoutine SISVAT_TS2 computes the Soil/Snow temperature and fluxes  |
  ! |   using the same method as in LMDZ for consistency.                    |
  ! |   The corresponding LMDZ routines are soil (soil.F90) and calcul_fluxs |
  ! |   (calcul_fluxs_mod.F90).                                              |
  ! +------------------------------------------------------------------------+
  ! |                                                                        |
  ! |                                                                        |
  ! |   PARAMETERS:  klonv: Total Number of columns =                        |
  ! |   ^^^^^^^^^^        = Total Number of grid boxes of surface type       |
  ! |                       (land ice for now)                               |
  ! |                                                                        |
  ! |   INPUT:   isnoSV   = total Nb of Ice/Snow Layers                      |
  ! |   ^^^^^    sol_SV   : Downward Solar Radiation                  [W/m2] |
  ! |            IRd_SV   : Surface Downward  Longwave   Radiation    [W/m2] |
  ! |            VV__SV   : SBL Top    Wind Speed                      [m/s] |
  ! |            TaT_SV   : SBL Top    Temperature                       [K] |
  ! |            QaT_SV   : SBL Top    Specific  Humidity            [kg/kg] |
  ! |            dzsnSV   : Snow Layer Thickness                         [m] |
  ! |            dt__SV   : Time Step                                    [s] |
  ! |                                                                        |
  ! |            SoSosv   : Absorbed Solar Radiation by Surfac.(Normaliz)[-] |
  ! |            Eso_sv   : Soil+Snow       Emissivity                   [-] |
  ! |   ?        rah_sv   : Aerodynamic Resistance for Heat            [s/m] |
  ! |                                                                        |
  ! |            dz1_SV    : "inverse" layer thickness (centered)      [1/m] |
  ! |            dz2_SV    : layer thickness (layer above (?))           [m] |
  ! |            AcoHSV    : coefficient for Enthalpy evolution, from atm.   |
  ! |            AcoHSV    : coefficient for Enthalpy evolution, from atm.   |
  ! |            AcoQSV    : coefficient for Humidity evolution, from atm.   |
  ! |            BcoQSV    : coefficient for Humidity evolution, from atm.   |
  ! |            ps__SV    : surface pressure                           [Pa] |
  ! |            p1l_SV    : 1st atmospheric layer pressure             [Pa] |
  ! |            cdH_SV    : drag coeff Energy (?)                           |
  ! |            rsolSV    : Radiation balance surface                [W/m2] |
  ! |            lambSV    : Coefficient for soil layer geometry         [-] |
  ! |                                                                        |
  ! |   INPUT /  TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
  ! |   OUTPUT:           & Snow     Temperatures (layers  1,2,...,nsno) [K] |
  ! |   ^^^^^^   rsolSV   : Radiation balance surface                 [W/m2] |
  ! |                                                                        |
  ! |   OUTPUT:  IRs_SV   : Soil      IR Radiation                    [W/m2] |
  ! |   ^^^^^^   HSs_sv   : Sensible  Heat Flux                       [W/m2] |
  ! |            HLs_sv   : Latent    Heat Flux                       [W/m2] |
  ! |            TsfnSV   : new surface temperature                      [K] |
  ! |            Evp_sv   : Evaporation                              [kg/m2] |
  ! |            dSdTSV   : Sensible Heat Flux temp. derivative     [W/m2/K] |
  ! |            dLdTSV   : Latent Heat Flux temp. derivative       [W/m2/K] |
  ! |                                                                        |
  ! |   ?        ETSo_0   : Snow/Soil Energy Power, before Forcing    [W/m2] |
  ! |   ?        ETSo_1   : Snow/Soil Energy Power, after  Forcing    [W/m2] |
  ! |   ?        ETSo_d   : Snow/Soil Energy Power         Forcing    [W/m2] |
  ! |                                                                        |
  ! |________________________________________________________________________|

  USE yoethf_mod_h
  USE VAR_SV
  USE VARdSV
  USE VARySV
  USE VARtSV
  USE VARxSV
  USE VARphy
  USE indice_sol_mod
  USE yomcst_mod_h
  USE comsoil_mod_h
IMPLICIT NONE


  ! +--Global Variables
  ! +  ================
  INCLUDE "FCTTRE.h"

  ! +--OUTPUT for Stand Alone NetCDF File
  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! #NC      real*8        SOsoKL(klonv)             ! Absorbed Solar Radiation
  ! #NC      real*8        IRsoKL(klonv)             ! Absorbed IR    Radiation
  ! #NC      real*8        HSsoKL(klonv)             ! Absorbed Sensible Heat Flux
  ! #NC      real*8        HLsoKL(klonv)             ! Absorbed Latent   Heat Flux
  ! #NC      real*8        HLs_KL(klonv)             ! Evaporation
  ! #NC      real*8        HLv_KL(klonv)             ! Transpiration
  ! #NC      common/DumpNC/SOsoKL,IRsoKL
  ! #NC     .             ,HSsoKL,HLsoKL
  ! #NC     .             ,HLs_KL,HLv_KL

  ! +--Internal Variables
  ! +  ==================

  integer :: ig,jk,isl
  real :: mu
  real :: Tsrf(klonv)               ! surface temperature as extrapolated from soil
  real :: mug(klonv)                 !hj coef top layers
  real :: ztherm_i(klonv),zdz2(klonv,-nsol:nsno),z1s
  real :: pfluxgrd(klonv), pcapcal(klonv), cal(klonv)
  real :: beta(klonv), dif_grnd(klonv)
  real :: C_coef(klonv,-nsol:nsno),D_coef(klonv,-nsol:nsno)

  REAL, DIMENSION(klonv)   :: zx_mh, zx_nh, zx_oh
  REAL, DIMENSION(klonv)   :: zx_mq, zx_nq, zx_oq
  REAL, DIMENSION(klonv)   :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
  REAL, DIMENSION(klonv)   :: zx_sl, zx_k1
  REAL, DIMENSION(klonv)   :: d_ts
  REAL                     :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
  REAL                     :: qsat_new, q1_new
   ! REAL, PARAMETER          :: t_grnd = 271.35, t_coup = 273.15
   ! REAL, PARAMETER          :: max_eau_sol = 150.0
  REAL, DIMENSION(klonv)   :: IRs__D, dIRsdT


  REAL :: t_grnd                      ! not used
  parameter(t_grnd = 271.35)       !
  REAL :: t_coup                      ! distinguish evap/sublimation
  parameter(t_coup = 273.15)       !
  REAL :: max_eau_sol
  parameter(max_eau_sol = 150.0)


     ! write(*,*)'T check'
  !
  !    DO  ig = 1,knonv
  !        DO  jk = 1,isnoSV(ig) !nsno
  !          IF (TsisSV(ig,jk) <= 1.) THEN          !hj check
  !            TsisSV(ig,jk) = TsisSV(ig,isnoSV(ig))
  !          ENDIF
  !
  !          IF (TsisSV(ig,jk) <= 1.) THEN          !hj check
  !            TsisSV(ig,jk) = 273.15
  !          ENDIF
  !        END DO
  !        END DO

  !!=======================================================================
  !! I. First part: corresponds to soil.F90 in LMDZ
  !!=======================================================================

  DO ig = 1,knonv
    DO jk =1,isnoSV(ig)
      dz2_SV(ig,jk)=dzsnSV(ig,jk)
  !! use arithmetic center between layers to derive dz1 for snow layers for simplicity:
      dz1_SV(ig,jk)=2./(dzsnSV(ig,jk)+dzsnSV(ig,jk-1))
    ENDDO
  ENDDO

  DO ig = 1,knonv
    ztherm_i(ig)   = inertie_lic
    IF (isnoSV(ig) > 0) ztherm_i(ig)   = inertie_sno
  ENDDO

  !!-----------------------------------------------------------------------
  !! 1)
  !! Calculation of Cgrf and Dgrd coefficients using soil temperature from
  !! previous time step.
  !!
  !! These variables are recalculated on the local compressed grid instead
  !! of saved in restart file.
  !!-----------------------------------------------------------------------
  DO ig=1,knonv
    DO jk=-nsol,nsno
      zdz2(ig,jk)=dz2_SV(ig,jk)/dt__SV                                !ptimestep
    ENDDO
  ENDDO

  DO ig=1,knonv
    z1s = zdz2(ig,-nsol)+dz1_SV(ig,-nsol+1)
    C_coef(ig,-nsol+1)=zdz2(ig,-nsol)*TsisSV(ig,-nsol)/z1s
    D_coef(ig,-nsol+1)=dz1_SV(ig,-nsol+1)/z1s
  ENDDO

  DO ig=1,knonv
    DO jk=-nsol+1,isnoSV(ig)-1,1
      z1s = 1./(zdz2(ig,jk)+dz1_SV(ig,jk+1)+dz1_SV(ig,jk)           &
            *(1.-D_coef(ig,jk)))
      C_coef(ig,jk+1)=                                              &
            (TsisSV(ig,jk)*zdz2(ig,jk)                              &
            +dz1_SV(ig,jk)*C_coef(ig,jk)) * z1s
      D_coef(ig,jk+1)=dz1_SV(ig,jk+1)*z1s
    ENDDO
  ENDDO

  !!-----------------------------------------------------------------------
  !! 2)
  !! Computation of the soil temperatures using the Cgrd and Dgrd
  !! coefficient computed above
  !!
  !!-----------------------------------------------------------------------
  !! Extrapolate surface Temperature                   !hj check
  mu=1./((2.**1.5-1.)/(2.**(0.5)-1.)-1.)

  ! IF (knonv>0) THEN
  !  DO ig=1,8
  !    write(*,*)ig,'sisvat: Tsis ',TsisSV(ig,isnoSV(ig))
  !    write(*,*)'max-1            ',TsisSV(ig,isnoSV(ig)-1)
  !    write(*,*)'max-2            ',TsisSV(ig,isnoSV(ig)-2)
  !    write(*,*)'0                ',TsisSV(ig,0)
  !!        write(*,*)min(max(isnoSV(ig),0),1),max(1-isnoSV(ig),0)
  !  ENDDO
  ! END IF

  DO ig=1,knonv
    IF (isnoSV(ig).GT.0) THEN
      IF (isnoSV(ig).GT.1) THEN
       mug(ig)=1./(1.+dzsnSV(ig,isnoSV(ig)-1)/dzsnSV(ig,isnoSV(ig))) !mu
      ELSE
       mug(ig) = 1./(1.+dzsnSV(ig,isnoSV(ig)-1)/dz_dSV(0)) !mu
      ENDIF
    ELSE
      mug(ig) = lambSV
    ENDIF

    IF (mug(ig)  .LE. 0.05) THEN
     write(*,*)'Attention mu low', mug(ig)
    ENDIF
    IF (mug(ig)  .GE. 0.98) THEN
     write(*,*)'Attention mu high', mug(ig)
    ENDIF

    Tsrf(ig)=(1.5*TsisSV(ig,isnoSV(ig))-0.5*TsisSV(ig,isnoSV(ig)-1))&
          *min(max(isnoSV(ig),0),1)+                             &
          ((mug(ig)+1)*TsisSV(ig,0)-mug(ig)*TsisSV(ig,-1))       &
          *max(1-isnoSV(ig),0)
  ENDDO



  !! Surface temperature
  DO ig=1,knonv
  TsisSV(ig,isnoSV(ig))=(mug(ig)*C_coef(ig,isnoSV(ig))+Tsf_SV(ig))/ &
        (mug(ig)*(1.-D_coef(ig,isnoSV(ig)))+1.)
  ENDDO

  !! Other temperatures
  DO ig=1,knonv
    DO jk=isnoSV(ig),-nsol+1,-1
      TsisSV(ig,jk-1)=C_coef(ig,jk)+D_coef(ig,jk)                   &
            *TsisSV(ig,jk)
    ENDDO
  ENDDO
   ! write(*,*)ig,'Tsis',TsisSV(ig,0)

   ! IF (indice == is_sic) THEN
   !   DO ig = 1,knonv
   !     TsisSV(ig,-nsol) = RTT - 1.8
   !   END DO
   ! ENDIF

  !C      !hj new 11 03 2010
    DO ig=1,knonv
      isl         = isnoSV(ig)
       ! dIRsdT(ig) = Eso_sv(ig)* SteBo   * 4.                        & ! - d(IR)/d(T)
  ! &                             * Tsf_SV(ig)                         & !T TsisSV(ig,isl)           !
  ! &                             * Tsf_SV(ig)                         & !TsisSV(ig,isl)           !
  ! &                             * Tsf_SV(ig) !TsisSV(ig,isl)           !
  !      IRs__D(ig) = dIRsdT(ig)* Tsf_SV(ig) * 0.75 !TsisSV(ig,isl) * 0.75    !:
      dIRsdT(ig) = Eso_sv(ig)* StefBo   * 4.                        & ! - d(IR)/d(T)
            * TsisSV(ig,isl)                     &
            * TsisSV(ig,isl)                     &
            * TsisSV(ig,isl)
      IRs__D(ig) = dIRsdT(ig)* TsisSV(ig,isl) * 0.75                  !:
     END DO
   ! !hj
  !!-----------------------------------------------------------------------
  !! 3)
  !! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil
  !! temperature
  !!-----------------------------------------------------------------------
  DO ig=1,knonv
    z1s = zdz2(ig,-nsol)+dz1_SV(ig,-nsol+1)
    C_coef(ig,-nsol+1) = zdz2(ig,-nsol)*TsisSV(ig,-nsol)/z1s
    D_coef(ig,-nsol+1) = dz1_SV(ig,-nsol+1)/z1s
  ENDDO

  DO ig=1,knonv
    DO jk=-nsol+1,isnoSV(ig)-1,1
      z1s = 1./(zdz2(ig,jk)+dz1_SV(ig,jk+1)+dz1_SV(ig,jk)           &
            *(1.-D_coef(ig,jk)))
      C_coef(ig,jk+1) = (TsisSV(ig,jk)*zdz2(ig,jk)+                 &
            dz1_SV(ig,jk)*C_coef(ig,jk)) * z1s
      D_coef(ig,jk+1) = dz1_SV(ig,jk+1)*z1s
    ENDDO
  ENDDO

  !!-----------------------------------------------------------------------
  !! 4)
  !! Computation of the surface diffusive flux from ground and
  !! calorific capacity of the ground
  !!-----------------------------------------------------------------------
  DO ig=1,knonv
  !! (pfluxgrd)
    pfluxgrd(ig) = ztherm_i(ig)*dz1_SV(ig,isnoSV(ig))*              &
          (C_coef(ig,isnoSV(ig))+(D_coef(ig,isnoSV(ig))-1.)           &
          *TsisSV(ig,isnoSV(ig)))
  !! (pcapcal)
    pcapcal(ig)  = ztherm_i(ig)*                                    &
          (dz2_SV(ig,isnoSV(ig))+dt__SV*(1.-D_coef(ig,isnoSV(ig)))    &
          *dz1_SV(ig,isnoSV(ig)))
    z1s = mug(ig)*(1.-D_coef(ig,isnoSV(ig)))+1.
    pcapcal(ig)  = pcapcal(ig)/z1s
    pfluxgrd(ig) = ( pfluxgrd(ig)                                   &
          + pcapcal(ig) * (TsisSV(ig,isnoSV(ig)) * z1s               &
          - mug(ig)* C_coef(ig,isnoSV(ig))                           &
          - Tsf_SV(ig))       /dt__SV )
  ENDDO


  cal(1:knonv) = RCPD / pcapcal(1:knonv)
  rsolSV(1:knonv)  = rsolSV(1:knonv) + pfluxgrd(1:knonv)
  !!=======================================================================
  !! II. Second part: corresponds to calcul_fluxs_mod.F90 in LMDZ
  !!=======================================================================

  Evp_sv = 0.
  ! #NC HSsoKL=0.
  ! #NC HLsoKL=0.
  dSdTSV = 0.
  dLdTSV = 0.

  beta(:) = 1.0
  dif_grnd(:) = 0.0

  !! zx_qs = qsat en kg/kg
  !!**********************************************************************x***************

  DO ig = 1,knonv
    IF (ps__SV(ig).LT.1.) THEN
       ! write(*,*)'ig',ig,'ps',ps__SV(ig)
      ps__SV(ig)=max(ps__SV(ig),1.e-8)
    ENDIF
    IF (p1l_SV(ig).LT.1.) THEN
       ! write(*,*)'ig',ig,'p1l',p1l_SV(ig)
      p1l_SV(ig)=max(p1l_SV(ig),1.e-8)
    ENDIF
    IF (TaT_SV(ig).LT.180.) THEN
       ! write(*,*)'ig',ig,'TaT',TaT_SV(ig)
      TaT_SV(ig)=max(TaT_SV(ig),180.)
    ENDIF
    IF (QaT_SV(ig).LT.1.e-8) THEN
       ! write(*,*)'ig',ig,'QaT',QaT_SV(ig)
      QaT_SV(ig)=max(QaT_SV(ig),1.e-8)
    ENDIF
    IF (Tsf_SV(ig).LT.100.) THEN
       ! write(*,*)'ig',ig,'Tsf',Tsf_SV(ig)
      Tsf_SV(ig)=max(Tsf_SV(ig),180.)
    ENDIF
    IF (Tsf_SV(ig).GT.500.) THEN
       ! write(*,*)'ig',ig,'Tsf',Tsf_SV(ig)
      Tsf_SV(ig)=min(Tsf_SV(ig),400.)
    ENDIF
      ! IF (Tsrf(ig).LT.1.) THEN
  !!          write(*,*)'ig',ig,'Tsrf',Tsrf(ig)
      !   Tsrf(ig)=max(Tsrf(ig),TaT_SV(ig)-20.)
      ! ENDIF
     IF (cdH_SV(ig).LT.1.e-10) THEN
        ! IF (ig.le.3)   write(*,*)'ig',ig,'cdH',cdH_SV(ig)
       cdH_SV(ig)=.5
     ENDIF
  ENDDO


  DO ig = 1,knonv
    zx_pkh(ig) = 1. ! (ps__SV(ig)/ps__SV(ig))**RKAPPA
    IF (thermcep) THEN
      zdelta=MAX(0.,SIGN(1.,rtt-Tsf_SV(ig)))
      zcvm5 = R5LES*LhvH2O*(1.-zdelta) + R5IES*LhsH2O*zdelta
      zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*QaT_SV(ig))
      zx_qs= r2es * FOEEW(Tsf_SV(ig),zdelta)/ps__SV(ig)
      zx_qs=MIN(0.5,zx_qs)
      ! !write(*,*)'zcor',retv*zx_qs
      zcor=1./(1.-retv*zx_qs)
      zx_qs=zx_qs*zcor
      zx_dq_s_dh = FOEDE(Tsf_SV(ig),zdelta,zcvm5,zx_qs,zcor)        &
            /LhvH2O / zx_pkh(ig)
    ELSE
      IF (Tsf_SV(ig).LT.t_coup) THEN
         zx_qs = qsats(Tsf_SV(ig)) / ps__SV(ig)
         zx_dq_s_dh = dqsats(Tsf_SV(ig),zx_qs)/LhvH2O               &
               / zx_pkh(ig)
      ELSE
         zx_qs = qsatl(Tsf_SV(ig)) / ps__SV(ig)
         zx_dq_s_dh = dqsatl(Tsf_SV(ig),zx_qs)/LhvH2O               &
               / zx_pkh(ig)
      ENDIF
    ENDIF
    zx_dq_s_dt(ig) = RCPD * zx_pkh(ig) * zx_dq_s_dh
    zx_qsat(ig) = zx_qs
     ! zx_coef(ig) = cdH_SV(ig) *                                     &
  ! &       (1.0+SQRT(u1lay(ig)**2+v1lay(ig)**2)) *                   &
  ! &       p1l_SV(ig)/(RD*t1lay(ig))
    zx_coef(ig) = cdH_SV(ig) *                                      &
          (1.0+VV__SV(ig)) *                                         &
          p1l_SV(ig)/(RD*TaT_SV(ig))

  ENDDO


  !! === Calcul de la temperature de surface ===
  !! zx_sl = chaleur latente d'evaporation ou de sublimation
  !!****************************************************************************

  DO ig = 1,knonv
    zx_sl(ig) = LhvH2O
    IF (Tsf_SV(ig) .LT. RTT) zx_sl(ig) = LhsH2O
    zx_k1(ig) = zx_coef(ig)
  ENDDO


  DO ig = 1,knonv
  !! Q
    zx_oq(ig) = 1. - (beta(ig) * zx_k1(ig) * BcoQSV(ig) * dt__SV)
    zx_mq(ig) = beta(ig) * zx_k1(ig) *                              &
          (AcoQSV(ig) - zx_qsat(ig) +                                &
          zx_dq_s_dt(ig) * Tsf_SV(ig))                               &
          / zx_oq(ig)
    zx_nq(ig) = beta(ig) * zx_k1(ig) * (-1. * zx_dq_s_dt(ig))       &
          / zx_oq(ig)

  !! H
    zx_oh(ig) = 1. - (zx_k1(ig) * BcoHSV(ig) * dt__SV)
    zx_mh(ig) = zx_k1(ig) * AcoHSV(ig) / zx_oh(ig)
    zx_nh(ig) = - (zx_k1(ig) * RCPD * zx_pkh(ig))/ zx_oh(ig)

  !! surface temperature
    TsfnSV(ig) = (Tsf_SV(ig) + cal(ig)/RCPD * zx_pkh(ig) * dt__SV * &
          (rsolSV(ig) + zx_mh(ig) + zx_sl(ig) * zx_mq(ig))           &
          + dif_grnd(ig) * t_grnd * dt__SV)/                         &
          ( 1. - dt__SV * cal(ig)/(RCPD * zx_pkh(ig)) *              &
          (zx_nh(ig) + zx_sl(ig) * zx_nq(ig))                        &
          + dt__SV * dif_grnd(ig))

  !hj rajoute 22 11 2010 tuning...
    TsfnSV(ig) = min(RTT+0.02,TsfnSV(ig))

    d_ts(ig) = TsfnSV(ig) - Tsf_SV(ig)


  !!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
  !!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
    Evp_sv(ig) = - zx_mq(ig) - zx_nq(ig) * TsfnSV(ig)
    HLs_sv(ig) = - Evp_sv(ig) * zx_sl(ig)
    HSs_sv(ig) = zx_mh(ig) + zx_nh(ig) * TsfnSV(ig)

  !! Derives des flux dF/dTs (W m-2 K-1):
    dSdTSV(ig) = zx_nh(ig)
    dLdTSV(ig) = zx_sl(ig) * zx_nq(ig)


  !hj  new 11 03 2010
    isl         = isnoSV(ig)
     ! TsisSV(ig,isl) = TsfnSV(ig)
    IRs_SV(ig) = IRs__D(ig)                                         & !
          - dIRsdT(ig) * TsfnSV(ig) !TsisSV(ig,isl)?      !

  ! hj
  ! #NC   SOsoKL(ig) = sol_SV(ig) * SoSosv(ig)              ! Absorbed Sol.
  ! #NC   IRsoKL(ig) =               IRs_SV(ig)                           & !Up Surf. IR
  ! #NC&        +     tau_sv(ig)      *IRd_SV(ig)*Eso_sv(ig)              & !Down Atm IR
  ! #NC&        -(1.0-tau_sv(ig)) *0.5*IRv_sv(ig)            ! Down Veg IR
  ! #NC   HLsoKL(ig) = HLs_sv(ig)
  ! #NC   HSsoKL(ig) = HSs_sv(ig)
  ! #NC   HLs_KL(ig) = Evp_sv(ig)

  !! Nouvelle valeure de l'humidite au dessus du sol
    qsat_new=zx_qsat(ig) + zx_dq_s_dt(ig) * d_ts(ig)
    q1_new = AcoQSV(ig) - BcoQSV(ig)* Evp_sv(ig)*dt__SV
    QaT_SV(ig)=q1_new*(1.-beta(ig)) + beta(ig)*qsat_new

  ENDDO

end subroutine sisvat_ts2
