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

C |                                                                        |
C | # OPTIONS: #GF: Saturation Front                                       |
C | # ^^^^^^^  #GH: Saturation Front allows Horton Runoff                  |
C | #          #GA: Soil Humidity Geometric Average                        |
C | #          #BP: Parameterization of Terrain Bumps                      |
C |                                                                        |
C |                                                                        |
C +------------------------------------------------------------------------+
 
 
 
 
C +--Global Variables
C +  ================
 
      use VARphy
      use VAR_SV
      use VARdSV
      use VAR0SV
      use VARxSV
      use VARySV

     
      IMPLICIT NONE

 
C +--OUTPUT
C +  ------
 
! 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
 
 
C +--Internal Variables
C +  ==================
 
      integer   isl   ,jsl   ,ist   ,ikl      !
      integer   ikm   ,ikp   ,ik0   ,ik1      !
      integer   ist__s,ist__w                 ! Soil/Water Body Identifier
c #BP real      z0soil                        ! Soil Surface Bumps Height  [m]
c #BP real      z_Bump                        !(Partly)Bumpy Layers Height [m]
c #BP real      z0Bump                        !        Bumpy Layers Height [m]
c #BP real      dzBump                        ! Lowest Bumpy Layer:
 
c #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
! ~~~~~~~~~~~~~~~~~~~
c #mw logical         mwopen                  ! IO   Switch
c #mw common/Sm_qSo_L/mwopen                  !
c #mw real     hourwr,timewr                  !
c #mw common/Sm_qSo_R/timewr                  !
c #mw real            Evapor(knonv)           !
 
 
C +--Internal DATA
C +  =============
 
c #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
 
 
C +--Gravitational Flow
C +  ==================
 
C +...    METHOD: Surface Water Flux saturates successively the soil layers
C +       ^^^^^^  from up to below, but is limited by infiltration capacity.
C +               Hydraulic Conductivity again contributes after this step,
C +               not redundantly because of a constant (saturated) profile.
 
C +--Flux  Limitor
C +  ^^^^^^^^^^^^^
           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
C +
          Khydav = ist__s    * Ks_dSV(ist)         ! DR97  Assumption
     .           + ist__w    * vK_dSV              ! Water Bodies
C +
          Wg_MAX = ro_Wat     *Dhydif              ! MAXimum  Infiltration
     .           *(etadSV(ist)-eta_SV(ikl,isl))    !          Rate
     .           /(dzAvSV(isl)*LSdzsv(ikl)    )    !
     .          +  ro_Wat     *Khydav              !
 
C +--Surface Horton RunOFF
C +  ^^^^^^^^^^^^^^^^^^^^^
          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
 
c #GF DO   isl=0,-nsol,-1
c #GF   DO ikl=1,knonv
c #GF     ist    = isotSV(ikl)                     ! Soil Type
c #GF     ist__s = min(ist, 1)                     ! 1 => Soil
c #GF     ist__w = 1 - ist__s                      ! 1 => Water Body
 
C +--Water Diffusion
C +  ^^^^^^^^^^^^^^^
c #GF     Dhydif = s1__SV(ist)
c #GF.               *max(epsi,eta_SV(ikl,isl))    ! Hydraulic Diffusivity
c #GF.                      **(bCHdSV(ist)+2.)     ! DR97, Eqn.(3.36)
c #GF     Dhydif = ist__s    * Dhydif              !
c #GF.           + ist__w    * vK_dSV              ! Water Bodies
 
C +--Water Conduction (without Horton Runoff)
C +  ^^^^^^^^^^^^^^^^
c #GF     Khyd_f =             Ks_dSV(ist)
C +...    Uses saturated K ==> Horton Runoff ~0    !
 
C +--Water Conduction (with    Horton Runoff)
C +  ^^^^^^^^^^^^^^^^
c #GH     ik0    = nkhy       *eta_SV(ikl,isl)
c #GH.                        /etadSV(ist)
c #GH     eta__f         =            1.
c #GH.   -aKdtSV3(ist,ik0)/(2. *dzAvSV(isl)
c #GH.                        *LSdzsv(ikl))
c #GH     eta__f         = max(eps_21,eta__f)
c #GH     eta__f         = min(etadSV(ist),
c #GH.                         eta_SV(ikl,isl) +
c #GH.   (aKdtSV3(ist,ik0)     *eta_SV(ikl,isl)
c #GH.   +bKdtSV3(ist,ik0))   /(dzAvSV(isl)
c #GH.                        *LSdzsv(ikl))
c #GH.                       / eta__f          )
c #GH     eta__f         = .5*(eta_SV(ikl,isl)
c #GH.                        +eta__f)
 
c #gh     eta__f         =     eta_SV(ikl,isl)
 
c #GH     ik0    = nkhy       *eta__f
c #GH.                        /etadSV(ist)
c #GH     Khyd_f =
c #GH.   (aKdtSV3(ist,ik0)     *eta__f
c #GH.   +bKdtSV3(ist,ik0))    /dt__SV
 
c #GF     Khydav = ist__s    * Khyd_f              ! DR97  Assumption
c #GF.           + ist__w    * vK_dSV              ! Water Bodies
 
C +--Gravitational Flow
C +  ^^^^^^^^^^^^^^^^^^
c #GF     Wg_MAX =                                 ! MAXimum  Infiltration
c #GF.             ro_Wat     *Dhydif              !          Rate
c #GF.           *(etadSV(ist)-eta_SV(ikl,isl))    !
c #GF.           /(dzAvSV(isl)*LSdzsv(ikl)    )    !
c #GF.          +  ro_Wat     *Khydav              !
c #GF   END DO
c #GF END DO
c #GF   DO ikl=1,knonv
c #GF     SoRnOF(ikl)     =    SoRnOF(ikl)         ! RunOFF Intensity
c #GF.                    +    drr_SV(ikl)         ! [kg/m2/s]
C +!!!    Inclure la possibilite de creer une mare sur un bedrock impermeable
c #GF     drr_SV(ikl) = 0.
c #GF   END DO
 
 
C +--Temperature Correction due to a changed Soil Energy Content
C +  ===========================================================
 
C +!!!    Mettre en oeuvre le couplage humidit?-?nergie
 
 
C +--Full Resolution of the Richard's Equation
C +  =========================================
 
C +...    METHOD: Water content evolution results from water fluxes
C +       ^^^^^^  at the layer boundaries
C +               Conductivity is approximated by a piecewise linear profile.
C +               Semi-Implicit Crank-Nicholson scheme is used.
C +              (Bruen, 1997, Sensitivity of hydrological processes
C +                            at the land-atmosphere interface.
C +                            Proc. Royal Irish Academy,  IGBP symposium
C +                            on global change and the Irish Environment.
C +                            Publ.: Maynooth)
 
C +                      - - - - - - - -   isl+1/2   - -  ^
C +                                                       |
C +   eta_SV(isl)        ---------------   isl     -----  +--dz_dSV(isl)  ^
C +                                                       |               |
C +   Dhydtz(isl) etaMid - - - - - - - -   isl-1/2   - -  v  dzmiSV(isl)--+
C +                                                                       |
C +   eta_SV(isl-1)      ---------------   isl-1   -----                  v
 
C +--Transfert       Coefficients
C +  ----------------------------
 
      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 !
c #GA     etaMid = sqrt(dz_dSV(isl)  *eta_SV(ikl,isl-1)   ! Idem, geometric
c #GA.                 *dz_dSV(isl-1)*eta_SV(ikl,isl)  )  !       average
c #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
 
 
C +--Tridiagonal Elimination: Set Up
C +  -------------------------------
 
C +--Soil/Snow Interior
C +  ^^^^^^^^^^^^^^^^^^
 
      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)
c #       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.
C +
          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))
cXF 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
 
C +
C +
C +--Tridiagonal Elimination
C +  =======================
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
C +
      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
C +
        DO ikl=      1,knonv
          eta_SV(ikl,-nsol) = Term_D(ikl,-nsol)/Aux__P(ikl,-nsol)
        END DO
C +
      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
 
C +--Backward Sweep
C +  ^^^^^^^^^^^^^^
      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
 
 
C +--Horton RunOFF Intensity
C +  =======================
 
      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               !
c #ED.                         +etamSV(isotSV(ikl))!
     .                         ,eta_SV(ikl,isl))   !
          eta_SV(ikl,isl) = min(eta_SV(ikl,isl)    !
     .                         ,etadSV(ist)    )   !
        END DO
      END DO
 
C +--IO, for Verification
C +  ~~~~~~~~~~~~~~~~~~~~
c #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
c #WR     write(6,6011) ikl,isl,eta_SV(ikl,isl)*1.e3,
c #WR.                  ikp,    aKdtSV3(ist,ikp),bKdtSV3(ist,ikp),
c #WR.                          Khydsv(ikl,isl)
 6011     format(2i3,f8.1,i3,3e12.3)
        END DO
      END DO
 
 
C +--Additional RunOFF Intensity
C +  ===========================
 
        DO ikl=1,knonv
          ist      =          isotSV(ikl)
          ik0      = nkhy  *  etaaux(ikl,-nsol  ) /etadSV(ist)
c #       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
 
C +--Full Run OFF: Update
C +  ~~~~~~~~~~~~~~~~~~~~
          RnofSV(ikl)   = RnofSV(ikl)   + SoRnOF(ikl)
          RuofSV(ikl,4) = RuofSV(ikl,4) + SoRnOF(ikl)
        END DO
 
 
C +--Temperature Correction due to a changed Soil Energy Content
C +  ===========================================================
 
C +!!!    Mettre en oeuvre le couplage humidit?-?nergie
 
 
C +--Bumps/Asperites Treatment
C +  =========================
 
C +--Average over Bump Depth (z0soil)
C +  --------------------------------
 
c #BP       z_Bump      = 0.
c #BP     DO ikl=1,knonv
c #BP       etBump(ikl) = 0.
c #BP     END DO
C +
c #BP DO     isl=0,-nsol,-1
c #BP       z0Bump      = z_Bump
c #BP       z_Bump      = z_Bump      +  dzAvSV(isl)
c #BP   IF (z_Bump.lt.z0soil)                                       THEN
c #BP     DO ikl=1,knonv
c #BP       etBump(ikl) = etBump(ikl) +  dzAvSV(isl)   *eta_SV(ikl,isl)
c #BP     END DO
c #BP   END IF
c #BP   IF (z_Bump.gt.z0soil.AND.z0Bump.lt.z0soil)                  THEN
c #BP     DO ikl=1,knonv
c #BP       etBump(ikl) = etBump(ikl) + (z0soil-z0Bump)*eta_SV(ikl,isl)
c #BP       etBump(ikl) = etBump(ikl) /  z0soil
c #BP     END DO
c #BP   END IF
c #BP END DO
 
 
C +--Correction
C +  ----------
 
c #BP       z_Bump      = 0.
c #BP DO     isl=0,-nsol,-1
c #BP       z0Bump =  z_Bump
c #BP       z_Bump =  z_Bump +dzAvSV(isl)
c #BP   IF (z_Bump.lt.z0soil)                                       THEN
c #BP     DO ikl=1,knonv
c #BP       eta_SV(ikl,isl) = etBump(ikl)
c #BP     END DO
c #BP   END IF
c #BP   IF (z_Bump.gt.z0soil.AND.z0Bump.lt.z0soil)                  THEN
c #BP       dzBump          =    z_Bump -  z0soil
c #BP     DO ikl=1,knonv
c #BP       eta_SV(ikl,isl) =(etBump(ikl)    *(dzAvSV(isl)-dzBump)
c #BP.                      + eta_SV(ikl,isl)*             dzBump)
c #BP.                      /                  dzAvSV(isl)
c #BP     END DO
c #BP   END IF
c #BP END DO
 
 
C +--Positive Definite
C +  =================
 
c #BP DO   isl= 0,-nsol,-1
c #BP   DO ikl= 1,knonv
c #BP     eta_SV(ikl,isl)   =          max(epsi,eta_SV(ikl,isl))
c #BP   END DO
c #BP END DO
 
 
C +--Water  Budget (OUT)
C +  ===================
 
! #m0   DO ikl=1,knonv
! #m0     Wats_d(ikl) = Wats_d(ikl)                    !
! #m0.                + drr_SV(ikl)     *zero          ! Precipitation is
C +                                      \______________ 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.                             !
c #mw     Evapor(ikl) = HLs_sv(ikl)     *dt__SV        !
c #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
