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

      IMPLICIT NONE

 
C +--OUTPUT
C +  ------
 
! #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
 
 
C +--Internal Variables
C +  ==================
 
      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)   !
 

 
C +--Internal DATA
C +  =============
 
      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
C +                                           ! including   Snow Melt  Energy
 
C +-- Initilialisation of local arrays
C +   ================================
        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

        

C +--Heat Conduction Coefficient (zero in the Layers over the highest one)
C +  ===========================
C +                             ---------------- isl    eta_SV, rho C (isl)
C +
C +--Soil                       ++++++++++++++++        etaMid,    mu (isl)
C +  ----
C +                             ---------------- 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
C +                                                       ! DR97 eq.3.31
          mu_eta =  ist__s *mu_eta +ist__w * vK_dSV       ! Water Bodies
C +                                                       ! 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
 
 
C +--Soil/Snow Interface
C +  -------------------
 
C +--Soil Contribution
C +  ^^^^^^^^^^^^^^^^^
           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
C +                                                       ! DR97 eq.3.31
          mu_eta =  ist__s *mu_eta +ist__w * vK_dSV       ! Water Bodies
 
C +--Snow Contribution
C +  ^^^^^^^^^^^^^^^^^
          mu_sno(ikl) =  CdidSV                           !
     .                 *(ro__SV(ikl,isl) /ro_Wat) ** 1.88 !
          mu_sno(ikl) =          max(epsi,mu_sno(ikl))    !
C +...    mu_sno :  Snow Heat Conductivity Coefficient [Wm/K]
C +                 (Yen 1981, CRREL Rep., 81-10)
 
C +--Combined Heat Conductivity
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
          mu__dz(ikl,isl) = 2./(dzsnSV(ikl,isl  )         ! Combined Heat
     .                         /mu_sno(ikl)               ! Conductivity
     .                         +LSdzsv(ikl)               !
     .                         *dz_dSV(    isl-1)/mu_eta) ! Coefficient
 
C +--Inverted Heat Capacity
C +  ^^^^^^^^^^^^^^^^^^^^^^
          dtC_sv(ikl,isl) = dt__SV/max(epsi,              ! dt / (dz X rho C)
     .    dzsnSV(ikl,isl) * ro__SV(ikl,isl) *Cn_dSV)      !
        END DO
 
 
C +--Snow
C +  ----
 
      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)
 
C +--Combined Heat Conductivity
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
          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             !
 
C +--Inverted Heat Capacity
C +  ^^^^^^^^^^^^^^^^^^^^^^
          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
 
 
C +--Uppermost Effective Layer: NO conduction
C +  ----------------------------------------
 
        DO ikl=1,knonv
          mu__dz(ikl,isnoSV(ikl)+1) = 0.0
        END DO
 
 
C +--Energy Budget (IN)
C +  ==================
 
! #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
 
 
C +--Tridiagonal Elimination: Set Up
C +  ===============================
 
C +--Soil/Snow Interior
C +  ^^^^^^^^^^^^^^^^^^
      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
 
C +--Soil  lowest Layer
C +  ^^^^^^^^^^^^^^^^^^
           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
 
C +--Snow highest Layer (dummy!)
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^

        !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
 
C +--Surface: UPwardIR Heat Flux
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
        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    !
 
C +--Surface: Richardson Number:   T Derivative
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c #RC     dRidTs(ikl) =-gravit      *    za__SV(ikl)
c #RC.                /(TaT_SV(ikl) *    VV__SV(ikl)
c #RC.                              *    VV__SV(ikl))
 
C +--Surface: Turbulent Heat Flux: Factors
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          f_HSHL(ikl) = rhT_SV(ikl)  /    rah_sv(ikl)           ! to  HS, HL
          f___HL(ikl) = f_HSHL(ikl) *    Lx_H2O(ikl)
 
C +--Surface: Sensible  Heat Flux: T Derivative
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          dSdTSV(ikl) = f_HSHL(ikl) *    Cp                    !#- d(HS)/d(T)
c #RC.         *(1.0  -(TsisSV(ikl,isl) -TaT_SV(ikl))          !#Richardson
c #RC.         * dRidTs(ikl)*dFh_sv(ikl)/rah_sv(ikl))          ! Nb. Correct.
          HS___D(ikl) = dSdTSV(ikl) *    TaT_SV(ikl)           !
 
C +--Surface: Latent    Heat Flux: Saturation Specific Humidity
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c         den_qs      =         TsisSV(ikl,isl)- 35.8          !
c         arg_qs      = 17.27 *(TsisSV(ikl,isl)-273.16)        !
c    .                                   / den_qs              !
c         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)

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


 
C +--Surface: Latent    Heat Flux: Surface    Relative Humidity
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
              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                !
c #VX         faceta(ikl)  = faceta(ikl)                       !
c #VX.                  /(1.+faceta(ikl)*dQa_SV(ikl))          !    Limitation
                                                               ! by Atm.Conten
c #??.        *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)          !
c    .                                 /(Ro_Wat*dz_dSV(0))     !
     .                    * dt_srf     /(Ro_Wat*dz_dSV(0))     !
cXF 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                                       !
 
C +--Surface: Latent    Heat Flux: Soil/Water Surface Contributions
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        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
 
c #DL     RHuSol(ikl) =(QaT_SV(ikl)                            !
c #DL.                 -HL___D(ikl)    / f___HL(ikl))          !
c #DL.                / qsatsg(ikl)                            !
 
C +--Surface: Latent    Heat Flux: T Derivative
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          dLdTSV(ikl) = 0.
c #DL     dLdTSV(ikl) = f___HL(ikl) * RHuSol(ikl) *dqs_dT(ikl) ! - d(HL)/d(T)
c #DL     HL___D(ikl) = HL___D(ikl)                            !
c #DL.                 +dLdTSV(ikl) * TsisSV(ikl,isl)          !
        END DO                                                 !
 
C +--Surface: Tridiagonal Matrix Set Up
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        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
 
 
C +--Tridiagonal Elimination
C +  =======================
 
C +--Forward  Sweep
C +  ^^^^^^^^^^^^^^
        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
 
C +--Backward Sweep
C +  ^^^^^^^^^^^^^^
        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))


c #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


 
C +--Temperature Limits (avoids problems in case of no Snow Layers)
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        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
 
C +--Update Surface    Fluxes
C +  ========================
        


        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
