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

      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
 
 
C +--Internal Variables
C +  ==================
 
      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
 
C +--OUTPUT of SISVAT Trace Statistics (see assignation in PHY_SISVAT)
C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      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
 
C +--Energy and Mass Budget
C +  ~~~~~~~~~~~~~~~~~~~~~~
! #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                 !
 
C +--Slush Diagnostic: IO
C +  ~~~~~~~~~~~~~~~~~~~~
! #vu logical         su_opn                 ! IO   Switch
! #vu common/SI_qSn_L/su_opn                 !
 
 
C +--DATA
C +  ====
 
      data      dzepsi/0.0001/                ! Minim. Snow Layer Thickness (!)
c #?? data      dz_Min/0.005/                 ! Minim. Snow Layer Thickness
c ... 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)
 
C +--Energy Budget (IN)
C +  ==================
 
! #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
 
 
C +--Water  Budget (IN)
C +  ==================
 
! #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
 
 
C +--Snow Melt Budget
C +  ================
 
! #m1   DO ikl=1,knonv
! #m1     SImelt(ikl) = 0.
! #m1     SIrnof(ikl) = rusnSV(ikl) + RnofSV(ikl) * dt__SV
! #m1   END DO
 
 
C +--Initialization
C +  ==============
 
      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
 
 
C +--Melting/Freezing Energy
C +  =======================
 
C +...REMARK: Snow liquid Water Temperature assumed = TfSnow
C +   ^^^^^^
        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
 
 
C +--Surficial Water Status
C +  ----------------------
 
        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
C +--Energy, store Previous Content
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          dTSnow      = TsisSV(ikl,isn) -          TfSnow
          EExcsv(ikl) = EExcsv(ikl)
     .                + ro__SV(ikl,isn) * Cn_dSV * dTSnow
     .                                           * dzsnSV(ikl,isn)
          TsisSV(ikl,isn) =                        TfSnow
 
C +--Water,  store Previous Content
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          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.
 
 
C +--Melting  if EExcsv > 0
C +  ======================
 
          EnMelt      =    max(zero,          EExcsv(ikl) )
 
C +--Energy Consumption
C +  ^^^^^^^^^^^^^^^^^^
          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)
 
C +--Water  Production
C +  ^^^^^^^^^^^^^^^^^
          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))
 
C +--Snow History
C +  ^^^^^^^^^^^^
          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)      )           !
 
 
C +--Freezing if EExcsv < 0
C +  ======================
 
          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)
 
 

C +--Snow Water Content
C +  ==================
 
C +--Percolation Velocity
C +  ^^^^^^^^^^^^^^^^^^^^
c #PW     SGDiam    = 1.6d-4
c #PW.              + 1.1d-13 *(ro__SV(ikl,isn)*ro__SV(ikl,isn)
c #PW.                         *ro__SV(ikl,isn)*ro__SV(ikl,isn))
 
C +--Pore   Volume [-]
C +  ^^^^^^^^^^^^^^^^^
          rosDry      =(1.     - eta_SV(ikl,isn))* ro__SV(ikl,isn) !
          PorVol      = 1.     - rosDry          / ro_Ice          !
          PorVol      =      max(PorVol          , zero  )         !
 
C +--Water  Retention
C +  ^^^^^^^^^^^^^^^^
          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))
 
C +--Pore Hole Close OFF
C +  ^^^^^^^^^^^^^^^^^^^
          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
 
cXF
          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
 
c         if (isn>1.and.isn<nsno     .and.
c    .      ro__SV(ikl,isn-1)>900    .and.
c    .      ro__SV(ikl,isn)  >roCdSV .and.
c    .      ro__SV(ikl,isn)  <900    .and.
c    .      TsisSV(ikl,isn)  >273.14 .and.
c    .      TsisSV(ikl,isn+1)<273.15 .and.
c    .      drr_SV(ikl)      >0)     then
c          TsisSV(ikl,isn)=273.14
c          PClose = 1
c         endif
 
cXF
          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
 
 
C +--Remove Zero-Thickness Layers
C +  ============================
 
 1000 CONTINUE
           isnitr =          0
      DO   ikl=1,knonv
           isnUpD =          0
           isinew =          0
cXF


        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
 
 
C +--New upper Limit of the non erodible Snow (istoSV .GT. 1)
C +  ========================================
 
      DO   ikl=1,knonv
           nh =     0
cXF
        DO isn=  isnoSV(ikl),1,-1
           nh =    nh + isn* min(istoSV(ikl,isn)-1,1)*max(0,1-nh)
        ENDDO
           zc =     0.
           zt =     0.
cXF
        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
 
 
C +--Energy Budget (OUT)
C +  ===================
 
! #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
 
 
C +--"Negative Heat" from supercooled rain
C +   ------------------------------------
 
      DO ikl=1,knonv
          EExcsv(ikl) = EExcsv(ikl) + EExdum(ikl)
 
 
C +--Surficial Water Run OFF
C +  -----------------------
 
          rusnew      = rusnSV(ikl) * SWf_SV(ikl)
 
          if(isnoSV(ikl)<=1 .OR. opt_runoff_ac) rusnew = 0.
          !if(ivgtSV(ikl)>=1) rusnew = 0.
 
c #EU                        rusnew = 0.
c #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
 
 
C +--Percolation down the Continental Ice Pack
C +  -----------------------------------------
 
        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
 
cXF 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
 
C +--Slush Formation (Activated. CAUTION: ADD RunOff Possibility before Activation)
C +  ---------------  ^^^^^^^  ^^^
 
      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
 
C +--Available Additional Pore   Volume [-]
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          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
 
C +--Impact of the Sublimation/Deposition on the Surface Mass Balance
C +  ================================================================
 
      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)
 
 
C +--Additional Energy
C +  -----------------
 
c #VH     AdEnrg = dzVap0 * ro__SV(ikl,isnoSV(ikl))           ! Water   Vapor
c #VH.            *C__Wat *(TsisSV(ikl,isnoSV(ikl)) -TfSnow)  ! Sensible Heat
 
c #aH     B_Enrg =(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
c #aH.            -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
c #aH.           /(1.          + dzVap0 /max(epsi,dzsnSV(ikl,isn)))
c #aH     eta_SV(ikl,isn) =
c #aH.           max(zero,unun +(B_Enrg
c #aH.                         -(TsisSV(ikl,isn) -TfSnow)*Cn_dSV)
c #aH.                          /Lf_H2O                          )
c #aH     TsisSV(ikl,isn) =    ( B_Enrg
c #aH.                         +(1.              -eta_SV(ikl,isn))
c #aH.                          *Lf_H2O                          )
c #aH.                         / Cn_dSV
c #aH.                         + TfSnow
 
! #e1     STOP "PLEASE add Energy (#aH) from deposition/sublimation"
 
 
C +--Update of the upper Snow layer Thickness
C +  ----------------------------------------
 
          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
 
 
C +--Energy Budget (OUT)
C +  ===================
 
! #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
 
 
C +--Snow/I Budget
C +  -------------
 
! #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
 
 
C +--Anticipated Disappearance of a rapidly Melting too thin Last Snow Layer
C +  =======================================================================
 
      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
 
C +--Water  Production
C +  ^^^^^^^^^^^^^^^^^
        drr_SV(ikl)   = drr_SV(ikl)
     .                + ro__SV(ikl,1) * z_Melt /dt__SV
      END DO
 
 
C +--Update Nb of Layers
C +  ===================
 
      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
 
 
C +--Water  Budget (OUT)
C +  ===================
 
! #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