



subroutine SISVAT_TSo
  ! #e1.                     (ETSo_0,ETSo_1,ETSo_d)

  ! +------------------------------------------------------------------------+
  ! | MAR          SISVAT_TSo                                06-10-2020  MAR |
  ! |   SubRoutine SISVAT_TSo computes the Soil/Snow Energy Balance          |
  ! +------------------------------------------------------------------------+
  ! |                                                                        |
  ! |   PARAMETERS:  knonv: Total Number of columns =                        |
  ! |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
  ! |                     X       Number of Mosaic Cell per grid box         |
  ! |                                                                        |
  ! |   INPUT:   isotSV   = 0,...,11:   Soil       Type                      |
  ! |   ^^^^^               0:          Water, Solid or Liquid               |
  ! |            isnoSV   = total Nb of Ice/Snow Layers                      |
  ! |            dQa_SV   = Limitation of  Water Vapor  Turbulent Flux       |
  ! |                                                                        |
  ! |   INPUT:   sol_SV   : Downward Solar Radiation                  [W/m2] |
  ! |   ^^^^^    IRd_SV   : Surface Downward  Longwave   Radiation    [W/m2] |
  ! |            za__SV   : SBL Top    Height                            [m] |
  ! |            VV__SV   : SBL Top    Wind Speed                      [m/s] |
  ! |            TaT_SV   : SBL Top    Temperature                       [K] |
  ! |            rhT_SV   : SBL Top    Air  Density                  [kg/m3] |
  ! |            QaT_SV   : SBL Top    Specific  Humidity            [kg/kg] |
  ! |            LSdzsv   : Vertical   Discretization Factor             [-] |
  ! |                     =    1. Soil                                       |
  ! |                     = 1000. Ocean                                      |
  ! |            dzsnSV   : Snow Layer Thickness                         [m] |
  ! |            ro__SV   : Snow/Soil  Volumic Mass                  [kg/m3] |
  ! |            eta_SV   : Soil Water Content                       [m3/m3] |
  ! |            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] |
  ! |            Lx_H2O   : Latent Heat of Vaporization/Sublimation   [J/kg] |
  ! |            sEX_sv   : Verticaly Integrated Extinction Coefficient  [-] |
  ! |                                                                        |
  ! |   INPUT /  TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
  ! |   OUTPUT:           & Snow     Temperatures (layers  1,2,...,nsno) [K] |
  ! |   ^^^^^^                                                               |
  ! |                                                                        |
  ! |   OUTPUT:  IRs_SV   : Soil      IR Radiation                    [W/m2] |
  ! |   ^^^^^^   HSs_sv   : Sensible  Heat Flux                       [W/m2] |
  ! |            HLs_sv   : Latent    Heat Flux                       [W/m2] |
  ! |            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] |
  ! |                                                                        |
  ! |   Internal Variables:                                                  |
  ! |   ^^^^^^^^^^^^^^^^^^                                                   |
  ! |                                                                        |
  ! |   METHOD: NO   Skin Surface Temperature                                |
  ! |   ^^^^^^  Semi-Implicit Crank Nicholson Scheme                         |
  ! |                                                                        |
  ! | # OPTIONS: #E0: Energy Budget Verification                             |
  ! | # ^^^^^^^  #kd: KDsvat Option:NO Flux  Limitor     on HL               |
  ! | #          #KD: KDsvat Option:Explicit Formulation of HL               |
  ! | #          #NC: OUTPUT for Stand Alone NetCDF File                     |
  ! |                                                                        |
  ! +------------------------------------------------------------------------+




  ! +--Global Variables
  ! +  ================

  use VARphy
  use VAR_SV
  use VARdSV
  use VARxSV
  use VARySV
  use VARtSV
  use VAR0SV


  IMPLICIT NONE


  ! +--OUTPUT
  ! +  ------

  ! #e1 real      ETSo_0(knonv)                 ! Soil/Snow Power, before Forcing
  ! #e1 real      ETSo_1(knonv)                 ! Soil/Snow Power, after  Forcing
  ! #e1 real      ETSo_d(knonv)                 ! Soil/Snow Power, Forcing


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

  integer :: ikl   ,isl   ,jsl   ,ist      !
  integer :: ist__s,ist__w                 ! Soil/Water  Body Identifier
  integer :: islsgn                        ! Soil/Snow Surfac.Identifier
  real :: eps__3                        ! Arbitrary    Low Number
  real :: etaMid,psiMid                 ! Layer Interface's Humidity
  real :: mu_eta                        !     Soil thermal Conductivity
  real :: mu_exp                        ! arg Soil thermal Conductivity
  real :: mu_min                        ! Min Soil thermal Conductivity
  real :: mu_max                        ! Max Soil thermal Conductivity
  real :: mu_sno(knonv),mu_aux          !     Snow thermal Conductivity
  real :: mu__dz(knonv,-nsol:nsno+1)    ! mu_(eta,sno)   / dz
  real :: dtC_sv(knonv,-nsol:nsno)      ! dt      / C
  real :: IRs__D(knonv)                 ! UpwardIR Previous Iter.Contr.
  real :: dIRsdT(knonv)                 ! UpwardIR           T Derivat.
  real :: f_HSHL(knonv)                 ! Factor common to HS and HL
  real :: dRidTs(knonv)                 ! d(Rib)/d(Ts)
  real :: HS___D(knonv)                 ! Sensible Heat Flux Atm.Contr.
  real :: f___HL(knonv)                 !
  real :: HL___D(knonv)                 ! Latent   Heat Flux Atm.Contr.
  REAL :: TSurf0(knonv),dTSurf          ! Previous Surface Temperature
  real :: qsatsg(knonv) !,den_qs,arg_qs ! Soil   Saturat. Spec. Humidity
  real :: dqs_dT(knonv)                 ! d(qsatsg)/dTv
  real :: Psi(   knonv)                 ! 1st Soil Layer Water Potential
  real :: RHuSol(knonv)                 ! Soil Surface Relative Humidity
  real :: etaSol                        ! Soil Surface          Humidity
  real :: d__eta                        ! Soil Surface Humidity Increm.
  real :: Elem_A,Elem_C                 !   Diagonal Coefficients
  real :: Diag_A(knonv,-nsol:nsno)      ! A Diagonal
  real :: Diag_B(knonv,-nsol:nsno)      ! B Diagonal
  real :: Diag_C(knonv,-nsol:nsno)      ! C Diagonal
  real :: Term_D(knonv,-nsol:nsno)      !   Independant Term
  real :: Aux__P(knonv,-nsol:nsno)      ! P Auxiliary Variable
  real :: Aux__Q(knonv,-nsol:nsno)      ! Q Auxiliary Variable
  real :: Ts_Min,Ts_Max                 ! Temperature Limits
  ! #e1 real      Exist0                        ! Existing Layer Switch
  real :: psat_wat, psat_ice, sp        ! computation of qsat

  integer :: nt_srf,it_srf,itEuBk          ! HL: Surface Scheme
  parameter(nt_srf=10)                     ! 10 before
  real :: agpsrf,xgpsrf,dt_srf,dt_ver   !
  real :: etaBAK(knonv)                 !
  real :: etaNEW(knonv)                 !
  real :: etEuBk(knonv)                 !
  real :: fac_dt(knonv),faceta(knonv)   !
  real :: PsiArg(knonv),SHuSol(knonv)   !



  ! +--Internal DATA
  ! +  =============

  data      eps__3 /   1.e-3   /          ! Arbitrary    Low Number
  data      mu_exp /  -0.4343  /          ! Soil Thermal Conductivity
  data      mu_min /   0.172   /          ! Min Soil Thermal Conductivity
  data      mu_max /   2.000   /          ! Max Soil Thermal Conductivity
  data      Ts_Min / 175.      /          ! Temperature            Minimum
  data      Ts_Max / 300.      /          ! Temperature Acceptable Maximum
  ! +                                           ! including   Snow Melt  Energy

  ! +-- Initilialisation of local arrays
  ! +   ================================
    DO ikl=1,knonv

      mu_sno(ikl)=0.
      mu__dz(ikl,:)=0.
      dtC_sv(ikl,:)=0.
      IRs__D(ikl)=0.
      dIRsdT(ikl)=0.
      f_HSHL(ikl)=0.
      dRidTs(ikl)=0.
      HS___D(ikl)=0.
      f___HL(ikl)=0.
      HL___D(ikl)=0.
      TSurf0(ikl)=0.
      qsatsg(ikl)=0.
      dqs_dT(ikl)=0.
      Psi(ikl)=0.
      RHuSol(ikl)=0.
      Diag_A(ikl,:)=0.
      Diag_B(ikl,:)=0.
      Diag_C(ikl,:)=0.
      Term_D(ikl,:)=0.
      Aux__P(ikl,:)=0.
      Aux__Q(ikl,:)=0.
      etaBAK(ikl)=0.
      etaNEW(ikl)=0.
      etEuBk(ikl)=0.
      fac_dt(ikl)=0.
      faceta(ikl)=0.
      PsiArg(ikl)=0.
      SHuSol(ikl)=0.

    END DO



  ! +--Heat Conduction Coefficient (zero in the Layers over the highest one)
  ! +  ===========================
  ! +                             ---------------- isl    eta_SV, rho C (isl)
  ! +
  ! +--Soil                       ++++++++++++++++        etaMid,    mu (isl)
  ! +  ----
  ! +                             ---------------- isl-1  eta_SV, rho C (isl-1)
       isl=-nsol
    DO ikl=1,knonv

      mu__dz(ikl,isl) = 0.

      dtC_sv(ikl,isl) = dtz_SV2(isl)  * dt__SV & ! dt / (dz X rho C)
            /((rocsSV(isotSV(ikl)) & ! [s / (m.J/m3/K)]
            +rcwdSV*eta_SV(ikl,isl)) & !
            *LSdzsv(ikl)            )      !
    END DO
  DO   isl=-nsol+1,0
    DO ikl=1,knonv
      ist    =      isotSV(ikl)                       ! Soil Type
      ist__s =  min(ist, 1)                           ! 1 => Soil
      ist__w =  1 - ist__s                            ! 1 => Water Body

      etaMid = 0.5*(dz_dSV(isl-1)*eta_SV(ikl,isl-1) & ! eta at layers
            +dz_dSV(isl)  *eta_SV(ikl,isl)  ) & !     interface
            /dzmiSV(isl)                       ! LSdzsv implicit !
      etaMid =  max(etaMid,epsi)
      psiMid =      psidSV(ist) &
            *(etadSV(ist)/etaMid)**bCHdSV(ist)
      mu_eta =      3.82      *(psiMid)**mu_exp       ! Soil Thermal
      mu_eta =  min(max(mu_eta, mu_min), mu_max)      ! Conductivity
  ! +                                                       ! DR97 eq.3.31
      mu_eta =  ist__s *mu_eta +ist__w * vK_dSV       ! Water Bodies
  ! +                                                       ! Correction
      mu__dz(ikl,isl) = mu_eta/(dzmiSV(isl) & !
            *LSdzsv(ikl))          !

      dtC_sv(ikl,isl) = dtz_SV2(isl)* dt__SV & ! dt / (dz X rho C)
            /((rocsSV(isotSV(ikl)) & !
            +rcwdSV*eta_SV(ikl,isl)) & !
            *LSdzsv(ikl)            )      !
    END DO
  END DO


  ! +--Soil/Snow Interface
  ! +  -------------------

  ! +--Soil Contribution
  ! +  ^^^^^^^^^^^^^^^^^
       isl=1
    DO ikl=1,knonv
      ist    =      isotSV(ikl)                       ! Soil Type
      ist__s =  min(ist, 1)                           ! 1 => Soil
      ist__w =  1 - ist__s                            ! 1 => Water Body
      psiMid =      psidSV(ist)                       ! Snow => Saturation
      mu_eta =      3.82      *(psiMid)**mu_exp       ! Soil Thermal
      mu_eta =  min(max(mu_eta, mu_min), mu_max)      ! Conductivity
  ! +                                                       ! DR97 eq.3.31
      mu_eta =  ist__s *mu_eta +ist__w * vK_dSV       ! Water Bodies

  ! +--Snow Contribution
  ! +  ^^^^^^^^^^^^^^^^^
      mu_sno(ikl) =  CdidSV & !
            *(ro__SV(ikl,isl) /ro_Wat) ** 1.88 !
      mu_sno(ikl) =          max(epsi,mu_sno(ikl))    !
  ! +...    mu_sno :  Snow Heat Conductivity Coefficient [Wm/K]
  ! +                 (Yen 1981, CRREL Rep., 81-10)

  ! +--Combined Heat Conductivity
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
      mu__dz(ikl,isl) = 2./(dzsnSV(ikl,isl  ) & ! Combined Heat
            /mu_sno(ikl) & ! Conductivity
            +LSdzsv(ikl) & !
            *dz_dSV(    isl-1)/mu_eta) ! Coefficient

  ! +--Inverted Heat Capacity
  ! +  ^^^^^^^^^^^^^^^^^^^^^^
      dtC_sv(ikl,isl) = dt__SV/max(epsi, & ! dt / (dz X rho C)
            dzsnSV(ikl,isl) * ro__SV(ikl,isl) *Cn_dSV)      !
    END DO


  ! +--Snow
  ! +  ----

  DO ikl=1,knonv
  DO   isl=1,min(nsno,isnoSV(ikl)+1)
      ro__SV(ikl,isl) = & !
            ro__SV(ikl ,isl) & !
            * max(0,min(isnoSV(ikl)-isl+1,1))            !

    END DO
  END DO

  DO ikl=1,knonv
  DO   isl=1,min(nsno,isnoSV(ikl)+1)

  ! +--Combined Heat Conductivity
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
      mu_aux      =  CdidSV & !
            *(ro__SV(ikl,isl) /ro_Wat) ** 1.88 !
      mu__dz(ikl,isl) = & !
            2.                        *mu_aux*mu_sno(ikl) & ! Combined Heat
            /max(epsi,dzsnSV(ikl,isl  )*mu_sno(ikl) & ! Conductivity
            +dzsnSV(ikl,isl-1)*mu_aux     )       ! For upper Layer
      mu_sno(ikl)     =            mu_aux             !

  ! +--Inverted Heat Capacity
  ! +  ^^^^^^^^^^^^^^^^^^^^^^
      dtC_sv(ikl,isl) = dt__SV/max(eps__3, & ! dt / (dz X rho C)
            dzsnSV(ikl,isl) * ro__SV(ikl,isl) *Cn_dSV)      !
    END DO
  END DO


  ! +--Uppermost Effective Layer: NO conduction
  ! +  ----------------------------------------

    DO ikl=1,knonv
      mu__dz(ikl,isnoSV(ikl)+1) = 0.0
    END DO


  ! +--Energy Budget (IN)
  ! +  ==================

  ! #e1   DO ikl=1,knonv
  ! #e1     ETSo_0(ikl) = 0.
  ! #e1   END DO
  ! #e1 DO   isl= -nsol,nsno
  ! #e1   DO ikl=1,knonv
  ! #e1     Exist0      = isl -           isnoSV(ikl)
  ! #e1     Exist0      = 1.  - max(zero,min(unun,Exist0))
  ! #e1     ETSo_0(ikl) = ETSo_0(ikl)
  ! #e1.                +(TsisSV(ikl,isl)-TfSnow)*Exist0
  ! #e1.                                 /dtC_sv(ikl,isl)
  ! #e1   END DO
  ! #e1 END DO


  ! +--Tridiagonal Elimination: Set Up
  ! +  ===============================

  ! +--Soil/Snow Interior
  ! +  ^^^^^^^^^^^^^^^^^^
  DO ikl=1,knonv
  DO   isl=-nsol+1,min(nsno-1,isnoSV(ikl)+1)

      Elem_A          =  dtC_sv(ikl,isl)         *mu__dz(ikl,isl)
      Elem_C          =  dtC_sv(ikl,isl)         *mu__dz(ikl,isl+1)
      Diag_A(ikl,isl) = -Elem_A  *Implic
      Diag_C(ikl,isl) = -Elem_C  *Implic
      Diag_B(ikl,isl) =  1.0d+0  -Diag_A(ikl,isl)-Diag_C(ikl,isl)
      Term_D(ikl,isl) =  Explic *(Elem_A         *TsisSV(ikl,isl-1) &
            +Elem_C         *TsisSV(ikl,isl+1)) &
            +(1.0d+0 -Explic *(Elem_A+Elem_C))*TsisSV(ikl,isl) &
            + dtC_sv(ikl,isl)           * sol_SV(ikl)    *SoSosv(ikl) &
            *(sEX_sv(ikl,isl+1) &
            -sEX_sv(ikl,isl  ))
    END DO
  END DO

  ! +--Soil  lowest Layer
  ! +  ^^^^^^^^^^^^^^^^^^
       isl= -nsol
    DO ikl=1,knonv
      Elem_A          =  0.
      Elem_C          =  dtC_sv(ikl,isl)         *mu__dz(ikl,isl+1)
      Diag_A(ikl,isl) =  0.
      Diag_C(ikl,isl) = -Elem_C  *Implic
      Diag_B(ikl,isl) =  1.0d+0  -Diag_A(ikl,isl)-Diag_C(ikl,isl)
      Term_D(ikl,isl) =  Explic * Elem_C         *TsisSV(ikl,isl+1) &
            +(1.0d+0 -Explic * Elem_C)        *TsisSV(ikl,isl) &
            + dtC_sv(ikl,isl)           * sol_SV(ikl)    *SoSosv(ikl) &
            *(sEX_sv(ikl,isl+1) &
            -sEX_sv(ikl,isl  ))
    END DO

  ! +--Snow highest Layer (dummy!)
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^

    ! !EV!isl=  min(isnoSV(1)+1,nsno)

    DO ikl=1,knonv
  ! EV try to calculate isl at the ikl grid point
      isl=  min(isnoSV(ikl)+1,nsno)

      Elem_A          =  dtC_sv(ikl,isl)  *mu__dz(ikl,isl)
      Elem_C          =  0.
      Diag_A(ikl,isl) = -Elem_A  *Implic
      Diag_C(ikl,isl) =  0.
      Diag_B(ikl,isl) =  1.0d+0  -Diag_A(ikl,isl)
      Term_D(ikl,isl) =  Explic * Elem_A  *TsisSV(ikl,isl-1) &
            +(1.0d+0 -Explic * Elem_A) *TsisSV(ikl,isl) &
            + dtC_sv(ikl,isl) * (sol_SV(ikl)      *SoSosv(ikl) &
            *(sEX_sv(ikl,isl+1) &
            -sEX_sv(ikl,isl  )))
    END DO

  ! +--Surface: UPwardIR Heat Flux
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
    DO ikl=1,knonv
      isl         = isnoSV(ikl)
      dIRsdT(ikl) = Eso_sv(ikl)* StefBo          * 4. & ! - d(IR)/d(T)
            * TsisSV(ikl,isl) & !
            * TsisSV(ikl,isl) & !
            * TsisSV(ikl,isl)           !
      IRs__D(ikl) = dIRsdT(ikl)* TsisSV(ikl,isl) * 0.75    !

  ! +--Surface: Richardson Number:   T Derivative
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  ! #RC     dRidTs(ikl) =-gravit      *    za__SV(ikl)
  ! #RC.                /(TaT_SV(ikl) *    VV__SV(ikl)
  ! #RC.                              *    VV__SV(ikl))

  ! +--Surface: Turbulent Heat Flux: Factors
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      f_HSHL(ikl) = rhT_SV(ikl)  /    rah_sv(ikl)           ! to  HS, HL
      f___HL(ikl) = f_HSHL(ikl) *    Lx_H2O(ikl)

  ! +--Surface: Sensible  Heat Flux: T Derivative
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      dSdTSV(ikl) = f_HSHL(ikl) *    Cp                    !#- d(HS)/d(T)
  ! #RC.         *(1.0  -(TsisSV(ikl,isl) -TaT_SV(ikl))          !#Richardson
  ! #RC.         * dRidTs(ikl)*dFh_sv(ikl)/rah_sv(ikl))          ! Nb. Correct.
      HS___D(ikl) = dSdTSV(ikl) *    TaT_SV(ikl)           !

  ! +--Surface: Latent    Heat Flux: Saturation Specific Humidity
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      ! den_qs      =         TsisSV(ikl,isl)- 35.8          !
      ! arg_qs      = 17.27 *(TsisSV(ikl,isl)-273.16)        !
  !    .                                   / den_qs              !
  !         qsatsg(ikl) = .0038 *        exp(arg_qs)             !
      !  sp = (pst_SV(ikl) + ptopSV) * 10.

      ! !sp=ps__SV(ikl)
      ! ! Etienne: in the formula herebelow sp should be in hPa, not
      ! ! in Pa so I divide by 100.
      sp=ps__SV(ikl)/100.
      psat_ice = 6.1070 * exp(6150. *(1./273.16 - &
            1./TsisSV(ikl,isl)))

      psat_wat = 6.1078 * exp (5.138*log(273.16   /TsisSV(ikl,isl))) &
            * exp (6827.*(1.         /273.16-1./TsisSV(ikl,isl)))

      if(TsisSV(ikl,isl)<=273.16) then
        qsatsg(ikl) = 0.622 * psat_ice / (sp - 0.378 * psat_ice)
      else
        qsatsg(ikl) = 0.622 * psat_wat / (sp - 0.378 * psat_wat)
      endif
      QsT_SV(ikl)=qsatsg(ikl)

      ! dqs_dT(ikl) = qsatsg(ikl)* 4099.2   /(den_qs *den_qs)!
      fac_dt(ikl) = f_HSHL(ikl)/(ro_Wat   * dz_dSV(0))     !
    END DO



  ! +--Surface: Latent    Heat Flux: Surface    Relative Humidity
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          xgpsrf       =   1.05                            !
          agpsrf       = dt__SV*(   1.0-xgpsrf        ) & !
                /(   1.0-xgpsrf**nt_srf)    !
          dt_srf       = agpsrf                            !
          dt_ver       = 0.

        DO ikl=1,knonv
          isl          =          isnoSV(ikl)
          ist          = max(0,isotSV(ikl)-100*isnoSV(ikl))! 0 if    H2O
          ist__s       = min(1,ist)
          etaBAK(ikl)  = max(epsi,eta_SV(ikl ,isl))        !
          etaNEW(ikl)  =          etaBAK(ikl)              !
          etEuBk(ikl)  =          etaNEW(ikl)              !
        END DO

    if(ist__s==1) then ! to reduce computer time
                                      ! !
    DO it_srf=1,nt_srf                                     !
          dt_ver       = dt_ver     +dt_srf                !
        DO ikl=1,knonv                                     !
          faceta(ikl)  = fac_dt(ikl)*dt_srf                !
  ! #VX         faceta(ikl)  = faceta(ikl)                       !
  ! #VX.                  /(1.+faceta(ikl)*dQa_SV(ikl))          !    Limitation
  !                                                              ! by Atm.Conten
  ! #??.        *max(0,sign(1.,qsatsg(ikl)-QaT_SV(ikl))))        ! NO Limitation
                                                           ! ! of Downw.Flux
        END DO                                             !
      DO itEuBk=1,2                                        !
        DO ikl=1,knonv
          ist    = max(0,isotSV(ikl)-100*isnoSV(ikl))      ! 0 if    H2O
                                                           ! !
          Psi(ikl) = & !
                psidSV(ist) & ! DR97, Eqn 3.34
                *(etadSV(ist) & !
                /max(etEuBk(ikl),epsi)) & !
                **bCHdSV(ist)                              !
          PsiArg(ikl) = 7.2E-5*Psi(ikl)                    !
          RHuSol(ikl) =   exp(-min(0.,PsiArg(ikl)))    !
          SHuSol(ikl) =     qsatsg(ikl)  *RHuSol(ikl)      ! DR97, Eqn 3.15
          etEuBk(ikl) = & !
                (etaNEW(ikl) + faceta(ikl)*(QaT_SV(ikl) & !
                -SHuSol(ikl) & !
                *(1.          -bCHdSV(ist) & !
                *PsiArg(ikl))       )) & !
                /(1.          + faceta(ikl)* SHuSol(ikl) & !
                *bCHdSV(ist) & !
                *PsiArg(ikl) & !
                /etaNEW(ikl))          !
          etEuBk(ikl) = etEuBk(ikl) & !
  !    .                                 /(Ro_Wat*dz_dSV(0))     !
                * dt_srf     /(Ro_Wat*dz_dSV(0))     !
  !XF 15/05/2017 BUG
        END DO                                             !
      END DO                                               !
        DO ikl=1,knonv                                     !
          etaNEW(ikl) =  max(etEuBk(ikl),epsi)             !
        END DO                                             !
          dt_srf      =      dt_srf         * xgpsrf       !
    END DO


    endif                                       !

  ! +--Surface: Latent    Heat Flux: Soil/Water Surface Contributions
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    DO ikl=1,knonv                                         !
      isl        =  isnoSV(ikl)                            !
      ist   = max(0,isotSV(ikl)-100*isnoSV(ikl))           ! 0 if    H2O
      ist__s= min(1,ist)                                   ! 1 if no H2O
      ist__w=     1-ist__s                                 ! 1 if    H2O
      d__eta     =  eta_SV(ikl,isl)-etaNEW(ikl)            !
      ! ! latent heat flux computation
      HL___D(ikl)=( ist__s *ro_Wat *dz_dSV(0) & ! Soil Contrib.
            *(etaNEW(ikl)    -etaBAK(ikl)) / dt__SV & !
            +ist__w         *f_HSHL(ikl) & ! H2O  Contrib.
            *(QaT_SV(ikl)    - qsatsg(ikl))        ) & !
            * Lx_H2O(ikl)                            ! common factor

  ! #DL     RHuSol(ikl) =(QaT_SV(ikl)                            !
  ! #DL.                 -HL___D(ikl)    / f___HL(ikl))          !
  ! #DL.                / qsatsg(ikl)                            !

  ! +--Surface: Latent    Heat Flux: T Derivative
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      dLdTSV(ikl) = 0.
  ! #DL     dLdTSV(ikl) = f___HL(ikl) * RHuSol(ikl) *dqs_dT(ikl) ! - d(HL)/d(T)
  ! #DL     HL___D(ikl) = HL___D(ikl)                            !
  ! #DL.                 +dLdTSV(ikl) * TsisSV(ikl,isl)          !
    END DO                                                 !

  ! +--Surface: Tridiagonal Matrix Set Up
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    DO ikl=1,knonv
      isl             =  isnoSV(ikl)
      TSurf0(ikl)     =  TsisSV(ikl,isl)

      Elem_A          =  dtC_sv(ikl,isl)*mu__dz(ikl,isl)
      Elem_C          =  0.
      Diag_A(ikl,isl) = -Elem_A *Implic
      Diag_C(ikl,isl) =  0.
      Diag_B(ikl,isl) =  1.0d+0 -Diag_A(ikl,isl)
      Diag_B(ikl,isl) =  Diag_B(ikl,isl) &
            + dtC_sv(ikl,isl) * (dIRsdT(ikl) & ! Upw. Sol IR
            +dSdTSV(ikl) & ! HS/Surf.Contr.
            +dLdTSV(ikl))                      ! HL/Surf.Contr.

      Term_D(ikl,isl) =  Explic *Elem_A *TsisSV(ikl,isl-1) &
            +(1.0d+0 -Explic *Elem_A)*TsisSV(ikl,isl)



      Term_D(ikl,isl) =  Term_D(ikl,isl) &
            + dtC_sv(ikl,isl) * (sol_SV(ikl)    *SoSosv(ikl) & ! Absorbed
            *(sEX_sv(ikl,isl+1) & ! Solar
            -sEX_sv(ikl,isl  )) & !
            +   IRd_SV(ikl)*Eso_sv(ikl) & ! Down Atm IR
            +IRs__D(ikl) & ! Upw. Sol IR
            +HS___D(ikl) & ! HS/Atmo.Contr.
            +HL___D(ikl)            )! HL/Atmo.Contr.

    END DO


  ! +--Tridiagonal Elimination
  ! +  =======================

  ! +--Forward  Sweep
  ! +  ^^^^^^^^^^^^^^
    DO ikl=  1,knonv
      Aux__P(ikl,-nsol) = Diag_B(ikl,-nsol)
      Aux__Q(ikl,-nsol) =-Diag_C(ikl,-nsol)/Aux__P(ikl,-nsol)
    END DO

    DO ikl=      1,knonv

    DO   isl=-nsol+1,min(nsno,isnoSV(ikl)+1)
      Aux__P(ikl,isl)   = Diag_A(ikl,isl)  *Aux__Q(ikl,isl-1) &
            +Diag_B(ikl,isl)
      Aux__Q(ikl,isl)   =-Diag_C(ikl,isl)  /Aux__P(ikl,isl)
    END DO
    END DO

    DO ikl=      1,knonv
      TsisSV(ikl,-nsol) = Term_D(ikl,-nsol)/Aux__P(ikl,-nsol)
    END DO

    DO ikl=      1,knonv
    DO   isl=-nsol+1,min(nsno,isnoSV(ikl)+1)
      TsisSV(ikl,isl)   =(Term_D(ikl,isl) &
            -Diag_A(ikl,isl)  *TsisSV(ikl,isl-1)) &
            /Aux__P(ikl,isl)


    END DO
    END DO

  ! +--Backward Sweep
  ! +  ^^^^^^^^^^^^^^
    DO ikl=     1,knonv
    DO   isl=min(nsno-1,isnoSV(ikl)+1),-nsol,-1


      TsisSV(ikl,isl)   = Aux__Q(ikl,isl)  *TsisSV(ikl,isl+1) &
            +TsisSV(ikl,isl)
      if(isl==0.and.isnoSV(ikl)==0) then

       TsisSV(ikl,isl)  = min(TaT_SV(ikl)+30,TsisSV(ikl,isl))
       TsisSV(ikl,isl)  = max(TaT_SV(ikl)-30,TsisSV(ikl,isl))


  ! #EU      TsisSV(ikl,isl)  = max(TaT_SV(ikl)-15.,TsisSV(ikl,isl))

      ! !XF 18/11/2018 to avoid ST reaching 70�C!!
      ! !It is an error compensation but does not work over tundra

      endif



    END DO

  END DO



  ! +--Temperature Limits (avoids problems in case of no Snow Layers)
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    DO ikl=     1,knonv
       isl              = isnoSV(ikl)

       dTSurf            = TsisSV(ikl,isl) -     TSurf0(ikl)
      TsisSV(ikl,isl)   = TSurf0(ikl) + sign(1.,dTSurf) & ! 180.0 dgC/hr
            * min(abs(dTSurf),5.e-2*dt__SV)         ! =0.05 dgC/s



    END DO

    DO ikl=     1,knonv
    DO   isl=min(nsno,isnoSV(ikl)+1),1      ,-1
      TsisSV(ikl,isl)   = max(Ts_Min,       TsisSV(ikl,isl))
      TsisSV(ikl,isl)   = min(Ts_Max,       TsisSV(ikl,isl))
    END DO

    END DO

  ! +--Update Surface    Fluxes
  ! +  ========================



    DO ikl=      1,knonv
      isl         = isnoSV(ikl)
      IRs_SV(ikl) = IRs__D(ikl) & !
            - dIRsdT(ikl) * TsisSV(ikl,isl)        !
      HSs_sv(ikl) = HS___D(ikl) & ! Sensible Heat
            - dSdTSV(ikl) * TsisSV(ikl,isl)        ! Downward > 0
      HLs_sv(ikl) = HL___D(ikl) & ! Latent   Heat
            - dLdTSV(ikl) * TsisSV(ikl,isl)        ! Downward > 0
    END DO

  return
end subroutine sisvat_tso
