subroutine SISVAT_BSn
 
C +------------------------------------------------------------------------+
C | MAR          SISVAT_BSn                               04-apr-2020  MAR |
C |   SubRoutine SISVAT_BSn treats Snow Erosion                            |
C |   (not deposition anymore since 2-jun 2018)                            |
C |                                                                        |
C |   SISVAT_bsn computes the snow erosion mass according to both the      |
C |   theoretical maximum erosion amount computed in SISVATesbl and the    |
C |   availability of snow (currently in the uppermost snow layer only)    |
C |                                                                        |
C |   Preprocessing  Option: SISVAT IO (not always a standard preprocess.) |
C |   ^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^                                     |
C |   FILE                 |      CONTENT                                  |
C |   ~~~~~~~~~~~~~~~~~~~~~+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
C | # stdout               | #sb: OUTPUT of Snow Erosion                   |
C |                        |      unit  6, SubRoutine  SISVAT_BSn **ONLY** |
C +------------------------------------------------------------------------+
 
 
 
 
C +--General Variables
C +  =================

      use VARphy
      use VAR_SV
      use VARdSV
      use VARxSV
      use VARySV
 
 
       IMPLICIT NONE

C +--Local Variables
C +  ===============
 
 
      integer  ikl   ,isn
      real     h_mmWE                        ! Eroded Snow Layer Min Thickness
      real     dbsaux(knonv)                 ! Drift Amount   (Dummy Variable)
      real     dzweqo,dzweqn,bsno_x          ! Conversion variables for erosion
      real     dz_new,rho_new
      real     snofOK                        ! Threshd Snow Fall
      real     Fac                           ! Correction factor for erosion
      real     densif                        ! Densification rate if erosion
 
C +--DATA
C +  ====
 
      data     h_mmWE  / 0.01e00  /          ! Eroded Snow Layer Min Thickness
 
C +--EROSION
C +  =======
 
      !DO isn = nsno,2,-1
      DO ikl = 1,knonv
 
        isn         = isnoSV(ikl)
        dzweqo      = dzsnSV(ikl,isn) *ro__SV(ikl,isn)      ! [kg/m2, mm w.e.]
 
        bsno_x      = min(0.,dbs_SV(ikl))
c       Fac         = min(1.,max(1-(ro__SV(ikl,isn)/700.),0.)**2)
c       Fac         = min(1.,max(1-(qsnoSV(ikl)*1000/30.),0.))
c       bsno_x      = bsno_x*Fac
 
        dzweqn      = dzweqo + bsno_x
        dzweqn      = max(dzweqn,h_mmWE)
        dzweqn      = min(dzweqn,dzweqo)
cXF
        dbs_SV(ikl) = dbs_SV(ikl)    +(dzweqo -dzweqn)
        dbs_Er(ikl) = dbs_Er(ikl)    +(dzweqo -dzweqn)
        dzsnSV(ikl,isn) =              dzweqn
     .                       /max(epsi,ro__SV(ikl,isn))
 
        ! Densification of the uppermost snow layer if erosion:
        if((dzweqo-dzweqn)>0                    .and.
     .     dzsnSV(ikl,isn)>0                    .and.
     .     ro__SV(ikl,max(1,isnoSV(ikl)))<roBdSV) then
 
        !characteristic time scale for drifting snow compaction set to 24h
        !linear densification rate [kg/m3/s] over 24h
        densif         = (450. - frsno) / (3600*24)
 
        !Attenuation of compaction rate from 450 to 500 kg/m3
        Fac         = 1-((ro__SV(ikl,max(1,isnoSV(ikl)))
     .                        -roBdSV)/(500.-roBdSV))
        Fac         = max(0.,min(1.,Fac))
 
        if (ro__SV(ikl,max(1,isnoSV(ikl)))>roBdSV) then
          densif=densif*Fac
        endif
 
        rho_new        = min(roBdSV,ro__SV(ikl,isn)+densif*dt__SV)
        dz_new         = dzsnSV(ikl,isn)*ro__SV(ikl,isn)/rho_new
        ro__SV(ikl,isn)=rho_new
        dzsnSV(ikl,isn)=dz_new
        endif
 
        if(dzsnSV(ikl,isn)>0 .and.dzsnSV(ikl,isn)<0.0001)then
        dbs_SV(ikl) = dbs_SV(ikl)+ dzsnSV(ikl,isn)*ro__SV(ikl,isn)
        dbs_Er(ikl) = dbs_Er(ikl)+ dzsnSV(ikl,isn)*ro__SV(ikl,isn)
        dzsnSV(ikl,isn) = 0
        ro__SV(ikl,isn) = 0
        isnoSV(ikl)     = max(0,isnoSV(ikl) - 1)
        endif
 
      END DO
      !END DO
 
      return
      END