

subroutine SISVAT_qSn &
        ( &
  ! #e1.                      EqSn_0,EqSn_1,EqSn_d
  ! #m1.                     ,SIsubl,SImelt,SIrnof
        )

  ! +------------------------------------------------------------------------+
  ! | MAR          SISVAT_qSn                           Fri 29-Jul-2011  MAR |
  ! |   SubRoutine SISVAT_qSn updates  the Snow Water Content                |
  ! +------------------------------------------------------------------------+
  ! |                                                                        |
  ! |   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                      |
  ! |   ^^^^^                                                                |
  ! |                                                                        |
  ! |   INPUT:   TaT_SV   : SBL Top    Temperature                       [K] |
  ! |   ^^^^^    dt__SV   : Time Step                                    [s] |
  ! |                                                                        |
  ! |   INPUT /  drr_SV   : Rain Intensity                         [kg/m2/s] |
  ! |   OUTPUT:  dzsnSV   : Snow Layer Thickness                         [m] |
  ! |   ^^^^^^   eta_SV   : Snow Water Content                       [m3/m3] |
  ! |            ro__SV   : Snow/Soil  Volumic Mass                  [kg/m3] |
  ! |            TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
  ! |                     & Snow     Temperatures (layers  1,2,...,nsno) [K] |
  ! |                                                                        |
  ! |   OUTPUT:  SWS_SV   : Surficial Water Status                           |
  ! |   ^^^^^^                                                               |
  ! |            EExcsv   : Snow Energy in Excess, initial Forcing    [J/m2] |
  ! |            EqSn_d   : Snow Energy in Excess, remaining          [J/m2] |
  ! |            EqSn_0   : Snow Energy, before Phase Change          [J/m2] |
  ! |            EqSn_1   : Snow Energy, after  Phase Change          [J/m2] |
  ! |            SIsubl   : Snow sublimed/deposed Mass             [mm w.e.] |
  ! |            SImelt   : Snow Melted           Mass             [mm w.e.] |
  ! |            SIrnof   : Surficial Water + Run OFF Change       [mm w.e.] |
  ! |                                                                        |
  ! |   Internal Variables:                                                  |
  ! |   ^^^^^^^^^^^^^^^^^^                                                   |
  ! |                                                                        |
  ! | # OPTIONS: #E0: IO for Verification: Energy       Budget               |
  ! | # ^^^^^^^                                                              |
  ! | #          #su: IO for Verification: Slush        Diagnostic           |
  ! |                                                                        |
  ! |                                                                        |
  ! +------------------------------------------------------------------------+




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

  use VARphy
  use VAR_SV
  use VARdSV
  use VAR0SV
  use VARxSV
  use VARySV
  use surface_data, only: is_ok_slush,opt_runoff_ac


  IMPLICIT NONE


  ! Energy          Budget
  ! ~~~~~~~~~~~~~~~~~~~~~~
  ! #e1 real     EqSn_d(knonv)                 ! Energy in Excess, initial
  ! #e1 real     EqSn_0(knonv)                 ! Snow Energy, befor Phase Change
  ! #vm real     EqSn01(knonv)                 ! Snow Energy, after Phase Change
  ! #vm real     EqSn02(knonv)                 ! Snow Energy, after Phase Change
                                         ! !              .AND. Last Melting
  ! #e1 real     EqSn_1(knonv)                 ! Snow Energy, after Phase Change
                                         ! !              .AND. Mass Redistr.
  ! Snow/Ice (Mass) Budget
  ! ~~~~~~~~~~~~~~~~~~~~~~
  ! #m1 real     SIsubl(knonv)                 ! Snow Deposed Mass
  ! #m1 real     SImelt(knonv)                 ! Snow Melted  Mass
  ! #m1 real     SIrnof(knonv)                 ! Local Surficial Water + Run OFF


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

  integer :: ikl   ,isn                    !
  integer :: nh                            ! Non erodible Snow: up.lay.Index
  integer :: LayrOK                        ! 1 (0)  if In(Above) Snow Pack
  integer :: k_face                        ! 1 (0)  if Crystal(no) faceted
  integer :: LastOK                        ! 1 ==>  1! Snow Layer
  integer :: NOLayr                        ! 1     Layer  Update
  integer :: noSnow(knonv)                 ! Nb of Layers Updater
  integer :: kSlush                        ! Slush Switch
  real :: dTSnow                        ! Temperature                  [C]
  real :: EExdum(knonv)                 ! Energy in Excess when no Snow
  real :: OKmelt                        ! 1 (0)  if        (no) Melting
  real :: EnMelt                        ! Energy in excess, for Melting
  real :: SnHLat                        ! Energy consumed   in  Melting
  real :: AdEnrg,B_Enrg                 ! Additional Energy from  Vapor
  real :: dzVap0,dzVap1                 ! Vaporized Thickness          [m]
  real :: dzMelt(knonv)                 ! Melted    Thickness          [m]
  real :: rosDry                        ! Snow volumic Mass if no Water in
  real :: PorVol                        ! Pore volume
  real :: PClose                        ! Pore Hole Close OFF Switch
  real :: SGDiam                        !      Snow Grain Diameter
  real :: SGDmax                        ! Max. Snow Grain Diameter
  real :: rWater                        ! Retained Water           [kg/m2]
  real :: drrNEW                        ! New available Water      [kg/m2]
  real :: rdzNEW                        ! Snow          Mass       [kg/m2]
  real :: rdzsno                        ! Snow          Mass       [kg/m2]
  real :: EnFrez                        ! Energy Release    in  Freezing
  real :: WaFrez                        ! Water  consumed   in  Melting
  real :: RapdOK                        ! 1. ==> Snow melts rapidly
  real :: ThinOK                        ! 1. ==> Snow Layer is thin
  real :: dzepsi                        ! Minim. Snow Layer Thickness (!)
  real :: dz_Min                        ! Minim. Snow Layer Thickness
  real :: z_Melt                        ! Last (thin) Layer Melting
  real :: rusnew                        ! Surficial Water Thickness   [mm]
  real :: zWater                        ! Max Slush Water Thickness   [mm]
  real :: zSlush                        !     Slush Water Thickness   [mm]
  real :: ro_new                        ! New Snow/ice  Density    [kg/m3]
  real :: zc,zt                         ! Non erod.Snow Thickness[mm w.e.]
  real :: rusnSV0(knonv)
  real :: Tsave

  ! +--OUTPUT of SISVAT Trace Statistics (see assignation in PHY_SISVAT)
  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  integer :: isnnew,isinew,isnUpD,isnitr

  ! OUTPUT in SISVAT at specified i,j,k,n (see assignation in PHY_SISVAT)
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! #wx integer             iSV_v1,jSV_v1,nSV_v1,kSV_v1,lSV_v1
  ! #wx common/SISVAT_EV/   iSV_v1,jSV_v1,nSV_v1,kSV_v1,lSV_v1

  ! +--Energy and Mass Budget
  ! +  ~~~~~~~~~~~~~~~~~~~~~~
  ! #vm real     WqSn_0(knonv)                 ! Snow Water+Forcing  Initial
  ! #vm real     WqSn_1(knonv)                 ! Snow Water+Forcing, Final
  ! #vm logical         emopen                 ! IO   Switch
  ! #vm common/Se_qSn_L/emopen                 !
  ! #vm integer         no_err                 !
  ! #vm common/Se_qSn_I/no_err                 !
  ! #vm real     hourer,timeer                 !
  ! #vm common/Se_qSn_R/timeer                 !

  ! +--Slush Diagnostic: IO
  ! +  ~~~~~~~~~~~~~~~~~~~~
  ! #vu logical         su_opn                 ! IO   Switch
  ! #vu common/SI_qSn_L/su_opn                 !


  ! +--DATA
  ! +  ====

  data      dzepsi/0.0001/                ! Minim. Snow Layer Thickness (!)
  ! #?? data      dz_Min/0.005/                 ! Minim. Snow Layer Thickness
  ! ... Warning: Too high for Col de Porte: precludes 1st snow (layer) apparition
  data      dz_Min/2.5e-3/                ! Minim. Snow Layer Thickness
  data      SGDmax/0.003/                 ! Maxim. Snow Grain Diameter  [m]
                                          ! ! (Rowe et al. 1995, JGR p.16268)

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

  ! #e1   DO ikl=1,knonv
  ! #e1     EqSn_0(ikl) = 0.
  ! #e1   END DO
  ! #e1 DO   isn=nsno,1,-1
  ! #e1   DO ikl=1,knonv
  ! #e1     EqSn_0(ikl) = EqSn_0(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
  ! #e1.                *(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
  ! #e1.                 -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
  ! #e1   END DO
  ! #e1 END DO


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

  ! #vm   DO ikl=1,knonv
  ! #vm     WqSn_0(ikl) = drr_SV(ikl) * dt__SV
  ! #vm.                 +rusnSV(ikl)
  ! #vm   END DO
  ! #vm DO   isn=nsno,1,-1
  ! #vm   DO ikl=1,knonv
  ! #vm     WqSn_0(ikl) = WqSn_0(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
  ! #vm   END DO
  ! #vm END DO


  ! +--Snow Melt Budget
  ! +  ================

  ! #m1   DO ikl=1,knonv
  ! #m1     SImelt(ikl) = 0.
  ! #m1     SIrnof(ikl) = rusnSV(ikl) + RnofSV(ikl) * dt__SV
  ! #m1   END DO


  ! +--Initialization
  ! +  ==============

  DO ikl=1,knonv
    noSnow(ikl)   = 0                   ! Nb of Layers Updater
    ispiSV(ikl)   = 0                   ! Pore Hole Close OFF Index
                                        ! ! (assumed to be the Top of
                                        ! !  the surimposed Ice Layer)
    zn5_SV(ikl)   = 0.
    rusnSV0(ikl)  = 0.

  END DO


  ! +--Melting/Freezing Energy
  ! +  =======================

  ! +...REMARK: Snow liquid Water Temperature assumed = TfSnow
  ! +   ^^^^^^
    DO ikl=1,knonv
      EExdum(ikl) = drr_SV(ikl)     * C__Wat *(TaT_SV(ikl)-TfSnow) &
            * dt__SV
      EExcsv(ikl) = EExdum(ikl)     *    min(1,isnoSV(ikl)) ! Snow exists
      EExdum(ikl) = EExdum(ikl)     -          EExcsv(ikl)  !
  ! #e1     EqSn_d(ikl) = EExcsv(ikl)                             !
    END DO


  ! +--Surficial Water Status
  ! +  ----------------------

    DO ikl=1,knonv
      SWS_SV(ikl) = max(zero,sign(unun,TfSnow &
            -TsisSV(ikl,isnoSV(ikl))))
    END DO

  DO ikl=1,knonv

  DO isn=min(nsno,isnoSV(ikl)+1),1,-1
  ! EV          DO isn=nsno,1,-1
  ! +--Energy, store Previous Content
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      dTSnow      = TsisSV(ikl,isn) -          TfSnow
      EExcsv(ikl) = EExcsv(ikl) &
            + ro__SV(ikl,isn) * Cn_dSV * dTSnow &
            * dzsnSV(ikl,isn)
      TsisSV(ikl,isn) =                        TfSnow

  ! +--Water,  store Previous Content
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      drr_SV(ikl) = drr_SV(ikl) &
            + ro__SV(ikl,isn)          * eta_SV(ikl,isn) &
            * dzsnSV(ikl,isn) &
            / dt__SV
      ro__SV(ikl,isn) = &
            ro__SV(ikl,isn) *(1.     - eta_SV(ikl,isn))
      eta_SV(ikl,isn) =  0.


  ! +--Melting  if EExcsv > 0
  ! +  ======================

      EnMelt      =    max(zero,          EExcsv(ikl) )

  ! +--Energy Consumption
  ! +  ^^^^^^^^^^^^^^^^^^
      SnHLat      = ro__SV(ikl,isn) * Lf_H2O
      dzMelt(ikl) = EnMelt      / max(SnHLat,    epsi )
      noSnow(ikl) = noSnow(ikl) &
            + max(zero  ,sign(unun,dzMelt(ikl) & !
            -dzsnSV(ikl ,isn))) & ! 1 if full Melt
            *min(1     , max(0 ,1+isnoSV(ikl)-isn))          ! 1 in the  Pack
      dzMelt(ikl) = &
            min(dzsnSV(ikl, isn),dzMelt(ikl))
      dzsnSV(ikl,isn) = &
            dzsnSV(ikl,isn) -dzMelt(ikl)
      zn5_SV(ikl) = zn5_SV(ikl)     +dzMelt(ikl)
      EExcsv(ikl) = EExcsv(ikl)     -dzMelt(ikl)*SnHLat
      wem_SV(ikl) = wem_SV(ikl)     -dzMelt(ikl)*ro__SV(ikl,isn)

  ! +--Water  Production
  ! +  ^^^^^^^^^^^^^^^^^
      drr_SV(ikl) = drr_SV(ikl) &
            + ro__SV(ikl,isn) * dzMelt(ikl)/dt__SV
  ! #m1     SImelt(ikl) = SImelt(ikl)
  ! #m1.                + ro__SV(ikl,isn) * dzMelt(ikl)
      OKmelt      =max(zero,sign(unun,drr_SV(ikl)-epsi))

  ! +--Snow History
  ! +  ^^^^^^^^^^^^
      k_face          =       min(    istoSV(ikl,isn),istdSV(1)) & ! = 1  if
            *max(0,2-istoSV(ikl,isn)          ) ! faceted
      istoSV(ikl,isn) = & !
            (1.-OKmelt) *               istoSV(ikl,isn) & !
            +     OKmelt  *((1-k_face) *  istdSV(2) & !
            +   k_face  *  istdSV(3)      )           !


  ! +--Freezing if EExcsv < 0
  ! +  ======================

      rdzsno      =          ro__SV(ikl,isn) * dzsnSV(ikl ,isn)
      LayrOK      = min(   1, max(0          , isnoSV(ikl)-isn+1))
      EnFrez      = min(zero,                  EExcsv(ikl))
      WaFrez      =   -(     EnFrez          * LayrOK / Lf_H2O)
      drrNEW      = max(zero,drr_SV(ikl)     - WaFrez / dt__SV)
      WaFrez      =    (     drr_SV(ikl)     - drrNEW)* dt__SV
      drr_SV(ikl) =          drrNEW
      EExcsv(ikl) =          EExcsv(ikl)     + WaFrez * Lf_H2O
      EnFrez      = min(zero,EExcsv(ikl))    * LayrOK
      rdzNEW      = WaFrez + rdzsno
      ro__SV(ikl,isn) =      rdzNEW /max(epsi, dzsnSV(ikl,isn))
      TsisSV(ikl,isn) =      TfSnow &
            + EnFrez /(Cn_dSV *max(epsi, rdzNEW)        )
      EExcsv(ikl) =          EExcsv(ikl)     - EnFrez
      wer_SV(ikl) = WaFrez &
            + wer_SV(ikl)



  ! +--Snow Water Content
  ! +  ==================

  ! +--Percolation Velocity
  ! +  ^^^^^^^^^^^^^^^^^^^^
  ! #PW     SGDiam    = 1.6d-4
  ! #PW.              + 1.1d-13 *(ro__SV(ikl,isn)*ro__SV(ikl,isn)
  ! #PW.                         *ro__SV(ikl,isn)*ro__SV(ikl,isn))

  ! +--Pore   Volume [-]
  ! +  ^^^^^^^^^^^^^^^^^
      rosDry      =(1.     - eta_SV(ikl,isn))* ro__SV(ikl,isn) !
      PorVol      = 1.     - rosDry          / ro_Ice          !
      PorVol      =      max(PorVol          , zero  )         !

  ! +--Water  Retention
  ! +  ^^^^^^^^^^^^^^^^
      rWater      = ws0dSV * PorVol * ro_Wat * dzsnSV(ikl,isn)
      drrNEW      = max(zero,drr_SV(ikl)     - rWater /dt__SV)
      rWater      =    (     drr_SV(ikl)     - drrNEW)*dt__SV
      drr_SV(ikl)     =      drrNEW
      rdzNEW          =      rWater &
            + rosDry          * dzsnSV(ikl,isn)
      eta_SV(ikl,isn) =      rWater / max(epsi,rdzNEW)
      ro__SV(ikl,isn) =      rdzNEW / max(epsi,dzsnSV(ikl,isn))

  ! +--Pore Hole Close OFF
  ! +  ^^^^^^^^^^^^^^^^^^^
      PClose = max(zero, &
            sign(unun,ro__SV(ikl,isn) &
            -roCdSV         ))
      ispiSV(ikl) =          ispiSV(ikl)      *(1.-PClose) &
            +      max(ispiSV(ikl),isn)    * Pclose
      PClose = max(0   , & ! Water under SuPer.Ice
            min (1   ,ispiSV(ikl) & ! contributes to
            -isn            ))   ! Surficial   Water

  !XF
      if(ro__SV(ikl,isn) >= roCdSV.and.ro__SV(ikl,1)<900) &
            PClose = min(0.50,PClose * &
            (1.-(ro_ice-ro__SV(ikl,isn))/(ro_ice-roCdSV)))

      PClose = max(0.,min(1.,PClose))

      if(isn==1) then
           PClose = 1
       ispiSV(ikl)= max(ispiSV(ikl),1)
      endif

      if(drr_SV(ikl)    >0  .and.TsisSV(ikl,isn)>273.14) then
       if((ro__SV(ikl,isn)>900.and.ro__SV(ikl,isn)<920).or. &
             ro__SV(ikl,isn)>950) then
         dzsnSV(ikl,isn) = dzsnSV(ikl,isn)*ro__SV(ikl,isn)/ro_ice
         ro__SV(ikl,isn) = ro_ice
         PClose          = 1
       endif
      endif

      ! if (isn>1.and.isn<nsno     .and.
  !    .      ro__SV(ikl,isn-1)>900    .and.
  !    .      ro__SV(ikl,isn)  >roCdSV .and.
  !    .      ro__SV(ikl,isn)  <900    .and.
  !    .      TsisSV(ikl,isn)  >273.14 .and.
  !    .      TsisSV(ikl,isn+1)<273.15 .and.
  !    .      drr_SV(ikl)      >0)     then
      !  TsisSV(ikl,isn)=273.14
      !  PClose = 1
      ! endif

  !XF
      rusnSV(ikl) =          rusnSV(ikl) &
            +          drr_SV(ikl) *dt__SV * PClose
      rusnSV0(ikl)=         rusnSV0(ikl) &
            +          drr_SV(ikl) *dt__SV * PClose
      drr_SV(ikl) =          drr_SV(ikl)      *(1.-PClose)

    END DO

  END DO


  ! +--Remove Zero-Thickness Layers
  ! +  ============================

 1000   CONTINUE
       isnitr =          0
  DO   ikl=1,knonv
       isnUpD =          0
       isinew =          0
  !XF


    DO isn=1,min(nsno-1,isnoSV(ikl))
       isnnew =(unun-max(zero  ,sign(unun,dzsnSV(ikl,isn)-dzepsi))) &
             *     max(0     , min(1   ,isnoSV(ikl) +1 -isn ))
       isnUpD =      max(isnUpD,          isnnew)
       isnitr =      max(isnitr,          isnnew)
       isinew =      isn*isnUpD *max(0, 1-isinew) & ! LowerMost  0-Layer
             +isinew       ! Index
       dzsnSV(ikl,isn) =                  dzsnSV(ikl,isn+isnnew)
       ro__SV(ikl,isn) =                  ro__SV(ikl,isn+isnnew)
       TsisSV(ikl,isn) =                  TsisSV(ikl,isn+isnnew)
       eta_SV(ikl,isn) =                  eta_SV(ikl,isn+isnnew)
       G1snSV(ikl,isn) =                  G1snSV(ikl,isn+isnnew)
       G2snSV(ikl,isn) =                  G2snSV(ikl,isn+isnnew)
       dzsnSV(ikl,isn+isnnew) =(1-isnnew)*dzsnSV(ikl,isn+isnnew)
       ro__SV(ikl,isn+isnnew) =(1-isnnew)*ro__SV(ikl,isn+isnnew)
       eta_SV(ikl,isn+isnnew) =(1-isnnew)*eta_SV(ikl,isn+isnnew)
       G1snSV(ikl,isn+isnnew) =(1-isnnew)*G1snSV(ikl,isn+isnnew)
       G2snSV(ikl,isn+isnnew) =(1-isnnew)*G2snSV(ikl,isn+isnnew)

    END DO
       isnoSV(ikl)   =   isnoSV(ikl)-isnUpD            ! Nb of Snow   Layer
       ispiSV(ikl)   =   ispiSV(ikl) & ! Nb of SuperI Layer
             -isnUpD *max(0,min(ispiSV(ikl)-isinew,1))        ! Update  if I=0

  END DO

  IF  (isnitr.GT.0)                                       GO TO 1000


  ! +--New upper Limit of the non erodible Snow (istoSV .GT. 1)
  ! +  ========================================

  DO   ikl=1,knonv
       nh =     0
  !XF
    DO isn=  isnoSV(ikl),1,-1
       nh =    nh + isn* min(istoSV(ikl,isn)-1,1)*max(0,1-nh)
    ENDDO
       zc =     0.
       zt =     0.
  !XF
    DO isn=1,isnoSV(ikl)
       zc =    zc +          dzsnSV(ikl,isn) *ro__SV(ikl,isn) &
             * max(0,min(1,nh+1-isn))
       zt =    zt +          dzsnSV(ikl,isn) *ro__SV(ikl,isn)
    END DO
       zWE_SV(ikl) =                 zt
       zWEcSV(ikl) = min(zWEcSV(ikl),zt)
       zWEcSV(ikl) = max(zWEcSV(ikl),zc)
  END DO


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

  ! #vm   DO ikl=1,knonv
  ! #vm     EqSn01(ikl) =-EqSn_0(ikl)
  ! #vm.                 -EExcsv(ikl)
  ! #vm   END DO
  ! #vm DO   isn=nsno,1,-1
  ! #vm   DO ikl=1,knonv
  ! #vm     EqSn01(ikl) = EqSn01(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
  ! #vm.                *(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
  ! #vm.                 -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
  ! #vm   END DO
  ! #vm END DO


  ! +--"Negative Heat" from supercooled rain
  ! +   ------------------------------------

  DO ikl=1,knonv
      EExcsv(ikl) = EExcsv(ikl) + EExdum(ikl)


  ! +--Surficial Water Run OFF
  ! +  -----------------------

      rusnew      = rusnSV(ikl) * SWf_SV(ikl)

      if(isnoSV(ikl)<=1 .OR. opt_runoff_ac) rusnew = 0.
      ! !if(ivgtSV(ikl)>=1) rusnew = 0.

  ! #EU                        rusnew = 0.
  ! #AC               rusnew = 0.

      RnofSV(ikl) = RnofSV(ikl) &
            +(rusnSV(ikl) - rusnew     ) / dt__SV
      RuofSV(ikl,1) = RuofSV(ikl,1) &
            +(rusnSV(ikl) - rusnew     ) / dt__SV
      RuofSV(ikl,4) = RuofSV(ikl,4) &
            +(rusnSV0(ikl)             ) / dt__SV
      rusnSV(ikl) = rusnew
  END DO


  ! +--Percolation down the Continental Ice Pack
  ! +  -----------------------------------------

    DO ikl=1,knonv
      drr_SV(ikl) = drr_SV(ikl) + rusnSV(ikl) &
            * (1-min(1,ispiSV(ikl)))/ dt__SV
      rusnSV(ikl) = rusnSV(ikl) &
            *    min(1,ispiSV(ikl))
    END DO

  !XF removal of too thin snowlayers if TT> 275.15 + bug if TT>> 273.15
    DO ikl=1,knonv
     zt=0.
     DO isn=1,isnoSV(ikl)
      zt=zt+dzsnSV(ikl,isn)
     ENDDO

     if(zt<0.005+(TaT_SV(ikl)-TfSnow)/1000..and. &
           isnoSV(ikl)             >0         .and. &
           TaT_SV(ikl)             >=TfSnow   .and. &
           istoSV(ikl,isnoSV(ikl)) >1       ) then
      DO isn=1,isnoSV(ikl)
       drr_SV(ikl)    = drr_SV(ikl) &
             + dzsnSV(ikl,isn)*ro__SV(ikl,isn) /dt__SV
       dzsnSV(ikl,isn)= 0.

      ENDDO
      isnoSV(ikl)     = 0
     endif
    ENDDO

  ! +--Slush Formation (Activated. CAUTION: ADD RunOff Possibility before Activation)
  ! +  ---------------  ^^^^^^^  ^^^

  IF (is_ok_slush) THEN

  DO  ikl=1,knonv
   DO isn=1,isnoSV(ikl)
      kSlush = min(1,max(0,isn+1-ispiSV(ikl)))        ! Slush Switch

  ! +--Available Additional Pore   Volume [-]
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      PorVol = 1. - ro__SV(ikl,isn) & ! [--]
            *(1. - eta_SV(ikl,isn))/ ro_Ice & !
            -      eta_SV(ikl,isn) & !
            *ro__SV(ikl,isn) / ro_Wat           !
      PorVol =  max(PorVol          , zero  )          !
      zWater =      dzsnSV(ikl,isn) * PorVol * 1000. & ! [mm] OR [kg/m2]
            * (1. -SWS_SV(ikl) & ! 0 <=> freezing
            *(1 -min(1,iabs(isn-isnoSV(ikl)))))  ! 1 <=> isn=isnoSV
      zSlush =  min(rusnSV(ikl)     , zWater)          ! [mm] OR [kg/m2]
      ro_new      =(dzsnSV(ikl,isn) * ro__SV(ikl,isn) & !
            +zSlush                           ) & !
            / max(dzsnSV(ikl,isn) , epsi           ) !
      if(ro_new<ro_Ice+20) then ! MAX 940kg/m3         !
       rusnSV(ikl)  = rusnSV(ikl)          - zSlush    ! [mm] OR [kg/m2]
       RuofSV(ikl,4)= max(0.,RuofSV(ikl,4) - zSlush/dt__SV)
       eta_SV(ikl,isn) =(ro_new - ro__SV(ikl,isn) & !
             *(1.     - eta_SV(ikl,isn))) & !
             / max (ro_new , epsi            )    !
       ro__SV(ikl,isn) =      ro_new                   !
      endif
    END DO
  END DO
  END IF

  ! +--Impact of the Sublimation/Deposition on the Surface Mass Balance
  ! +  ================================================================

  DO ikl=1,knonv
     isn                     = isnoSV(ikl)
      dzVap0                  =                   dt__SV &
            * HLs_sv(ikl)         * min(isn             , 1   ) &
            /(Lx_H2O(ikl)         * max(ro__SV(ikl,isn) , epsi))
      NOLayr=min(zero,sign(unun,dzsnSV(ikl,isn) + dzVap0))
      dzVap1=min(zero,          dzsnSV(ikl,isn) + dzVap0)


  ! +--Additional Energy
  ! +  -----------------

  ! #VH     AdEnrg = dzVap0 * ro__SV(ikl,isnoSV(ikl))           ! Water   Vapor
  ! #VH.            *C__Wat *(TsisSV(ikl,isnoSV(ikl)) -TfSnow)  ! Sensible Heat

  ! #aH     B_Enrg =(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
  ! #aH.            -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
  ! #aH.           /(1.          + dzVap0 /max(epsi,dzsnSV(ikl,isn)))
  ! #aH     eta_SV(ikl,isn) =
  ! #aH.           max(zero,unun +(B_Enrg
  ! #aH.                         -(TsisSV(ikl,isn) -TfSnow)*Cn_dSV)
  ! #aH.                          /Lf_H2O                          )
  ! #aH     TsisSV(ikl,isn) =    ( B_Enrg
  ! #aH.                         +(1.              -eta_SV(ikl,isn))
  ! #aH.                          *Lf_H2O                          )
  ! #aH.                         / Cn_dSV
  ! #aH.                         + TfSnow

  ! #e1     STOP "PLEASE add Energy (#aH) from deposition/sublimation"


  ! +--Update of the upper Snow layer Thickness
  ! +  ----------------------------------------

      dzsnSV(ikl,isn) = &
            max(zero,  dzsnSV(ikl,isnoSV(ikl)) + dzVap0)
      isnoSV(ikl)     = isnoSV(ikl)             + NOLayr
      isn             = isnoSV(ikl)
      dzsnSV(ikl,isn) = dzsnSV(ikl,isn) + dzVap1
      wes_SV(ikl)     = ro__SV(ikl,isn) * dzVap0

  END DO


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

  ! #vm   DO ikl=1,knonv
  ! #vm     EqSn02(ikl) =-EqSn_0(ikl)
  ! #vm.                 -EExcsv(ikl)
  ! #vm   END DO
  ! #vm DO   isn=nsno,1,-1
  ! #vm   DO ikl=1,knonv
  ! #vm     EqSn02(ikl) = EqSn01(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
  ! #vm.                *(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
  ! #vm.                 -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
  ! #vm   END DO
  ! #vm END DO


  ! +--Snow/I Budget
  ! +  -------------

  ! #m1   DO ikl=1,knonv
  ! #m1     SIsubl(ikl) = dt__SV*HLs_sv(ikl)*min(isnoSV(ikl),1)
  ! #m1.                        /Lx_H2O(ikl)
  ! #m1     SIrnof(ikl) = rusnSV(ikl) + RnofSV(ikl) * dt__SV
  ! #m1.                - SIrnof(ikl)
  ! #m1   END DO


  ! +--Anticipated Disappearance of a rapidly Melting too thin Last Snow Layer
  ! +  =======================================================================

  DO ikl=1,knonv
    LastOK = min(1   , max(0   ,iiceSV(ikl)-isnoSV(ikl)+2) &
          *min(1   ,isnoSV(ikl)-iiceSV(ikl)) &
          +min(1   ,isnoSV(ikl))              )
    RapdOK = max(zero,sign(unun,dzMelt(ikl)-epsi         ))
    ThinOK = max(zero,sign(unun,dz_Min     -dzsnSV(ikl,1)))
    z_Melt = LastOK     *RapdOK*ThinOK
    noSnow(ikl)   = noSnow(ikl)   + z_Melt
    z_Melt        =                 z_Melt *dzsnSV(ikl,1)
    dzsnSV(ikl,1) = dzsnSV(ikl,1) - z_Melt
    EExcsv(ikl)   = EExcsv(ikl)   - z_Melt *ro__SV(ikl,1) &
          *(1.     -eta_SV(ikl,1))*Lf_H2O

  ! +--Water  Production
  ! +  ^^^^^^^^^^^^^^^^^
    drr_SV(ikl)   = drr_SV(ikl) &
          + ro__SV(ikl,1) * z_Melt /dt__SV
  END DO


  ! +--Update Nb of Layers
  ! +  ===================

  DO ikl=1,knonv
    isnoSV(ikl)   = isnoSV(ikl) &
          * min(1,iabs(isnoSV(ikl)-noSnow(ikl)))
  END DO


  ! Energy Budget (OUT)
  ! ===================

  ! #e1   DO ikl=1,knonv
  ! #e1     EqSn_1(ikl) = 0.
  ! #e1   END DO
  ! #e1 DO   isn=nsno,1,-1
  ! #e1   DO ikl=1,knonv
  ! #e1     EqSn_1(ikl) = EqSn_1(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
  ! #e1.                *(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
  ! #e1.                 -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
  ! #e1   END DO
  ! #e1 END DO


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

  ! #vm   DO ikl=1,knonv
  ! #vm     WqSn_0(ikl) = WqSn_0(ikl)
  ! #vm.                + HLs_sv(ikl)    * dt__SV
  ! #vm.             *min(isnoSV(ikl),1) / Lx_H2O(ikl)
  ! #vm     WqSn_1(ikl) = drr_SV(ikl)    * dt__SV
  ! #vm.                + rusnSV(ikl)
  ! #vm.                + RnofSV(ikl)    * dt__SV
  ! #vm   END DO
  ! #vm DO   isn=nsno,1,-1
  ! #vm   DO ikl=1,knonv
  ! #vm     WqSn_1(ikl) = WqSn_1(ikl)
  ! #vm.                + ro__SV(ikl,isn)* dzsnSV(ikl,isn)
  ! #vm   END DO
  ! #vm END DO


  return
end subroutine sisvat_qsn
