subroutine SISVAT_BSn

  ! +------------------------------------------------------------------------+
  ! | MAR          SISVAT_BSn                               04-apr-2020  MAR |
  ! |   SubRoutine SISVAT_BSn treats Snow Erosion                            |
  ! |   (not deposition anymore since 2-jun 2018)                            |
  ! |                                                                        |
  ! |   SISVAT_bsn computes the snow erosion mass according to both the      |
  ! |   theoretical maximum erosion amount computed in inlandsis and the     |
  ! |   availability of snow (currently in the uppermost snow layer only)    |
  ! |                                                                        |
  ! +------------------------------------------------------------------------+




  ! +--General Variables
  ! +  =================

  use VARphy
  use VAR_SV
  use VARdSV
  use VARxSV
  use VARySV


   IMPLICIT NONE

  ! +--Local Variables
  ! +  ===============


  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

  ! +--DATA
  ! +  ====

  data     h_mmWE  / 0.01e00  /          ! Eroded Snow Layer Min Thickness

  ! +--EROSION
  ! +  =======

  ! !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))
    ! Fac         = min(1.,max(1-(ro__SV(ikl,isn)/700.),0.)**2)
    ! Fac         = min(1.,max(1-(qsnoSV(ikl)*1000/30.),0.))
    ! bsno_x      = bsno_x*Fac

    dzweqn      = dzweqo + bsno_x
    dzweqn      = max(dzweqn,h_mmWE)
    dzweqn      = min(dzweqn,dzweqo)
  !XF
    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 SUBROUTINE SISVAT_BSn