

subroutine SISVAT_qSo
  ! #m0.                     (Wats_0,Wats_1,Wats_d)

  ! +------------------------------------------------------------------------+
  ! | MAR          SISVAT_qSo                                 6-04-2001  MAR |
  ! |   SubRoutine SISVAT_qSo computes the Soil      Water  Balance          |
  ! +------------------------------------------------------------------------+
  ! |                                                                        |
  ! |   PARAMETERS:  knonv: Total Number of columns =                        |
  ! |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
  ! |                     X       Number of Mosaic Cell per grid box         |
  ! |                                                                        |
  ! |   INPUT:   isnoSV   = total Nb of Ice/Snow Layers                      |
  ! |   ^^^^^    isotSV   = 0,...,11:   Soil       Type                      |
  ! |                       0:          Water, Solid or Liquid               |
  ! |                                                                        |
  ! |   INPUT:   rhT_SV   : SBL Top    Air  Density                  [kg/m3] |
  ! |   ^^^^^    drr_SV   : Rain   Intensity                       [kg/m2/s] |
  ! |            LSdzsv   : Vertical   Discretization Factor             [-] |
  ! |                     =    1. Soil                                       |
  ! |                     = 1000. Ocean                                      |
  ! |            dt__SV   : Time   Step                                  [s] |
  ! |                                                                        |
  ! |            Lx_H2O   : Latent Heat of Vaporization/Sublimation   [J/kg] |
  ! |            HLs_sv   : Latent Heat  Flux                         [W/m2] |
  ! |                                                                        |
  ! |   INPUT /  eta_SV   : Water      Content                       [m3/m3] |
  ! |   OUTPUT:  Khydsv   : Soil   Hydraulic    Conductivity           [m/s] |
  ! |   ^^^^^^                                                               |
  ! |                                                                        |
  ! |   OUTPUT:  RnofSV   : RunOFF Intensity                       [kg/m2/s] |
  ! |   ^^^^^^   Wats_0   : Soil Water,  before Forcing                 [mm] |
  ! |            Wats_1   : Soil Water,  after  Forcing                 [mm] |
  ! |            Wats_d   : Soil Water          Forcing                 [mm] |
  ! |                                                                        |
  ! |   Internal Variables:                                                  |
  ! |   ^^^^^^^^^^^^^^^^^^                                                   |
  ! |            z_Bump   : (Partly)Bumpy Layers Height                  [m] |
  ! |            z0Bump   :         Bumpy Layers Height                  [m] |
  ! |            dzBump   :  Lowest Bumpy Layer:                         [m] |
  ! |            etBump   :         Bumps Layer Averaged Humidity    [m3/m3] |
  ! |            etaMid   : Layer Interface's Humidity               [m3/m3] |
  ! |            eta__f   : Layer             Humidity  (Water Front)[m3/m3] |
  ! |            Dhyd_f   : Soil  Hydraulic Diffusivity (Water Front) [m2/s] |
  ! |            Dhydif   : Soil  Hydraulic Diffusivity               [m2/s] |
  ! |            WgFlow   : Water         gravitational     Flux   [kg/m2/s] |
  ! |            Wg_MAX   : Water MAXIMUM gravitational     Flux   [kg/m2/s] |
  ! |            SatRat   : Water         Saturation        Flux   [kg/m2/s] |
  ! |            WExces   : Water         Saturation Excess Flux   [kg/m2/s] |
  ! |            Dhydtz   : Dhydif * dt / dz                             [m] |
  ! |            FreeDr   : Free Drainage Fraction                       [-] |
  ! |            Elem_A   : A Diagonal Coefficient                           |
  ! |            Elem_C   : C Diagonal Coefficient                           |
  ! |            Diag_A   : A Diagonal                                       |
  ! |            Diag_B   : B Diagonal                                       |
  ! |            Diag_C   : C Diagonal                                       |
  ! |            Term_D   :   Independant Term                               |
  ! |            Aux__P   : P Auxiliary Variable                             |
  ! |            Aux__Q   : Q Auxiliary Variable                             |
  ! |                                                                        |
  ! |   TUNING PARAMETER:                                                    |
  ! |   ^^^^^^^^^^^^^^^^                                                     |
  ! |            z0soil   : Soil Surface averaged Bumps Height           [m] |
  ! |                                                                        |
  ! |   METHOD: NO   Skin Surface Humidity                                   |
  ! |   ^^^^^^  Semi-Implicit Crank Nicholson Scheme                         |
  ! |           (Partial) free Drainage, Water Bodies excepted (Lakes, Sea)  |
  ! |                                                                        |

  ! |                                                                        |
  ! | # OPTIONS: #GF: Saturation Front                                       |
  ! | # ^^^^^^^  #GH: Saturation Front allows Horton Runoff                  |
  ! | #          #GA: Soil Humidity Geometric Average                        |
  ! | #          #BP: Parameterization of Terrain Bumps                      |
  ! |                                                                        |
  ! |                                                                        |
  ! +------------------------------------------------------------------------+




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

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


  IMPLICIT NONE


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

  ! Water (Mass) Budget
  ! ~~~~~~~~~~~~~~~~~~~
  ! #m0 real      Wats_0(knonv)                 ! Soil Water,  before forcing
  ! #m0 real      Wats_1(knonv)                 ! Soil Water,  after  forcing
  ! #m0 real      Wats_d(knonv)                 ! Soil Water          forcing


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

  integer :: isl   ,jsl   ,ist   ,ikl      !
  integer :: ikm   ,ikp   ,ik0   ,ik1      !
  integer :: ist__s,ist__w                 ! Soil/Water Body Identifier
  ! #BP real      z0soil                        ! Soil Surface Bumps Height  [m]
  ! #BP real      z_Bump                        !(Partly)Bumpy Layers Height [m]
  ! #BP real      z0Bump                        !        Bumpy Layers Height [m]
  ! #BP real      dzBump                        ! Lowest Bumpy Layer:

  ! #BP real      etBump(knonv)                 ! Bumps Layer Averaged Humidity
  real :: etaMid                        ! Layer Interface's Humidity
  real :: Dhydif                        ! Hydraulic Diffusivity   [m2/s]
  real :: eta__f                        ! Water Front Soil Water Content
  real :: Khyd_f                        ! Water Front Hydraulic Conduct.
  real :: Khydav                        ! Hydraulic Conductivity   [m/s]
  real :: WgFlow                        ! Water gravitat. Flux [kg/m2/s]
  real :: Wg_MAX                        ! Water MAX.grav. Flux [kg/m2/s]
  real :: SatRat                        ! Saturation      Flux [kg/m2/s]
  real :: WExces                        ! Saturat. Excess Flux [kg/m2/s]
  real :: SoRnOF(knonv)                 ! Soil     Run    OFF
  real :: Dhydtz(knonv,-nsol:0)         ! Dhydif * dt / dz           [m]
  real :: Elem_A,Elem_B,Elem_C          !   Diagonal Coefficients
  real :: Diag_A(knonv,-nsol:0)         ! A Diagonal
  real :: Diag_B(knonv,-nsol:0)         ! B Diagonal
  real :: Diag_C(knonv,-nsol:0)         ! C Diagonal
  real :: Term_D(knonv,-nsol:0)         !   Independant Term
  real :: Aux__P(knonv,-nsol:0)         ! P Auxiliary Variable
  real :: Aux__Q(knonv,-nsol:0)         ! Q Auxiliary Variable
  real :: etaaux(knonv,-nsol:-nsol+1)   ! Soil Water Content     [m3/m3]
  real :: FreeDr                        ! Free Drainage Fraction (actual)
  real :: FreeD0                        ! Free Drainage Fraction (1=Full)
  real :: aKdtSV3( 0:nsot, 0:nkhy)      ! Khyd=a*eta+b: a * dt
  real :: bKdtSV3( 0:nsot, 0:nkhy)      ! Khyd=a*eta+b: b * dt

  ! Water (Mass) Budget
  ! ~~~~~~~~~~~~~~~~~~~
  ! #mw logical         mwopen                  ! IO   Switch
  ! #mw common/Sm_qSo_L/mwopen                  !
  ! #mw real     hourwr,timewr                  !
  ! #mw common/Sm_qSo_R/timewr                  !
  ! #mw real            Evapor(knonv)           !


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

  ! #BP data      z0soil/0.020/                 ! Soil Surface Bumps Height  [m]
  data      FreeD0/1.000/                 ! Free Drainage Fraction (1=Full)

  aKdtSV3=aKdtSV2*dt__SV
  bKdtSV3=bKdtSV2*dt__SV

  ! Water  Budget (IN)
  ! ==================

  ! #m0   DO ikl=1,knonv
  ! #m0     Wats_0(ikl) = 0.                    ! OLD RunOFF Contrib.
  ! #m0     Wats_d(ikl) = drr_SV(ikl)           ! Water Surface Forc.
  ! #m0   END DO

  ! #m0      isl= -nsol
  ! #m0   DO ikl=1,knonv
  ! #m0     Wats_0(ikl) = Wats_0(ikl)
  ! #m0.      + ro_Wat *( eta_SV(ikl,isl)   *dz78SV(isl)
  ! #m0.                + eta_SV(ikl,isl+1) *dz_8SV(isl) ) * LSdzsv(ikl)
  ! #m0   END DO

  ! #m0 DO   isl= -nsol+1,-1
  ! #m0   DO ikl=1,knonv
  ! #m0     Wats_0(ikl) = Wats_0(ikl)
  ! #m0.      + ro_Wat *( eta_SV(ikl,isl)   *dz34SV(isl)
  ! #m0.                +(eta_SV(ikl,isl-1)
  ! #m0.                 +eta_SV(ikl,isl+1))*dz_8SV(isl) ) * LSdzsv(ikl)
  ! #m0   END DO
  ! #m0 END DO

  ! #m0      isl=  0
  ! #m0   DO ikl=1,knonv
  ! #m0     Wats_0(ikl) = Wats_0(ikl)
  ! #m0.      + ro_Wat *( eta_SV(ikl,isl)   *dz78SV(isl)
  ! #m0.                + eta_SV(ikl,isl-1) *dz_8SV(isl) ) * LSdzsv(ikl)
  ! #m0   END DO


  ! +--Gravitational Flow
  ! +  ==================

  ! +...    METHOD: Surface Water Flux saturates successively the soil layers
  ! +       ^^^^^^  from up to below, but is limited by infiltration capacity.
  ! +               Hydraulic Conductivity again contributes after this step,
  ! +               not redundantly because of a constant (saturated) profile.

  ! +--Flux  Limitor
  ! +  ^^^^^^^^^^^^^
       isl=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
      Dhydif = s1__SV(ist) &
            *max(epsi,eta_SV(ikl,isl)) & ! Hydraulic Diffusivity
            **(bCHdSV(ist)+2.)     ! DR97, Eqn.(3.36)
      Dhydif = ist__s    * Dhydif & !
            + ist__w    * vK_dSV              ! Water Bodies
  ! +
      Khydav = ist__s    * Ks_dSV(ist) & ! DR97  Assumption
            + ist__w    * vK_dSV              ! Water Bodies
  ! +
      Wg_MAX = ro_Wat     *Dhydif & ! MAXimum  Infiltration
            *(etadSV(ist)-eta_SV(ikl,isl)) & !          Rate
            /(dzAvSV(isl)*LSdzsv(ikl)    ) & !
            +  ro_Wat     *Khydav              !

  ! +--Surface Horton RunOFF
  ! +  ^^^^^^^^^^^^^^^^^^^^^
      SoRnOF(ikl) = &
            max(zero,drr_SV(ikl)-Wg_MAX)
    RuofSV(ikl,1) = RuofSV(ikl,1) +    SoRnOF(ikl)
      drr_SV(ikl) =        drr_SV(ikl)-SoRnOF(ikl)
    RuofSV(ikl,2) = RuofSV(ikl,2) +max(0.,drr_SV(ikl))
    END DO

  ! #GF DO   isl=0,-nsol,-1
  ! #GF   DO ikl=1,knonv
  ! #GF     ist    = isotSV(ikl)                     ! Soil Type
  ! #GF     ist__s = min(ist, 1)                     ! 1 => Soil
  ! #GF     ist__w = 1 - ist__s                      ! 1 => Water Body

  ! +--Water Diffusion
  ! +  ^^^^^^^^^^^^^^^
  ! #GF     Dhydif = s1__SV(ist)
  ! #GF.               *max(epsi,eta_SV(ikl,isl))    ! Hydraulic Diffusivity
  ! #GF.                      **(bCHdSV(ist)+2.)     ! DR97, Eqn.(3.36)
  ! #GF     Dhydif = ist__s    * Dhydif              !
  ! #GF.           + ist__w    * vK_dSV              ! Water Bodies

  ! +--Water Conduction (without Horton Runoff)
  ! +  ^^^^^^^^^^^^^^^^
  ! #GF     Khyd_f =             Ks_dSV(ist)
  ! +...    Uses saturated K ==> Horton Runoff ~0    !

  ! +--Water Conduction (with    Horton Runoff)
  ! +  ^^^^^^^^^^^^^^^^
  ! #GH     ik0    = nkhy       *eta_SV(ikl,isl)
  ! #GH.                        /etadSV(ist)
  ! #GH     eta__f         =            1.
  ! #GH.   -aKdtSV3(ist,ik0)/(2. *dzAvSV(isl)
  ! #GH.                        *LSdzsv(ikl))
  ! #GH     eta__f         = max(eps_21,eta__f)
  ! #GH     eta__f         = min(etadSV(ist),
  ! #GH.                         eta_SV(ikl,isl) +
  ! #GH.   (aKdtSV3(ist,ik0)     *eta_SV(ikl,isl)
  ! #GH.   +bKdtSV3(ist,ik0))   /(dzAvSV(isl)
  ! #GH.                        *LSdzsv(ikl))
  ! #GH.                       / eta__f          )
  ! #GH     eta__f         = .5*(eta_SV(ikl,isl)
  ! #GH.                        +eta__f)

  ! #gh     eta__f         =     eta_SV(ikl,isl)

  ! #GH     ik0    = nkhy       *eta__f
  ! #GH.                        /etadSV(ist)
  ! #GH     Khyd_f =
  ! #GH.   (aKdtSV3(ist,ik0)     *eta__f
  ! #GH.   +bKdtSV3(ist,ik0))    /dt__SV

  ! #GF     Khydav = ist__s    * Khyd_f              ! DR97  Assumption
  ! #GF.           + ist__w    * vK_dSV              ! Water Bodies

  ! +--Gravitational Flow
  ! +  ^^^^^^^^^^^^^^^^^^
  ! #GF     Wg_MAX =                                 ! MAXimum  Infiltration
  ! #GF.             ro_Wat     *Dhydif              !          Rate
  ! #GF.           *(etadSV(ist)-eta_SV(ikl,isl))    !
  ! #GF.           /(dzAvSV(isl)*LSdzsv(ikl)    )    !
  ! #GF.          +  ro_Wat     *Khydav              !
  ! #GF   END DO
  ! #GF END DO
  ! #GF   DO ikl=1,knonv
  ! #GF     SoRnOF(ikl)     =    SoRnOF(ikl)         ! RunOFF Intensity
  ! #GF.                    +    drr_SV(ikl)         ! [kg/m2/s]
  ! +!!!    Inclure la possibilite de creer une mare sur un bedrock impermeable
  ! #GF     drr_SV(ikl) = 0.
  ! #GF   END DO


  ! +--Temperature Correction due to a changed Soil Energy Content
  ! +  ===========================================================

  ! +!!!    Mettre en oeuvre le couplage humidit?-?nergie


  ! +--Full Resolution of the Richard's Equation
  ! +  =========================================

  ! +...    METHOD: Water content evolution results from water fluxes
  ! +       ^^^^^^  at the layer boundaries
  ! +               Conductivity is approximated by a piecewise linear profile.
  ! +               Semi-Implicit Crank-Nicholson scheme is used.
  ! +              (Bruen, 1997, Sensitivity of hydrological processes
  ! +                            at the land-atmosphere interface.
  ! +                            Proc. Royal Irish Academy,  IGBP symposium
  ! +                            on global change and the Irish Environment.
  ! +                            Publ.: Maynooth)

  ! +                      - - - - - - - -   isl+1/2   - -  ^
  ! +                                                       |
  ! +   eta_SV(isl)        ---------------   isl     -----  +--dz_dSV(isl)  ^
  ! +                                                       |               |
  ! +   Dhydtz(isl) etaMid - - - - - - - -   isl-1/2   - -  v  dzmiSV(isl)--+
  ! +                                                                       |
  ! +   eta_SV(isl-1)      ---------------   isl-1   -----                  v

  ! +--Transfert       Coefficients
  ! +  ----------------------------

  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 =     (dz_dSV(isl)  *eta_SV(ikl,isl-1) & ! eta at layers
            +dz_dSV(isl-1)*eta_SV(ikl,isl)  ) & !     interface
            /(2.0* dzmiSV(isl))                      ! LSdzsv implicit !
  ! #GA     etaMid = sqrt(dz_dSV(isl)  *eta_SV(ikl,isl-1)   ! Idem, geometric
  ! #GA.                 *dz_dSV(isl-1)*eta_SV(ikl,isl)  )  !       average
  ! #GA.           /(2.0* dzmiSV(isl))                      ! (Vauclin&al.1979)
      Dhydif          =    s1__SV(ist) & ! Hydraul.Diffusi.
            *(etaMid         **(   bCHdSV(ist)+2.))           ! DR97, Eqn.(3.36)
      Dhydtz(ikl,isl) =    Dhydif*dt__SV & !
            /(dzmiSV(isl) & !
            *LSdzsv(ikl))        !
      Dhydtz(ikl,isl) =    Dhydtz(ikl,isl) * ist__s & ! Soil
            +0.5*dzmiSV(isl)*LSdzsv(ikl)     * ist__w   ! Water bodies

    END DO
  END DO
       isl=-nsol
    DO ikl=1,knonv
      Dhydtz(ikl,isl) =    0.0                        !
    END DO


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

  ! +--Soil/Snow Interior
  ! +  ^^^^^^^^^^^^^^^^^^

  DO   isl=0,-nsol,-1
    DO ikl=1,knonv
     ist             = isotSV(ikl)
     eta_SV(ikl,isl) = max(epsi,           eta_SV(ikl,isl))
    END DO
  END DO

  DO   isl=-nsol,-nsol+1
    DO ikl=1,knonv
      etaaux(ikl,isl) =  eta_SV(ikl,isl)
    END DO
  END DO

  DO   isl=-nsol+1,-1
    DO ikl=1,knonv
      ist      =         isotSV(ikl)
      ikm      = nkhy *  eta_SV(ikl,isl-1) / etadSV(ist)
      ik0      = nkhy *  eta_SV(ikl,isl)   / etadSV(ist)
      ikp      = nkhy *  eta_SV(ikl,isl+1) / etadSV(ist)

      if(ikm<0.or.ik0<0.or.ikp<0)then
       print *,"CRASH1 in sisvat_qso.f on pixel (i,j,n)", &
             ii__SV(ikl),jj__SV(ikl),nn__SV(ikl)
       print *,"decrease your time step or increase ntphys "// &
             "and ntdiff in time_steps.f"
       stop
      endif


      Elem_A   =         Dhydtz(ikl,isl) &
            -  aKdtSV3(ist,ikm)* dziiSV(isl)  *LSdzsv(ikl)
      Elem_B   =      - (Dhydtz(ikl,isl) &
            +Dhydtz(ikl,isl+1) &
            -aKdtSV3(ist,ik0)*(dziiSV(isl+1) &
            -dzi_SV(isl) )*LSdzsv(ikl))
      Elem_C   =         Dhydtz(ikl,isl+1) &
            +  aKdtSV3(ist,ikp)* dzi_SV(isl+1)*LSdzsv(ikl)
      Diag_A(ikl,isl) =  dz_8SV(isl)        *LSdzsv(ikl) &
            -Implic            * Elem_A
      Diag_B(ikl,isl) =  dz34SV(isl)        *LSdzsv(ikl) &
            -Implic            * Elem_B
      Diag_C(ikl,isl) =  dz_8SV(isl)        *LSdzsv(ikl) &
            -Implic            * Elem_C

      Term_D(ikl,isl) = (dz_8SV(isl) *LSdzsv(ikl) &
            +Explic      *Elem_A     )*eta_SV(ikl,isl-1) &
            + (dz34SV(isl) *LSdzsv(ikl) &
            +Explic      *Elem_B     )*eta_SV(ikl,isl) &
            + (dz_8SV(isl) *LSdzsv(ikl) &
            +Explic      *Elem_C     )*eta_SV(ikl,isl+1) &
            + (bKdtSV3(ist,ikp)* dzi_SV(isl+1) &
            +bKdtSV3(ist,ik0)*(dziiSV(isl+1) &
            -dzi_SV(isl)  ) &
            -bKdtSV3(ist,ikm)* dziiSV(isl)   ) &
            * LSdzsv(ikl)
    END DO
  END DO

       isl=-nsol
    DO ikl=1,knonv
      ist      =         isotSV(ikl)
  ! #       FreeDr   =         FreeD0            *  min(ist,1)
      FreeDr   =         iWaFSV(ikl)       *  min(ist,1)
      ik0      = nkhy *  eta_SV(ikl,isl  ) / etadSV(ist)
      ikp      = nkhy *  eta_SV(ikl,isl+1) / etadSV(ist)

      if(ik0<0.or.ikp<0)then
       print *,"CRASH2 in sisvat_qso.f on pixel (i,j,n)", &
             ii__SV(ikl),jj__SV(ikl),nn__SV(ikl)
       print *,"decrease your time step or increase ntphys "// &
             "and ntdiff in time_steps.f"
       stop
      endif

      Elem_A   =         0.
      Elem_B   =      - (Dhydtz(ikl,isl+1) &
            -aKdtSV3(ist,ik0)*(dziiSV(isl+1)*LSdzsv(ikl) &
            -FreeDr                  ))
      Elem_C   =         Dhydtz(ikl,isl+1) &
            +  aKdtSV3(ist,ikp)* dzi_SV(isl+1)*LSdzsv(ikl)
      Diag_A(ikl,isl) =  0.
      Diag_B(ikl,isl) =  dz78SV(isl) *LSdzsv(ikl) &
            -Implic      *Elem_B
      Diag_C(ikl,isl) =  dz_8SV(isl) *LSdzsv(ikl) &
            -Implic      *Elem_C

      Term_D(ikl,isl) = (dz78SV(isl) *LSdzsv(ikl) &
            +Explic      *Elem_B     )*eta_SV(ikl,isl) &
            + (dz_8SV(isl) *LSdzsv(ikl) &
            +Explic      *Elem_C     )*eta_SV(ikl,isl+1) &
            + (bKdtSV3(ist,ikp)* dzi_SV(isl+1)*LSdzsv(ikl) &
            +bKdtSV3(ist,ik0)*(dziiSV(isl+1)*LSdzsv(ikl) &
            -FreeDr                  ))
    END DO

       isl=0
    DO ikl=1,knonv
      ist      =         isotSV(ikl)
      ikm      = nkhy *  eta_SV(ikl,isl-1) / etadSV(ist)
      ik0      = nkhy *  eta_SV(ikl,isl)   / etadSV(ist)
      Elem_A   =         Dhydtz(ikl,isl) &
            -  aKdtSV3(ist,ikm)* dziiSV(isl)*LSdzsv(ikl)
      Elem_B   =      - (Dhydtz(ikl,isl) &
            +aKdtSV3(ist,ik0)* dzi_SV(isl)*LSdzsv(ikl))
      Elem_C   =         0.
      Diag_A(ikl,isl) =  dz_8SV(isl) *LSdzsv(ikl) &
            -  Implic      *Elem_A
      Diag_B(ikl,isl) =  dz78SV(isl) *LSdzsv(ikl) &
            -  Implic      *Elem_B
      Diag_C(ikl,isl) =  0.
  ! +
      Term_D(ikl,isl) = (dz_8SV(isl) *LSdzsv(ikl) &
            +Explic      *Elem_A     )*eta_SV(ikl,isl-1) &
            + (dz78SV(isl) *LSdzsv(ikl) &
            +Explic      *Elem_B     )*eta_SV(ikl,isl) &
            - (bKdtSV3(ist,ik0)* dzi_SV(isl) &
            +bKdtSV3(ist,ikm)* dziiSV(isl))*LSdzsv(ikl) &
            + dt__SV *(HLs_sv(ikl)    *     (1-min(1,isnoSV(ikl))) &
            / (ro_Wat *dz_dSV(0) * Lx_H2O(ikl)) &
  !XF bug 17/05/2017
            +drr_SV(ikl))/ro_Wat
    END DO

    DO ikl=1,knonv
     drr_SV(ikl)=0. ! drr is included in the 1st soil layer
    ENDDO

  ! +
  ! +
  ! +--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   isl=-nsol+1,0
    DO ikl=      1,knonv
      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
      eta_SV(ikl,-nsol) = Term_D(ikl,-nsol)/Aux__P(ikl,-nsol)
    END DO
  ! +
  DO   isl=-nsol+1,0
    DO ikl=      1,knonv
      eta_SV(ikl,isl)   =(Term_D(ikl,isl) &
            -Diag_A(ikl,isl)  *eta_SV(ikl,isl-1)) &
            /Aux__P(ikl,isl)
    END DO
  END DO

  ! +--Backward Sweep
  ! +  ^^^^^^^^^^^^^^
  DO   isl=-1,-nsol,-1
    DO ikl= 1,knonv
      eta_SV(ikl,isl)   = Aux__Q(ikl,isl)  *eta_SV(ikl,isl+1) &
            +eta_SV(ikl,isl)
    END DO
  END DO


  ! +--Horton RunOFF Intensity
  ! +  =======================

  DO   isl=0,-nsol,-1
    DO ikl=1,knonv
      ist    =   isotSV(ikl)                   ! Soil Type
      SatRat =  (eta_SV(ikl,isl)-etadSV(ist)) & ! OverSaturation Rate
            *ro_Wat         *dzAvSV(isl) & !
            *LSdzsv(ikl) & !
            /dt__SV        !
      SoRnOF(ikl)     =          SoRnOF(ikl) & !
            + max(zero,SatRat)       !
      RuofSV(ikl,3)   = RuofSV(ikl,3) + &
            max(zero,SatRat)
      eta_SV(ikl,isl) = max(epsi & !
  ! #ED.                         +etamSV(isotSV(ikl))!
            ,eta_SV(ikl,isl))   !
      eta_SV(ikl,isl) = min(eta_SV(ikl,isl) & !
            ,etadSV(ist)    )   !
    END DO
  END DO

  ! +--IO, for Verification
  ! +  ~~~~~~~~~~~~~~~~~~~~
  ! #WR     write(6,6010)
 6010   format(/,1x)
  DO   isl= 0,-nsol,-1
    DO ikl= 1,knonv
      ist      =          isotSV(ikl)
      ikp      = nkhy  *  eta_SV(ikl,isl)  /etadSV(ist)
      Khydsv(ikl,isl)   =(aKdtSV3(ist,ikp)  *eta_SV(ikl,isl) &
            +bKdtSV3(ist,ikp)) *2.0/dt__SV
  ! #WR     write(6,6011) ikl,isl,eta_SV(ikl,isl)*1.e3,
  ! #WR.                  ikp,    aKdtSV3(ist,ikp),bKdtSV3(ist,ikp),
  ! #WR.                          Khydsv(ikl,isl)
 6011   format(2i3,f8.1,i3,3e12.3)
    END DO
  END DO


  ! +--Additional RunOFF Intensity
  ! +  ===========================

    DO ikl=1,knonv
      ist      =          isotSV(ikl)
      ik0      = nkhy  *  etaaux(ikl,-nsol  ) /etadSV(ist)
  ! #       FreeDr   =          FreeD0            *  min(ist,1)
      FreeDr   =          iWaFSV(ikl)       *  min(ist,1)
      SoRnOF(ikl) =  SoRnOF(ikl) &
            + (aKdtSV3(ist,ik0)*(etaaux(ikl,-nsol)*Explic &
            +eta_SV(ikl,-nsol)*Implic) &
            + bKdtSV3(ist,ik0)                           ) &
            * FreeDr          *ro_Wat           /dt__SV
    RuofSV(ikl,3) = RuofSV(ikl,3) &
          + (aKdtSV3(ist,ik0)*(etaaux(ikl,-nsol)*Explic &
          +eta_SV(ikl,-nsol)*Implic) &
          + bKdtSV3(ist,ik0)                           ) &
          * FreeDr          *ro_Wat           /dt__SV

  ! +--Full Run OFF: Update
  ! +  ~~~~~~~~~~~~~~~~~~~~
      RnofSV(ikl)   = RnofSV(ikl)   + SoRnOF(ikl)
      RuofSV(ikl,4) = RuofSV(ikl,4) + SoRnOF(ikl)
    END DO


  ! +--Temperature Correction due to a changed Soil Energy Content
  ! +  ===========================================================

  ! +!!!    Mettre en oeuvre le couplage humidit?-?nergie


  ! +--Bumps/Asperites Treatment
  ! +  =========================

  ! +--Average over Bump Depth (z0soil)
  ! +  --------------------------------

  ! #BP       z_Bump      = 0.
  ! #BP     DO ikl=1,knonv
  ! #BP       etBump(ikl) = 0.
  ! #BP     END DO
  ! +
  ! #BP DO     isl=0,-nsol,-1
  ! #BP       z0Bump      = z_Bump
  ! #BP       z_Bump      = z_Bump      +  dzAvSV(isl)
  ! #BP   IF (z_Bump.lt.z0soil)                                       THEN
  ! #BP     DO ikl=1,knonv
  ! #BP       etBump(ikl) = etBump(ikl) +  dzAvSV(isl)   *eta_SV(ikl,isl)
  ! #BP     END DO
  ! #BP   END IF
  ! #BP   IF (z_Bump.gt.z0soil.AND.z0Bump.lt.z0soil)                  THEN
  ! #BP     DO ikl=1,knonv
  ! #BP       etBump(ikl) = etBump(ikl) + (z0soil-z0Bump)*eta_SV(ikl,isl)
  ! #BP       etBump(ikl) = etBump(ikl) /  z0soil
  ! #BP     END DO
  ! #BP   END IF
  ! #BP END DO


  ! +--Correction
  ! +  ----------

  ! #BP       z_Bump      = 0.
  ! #BP DO     isl=0,-nsol,-1
  ! #BP       z0Bump =  z_Bump
  ! #BP       z_Bump =  z_Bump +dzAvSV(isl)
  ! #BP   IF (z_Bump.lt.z0soil)                                       THEN
  ! #BP     DO ikl=1,knonv
  ! #BP       eta_SV(ikl,isl) = etBump(ikl)
  ! #BP     END DO
  ! #BP   END IF
  ! #BP   IF (z_Bump.gt.z0soil.AND.z0Bump.lt.z0soil)                  THEN
  ! #BP       dzBump          =    z_Bump -  z0soil
  ! #BP     DO ikl=1,knonv
  ! #BP       eta_SV(ikl,isl) =(etBump(ikl)    *(dzAvSV(isl)-dzBump)
  ! #BP.                      + eta_SV(ikl,isl)*             dzBump)
  ! #BP.                      /                  dzAvSV(isl)
  ! #BP     END DO
  ! #BP   END IF
  ! #BP END DO


  ! +--Positive Definite
  ! +  =================

  ! #BP DO   isl= 0,-nsol,-1
  ! #BP   DO ikl= 1,knonv
  ! #BP     eta_SV(ikl,isl)   =          max(epsi,eta_SV(ikl,isl))
  ! #BP   END DO
  ! #BP END DO


  ! +--Water  Budget (OUT)
  ! +  ===================

  ! #m0   DO ikl=1,knonv
  ! #m0     Wats_d(ikl) = Wats_d(ikl)                    !
  ! #m0.                + drr_SV(ikl)     *zero          ! Precipitation is
  ! +                                      \______________ already included
  ! #m0.                + HLs_sv(ikl)
  ! #m0.          *(1-min(isnoSV(ikl),1)) /Lx_H2O(ikl)   ! Evaporation
  ! #m0.                - SoRnOF(ikl)                    ! Soil RunOFF Contrib.
  ! #m0     Wats_1(ikl) = 0.                             !
  ! #mw     Evapor(ikl) = HLs_sv(ikl)     *dt__SV        !
  ! #mw.          *(1-min(isnoSV(ikl),1)) /Lx_H2O(ikl)   !
  ! #m0   END DO

  ! #m0 DO   isl= -nsol,0
  ! #m0   DO ikl=1,knonv
  ! #m0     Wats_d(ikl) = Wats_d(ikl)                    !
  ! #m0   END DO
  ! #m0 END DO
  ! #m0   DO ikl=1,knonv
  ! #m0     Wats_d(ikl) = Wats_d(ikl)     *dt__SV        !
  ! #m0   END DO

  ! #m0      isl= -nsol
  ! #m0   DO ikl=1,knonv
  ! #m0     Wats_1(ikl) = Wats_1(ikl)
  ! #m0.      + ro_Wat *( eta_SV(ikl,isl)   *dz78SV(isl)
  ! #m0.                + eta_SV(ikl,isl+1) *dz_8SV(isl) ) *LSdzsv(ikl)
  ! #m0   END DO

  ! #m0 DO   isl= -nsol+1,-1
  ! #m0   DO ikl=1,knonv
  ! #m0     Wats_1(ikl) = Wats_1(ikl)
  ! #m0.      + ro_Wat *( eta_SV(ikl,isl)   *dz34SV(isl)
  ! #m0.                +(eta_SV(ikl,isl-1)
  ! #m0.                 +eta_SV(ikl,isl+1))*dz_8SV(isl) ) *LSdzsv(ikl)
  ! #m0   END DO
  ! #m0 END DO

  ! #m0      isl=  0
  ! #m0   DO ikl=1,knonv
  ! #m0     Wats_1(ikl) = Wats_1(ikl)
  ! #m0.      + ro_Wat *( eta_SV(ikl,isl)   *dz78SV(isl)
  ! #m0.                + eta_SV(ikl,isl-1) *dz_8SV(isl) ) *LSdzsv(ikl)
  ! #m0   END DO


  return
end subroutine sisvat_qso
