

subroutine SISVAT_zSn

  ! +------------------------------------------------------------------------+
  ! | MAR          SISVAT_zSn                                12-07-2019  MAR |
  ! |   SubRoutine SISVAT_zSn manages the Snow Pack vertical Discretization  |
  ! |                                                                        |
  ! +------------------------------------------------------------------------+
  ! |                                                                        |
  ! |   PARAMETERS:  knonv: Total Number of columns =                        |
  ! |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
  ! |                     X       Number of Mosaic Cell per grid box         |
  ! |                                                                        |
  ! |   INPUT /  NLaysv   = New             Snow Layer  Switch               |
  ! |   OUTPUT:  isnoSV   = total Nb of Ice/Snow Layers                      |
  ! |   ^^^^^^   ispiSV   = 0,...,nsno: Uppermost Superimposed Ice Layer     |
  ! |            iiceSV   = total Nb of Ice      Layers                      |
  ! |            istoSV   = 0,...,5 :   Snow     History (see istdSV data)   |
  ! |                                                                        |
  ! |   INPUT /  TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
  ! |   OUTPUT:           & Snow     Temperatures (layers  1,2,...,nsno) [K] |
  ! |   ^^^^^^   ro__SV   : Soil/Snow Volumic Mass                   [kg/m3] |
  ! |            eta_SV   : Soil/Snow Water   Content                [m3/m3] |
  ! |            dzsnSV   : Snow Layer        Thickness                  [m] |
  ! |            G1snSV   : Dendricity (<0) or Sphericity (>0) of Snow Layer |
  ! |            G2snSV   : Sphericity (>0) or Size            of Snow Layer |
  ! |            agsnSV   : Snow       Age                             [day] |
  ! |                                                                        |
  ! |   METHOD:  1) Agregate the thinest Snow Layer                          |
  ! |   ^^^^^^      if a new Snow Layer has been precipitated   (NLaysv = 1) |
  ! |            2) Divide   a too thick Snow Layer except                   |
  ! |               if the maximum Number of Layer is reached                |
  ! |               in this case forces                          NLay_s = 1  |
  ! |            3) Agregate the thinest Snow Layer                          |
  ! |               in order to divide a too thick Snow Layer                |
  ! |               at next Time Step when                       NLay_s = 1  |
  ! |                                                                        |
  ! +------------------------------------------------------------------------+




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


  use VARphy
  use VAR_SV
  use VARdSV
  use VAR0SV
  use VARxSV
  use VARySV
  use surface_data, only: ok_zsn_ii

  IMPLICIT NONE


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

  integer :: ikl   ,isn   ,i               !






  integer :: NLay_s(knonv)                 ! Split Snow Layer         Switch
  integer :: isagr1(knonv)                 ! 1st     Layer History
  integer :: isagr2(knonv)                 ! 2nd     Layer History
  integer :: LstLay                        ! 0 ====> isnoSV = 1
  integer :: isno_n                        ! Snow Normal.Profile
  integer :: iice_n                        ! Ice  Normal.Profile
  integer :: iiceOK                        ! Ice         Switch
  integer :: icemix                        ! 0 ====> Agregated Snow+Ice=Snow
  ! +                                           ! 1                          Ice
  integer :: isn1  (knonv)                 ! 1st layer to stagger
  real :: staggr                        !              stagger  Switch

  real :: WEagre(knonv)                 ! Snow Water Equivalent Thickness
  real :: dzthin(knonv)                 ! Thickness of the thinest layer
  real :: OKthin                        ! Swich ON  a  new thinest layer
  real :: dz_dif                        ! difference from ideal discret.
  real :: thickL                        ! Thick Layer          Indicator
  real :: OK_ICE                        ! Swich ON   uppermost Ice Layer

  real :: Agrege(knonv)                 ! 1. when Agregation constrained
  real :: dzepsi                        ! Min Single Snw Layer Thickness
  real :: dzxmin                        ! Min Acceptable Layer Thickness
  real :: dz_min                        ! Min            Layer Thickness
  real :: dz_max                        ! Max            Layer Thickness
  real :: dzagr1(knonv)                 ! 1st     Layer Thickness
  real :: dzagr2(knonv)                 ! 2nd     Layer Thickness
  real :: T_agr1(knonv)                 ! 1st     Layer Temperature
  real :: T_agr2(knonv)                 ! 2nd     Layer Temperature
  real :: roagr1(knonv)                 ! 1st     Layer Density
  real :: roagr2(knonv)                 ! 2nd     Layer Density
  real :: etagr1(knonv)                 ! 1st     Layer Water Content
  real :: etagr2(knonv)                 ! 2nd     Layer Water Content
  real :: G1agr1(knonv)                 ! 1st     Layer Dendricity/Spher.
  real :: G1agr2(knonv)                 ! 2nd     Layer Dendricity/Spher.
  real :: G2agr1(knonv)                 ! 1st     Layer Sphericity/Size
  real :: G2agr2(knonv)                 ! 2nd     Layer Sphericity/Size
  real :: agagr1(knonv)                 ! 1st     Layer Age
  real :: agagr2(knonv)                 ! 2nd     Layer Age


  ! +--DATA
  ! +  ====

  data      icemix /   0    /             ! 0 ====> Agregated Snow+Ice=Snow
  data      dzepsi /   0.0020/            ! Min single Layer Thickness
  data      dzxmin /   0.0025/            ! Min accept.Layer Thickness
  ! #EU data      dz_min /   0.0050/            ! Min Local  Layer Thickness < SMn
  data      dz_min /   0.0040/            ! Min Local  Layer Thickness < SMn
  data      dz_max /   0.0300/            ! Min Gener. Layer Thickness
  ! +   CAUTION:  dz_max > dz_min*2 is required ! Otherwise re-agregation is
  ! +                                           ! activated  after splitting





  ! +--Constrains Agregation         of too thin  Layers
  ! +  =================================================

  ! +--Search the thinest  non-zero Layer
  ! +  ----------------------------------

    DO ikl=1,knonv
      if(isnoSV(ikl)<=2)             dz_min=max(0.0050,dz_min)

                                      dzepsi=0.0015
      if(ro__SV(ikl,isnoSV(ikl))>920) dzepsi=0.0020

      dzthin(ikl) = 0.                              ! Arbitrary unrealistic
    END DO                                          !       Layer Thickness
  !XF
  DO ikl=1,knonv
  DO   isn=1,isnoSV(ikl)-3 ! no agregation of 3 first snowlayers
                           ! ! XF 04/07/2019

      isno_n    =             isnoSV(ikl)-isn+1     ! Snow Normal.Profile
      iice_n    =             iiceSV(ikl)-isn       ! Ice  Normal.Profile
      iiceOK    = min(1,max(0,iice_n         +1))   ! Ice         Switch
  ! #vz     dz_ref(isn) =                                 !
  ! #vz.          dz_min *((1-iiceOK)*isno_n*isno_n       ! Theoretical Profile
  ! #vz.                 +    iiceOK *    2**iice_n)      !
  ! #vz.               /max(1,isnoSV(ikl))                !
      dz_dif      = max(zero, & ! Actual      Profile
            dz_min & !
            *((1-iiceOK)*isno_n*isno_n & ! Theoretical Profile
            +    iiceOK *2.   **iice_n) & !
            - dzsnSV(ikl, isn)                    )   ! Actual      Profile
  ! #vz     dzwdif(isn) =     dz_dif                      !
      OKthin      = max(zero, & !
            sign(unun, & !
            dz_dif-dzthin(ikl))) & ! 1.=> New thinest Lay.
            * max(0, & ! 1 => .le. isnoSV
            min(1, & ! 1 => isn is in the
            isnoSV(ikl)-isn +1 )) & !          Snow Pack
            * min(unun, & !
  !
  !                   1st additional Condition to accept OKthin
            max(zero, & ! combination
            sign(unun,G1snSV(ikl,      isn  ) & ! G1 with same
            *G1snSV(ikl,max(1,isn-1)))) & !  sign => OK
  !
  !                   2nd additional Condition to accept OKthin
            + max(zero, & ! G1>0
            sign(unun,G1snSV(ikl,      isn   ))) & !  =>OK
  !
  !                   3rd additional Condition to accept OKthin
            + max(zero, & ! dz too small
            sign(unun,dzxmin & !  =>OK
            -dzsnSV(ikl,      isn   ))))!

      i_thin(ikl) =    (1. - OKthin)  * i_thin(ikl) & ! Update   thinest Lay.
            + OKthin   * isn         !                Index
      dzthin(ikl) =    (1. - OKthin)  * dzthin(ikl) & !
            + OKthin   * dz_dif      !
    END DO
  END DO



  DO ikl=1,knonv
  DO   isn=1,isnoSV(ikl)
      OKthin =      max(zero, & !
            sign(unun, & !
            dz_min & !
            -dzsnSV(ikl,isn))) & !
            * max(zero, & ! ON if dz > 0
            sign(unun, & !
            dzsnSV(ikl,isn)-epsi)) & !
            *min(1,max(0, & ! Multiple Snow    Lay.
            min (1, & ! Switch = 1
            isnoSV(ikl) & !   if isno > iice + 1
            -iiceSV(ikl)-1)) & !
  ! +                                                     !
            +int(max(zero, & !
            sign(unun, & !
            dzepsi & ! Minimum accepted for
            -dzsnSV(ikl,isn)))) & ! 1 Snow Layer over Ice
            *int(max(zero, & ! ON if dz > 0
            sign(unun, & !
            dzsnSV(ikl,isn)-epsi))) & !
            *(1 -min (abs(isnoSV(ikl) & ! Switch = 1
            -iiceSV(ikl)-1),1)) & !   if isno = iice + 1
  ! +                                                     !
            +max(0, & ! Ice
            min (1, & ! Switch
            iiceSV(ikl)+1-isn))) & !
            *min(unun, & !
            max(zero, & ! combination
            sign(unun,G1snSV(ikl,      isn  ) & ! G1>0 + G1<0
            *G1snSV(ikl,max(1,isn-1)))) & ! NO
            + max(zero, & !
            sign(unun,G1snSV(ikl,      isn   ))) & !
            + max(zero, & !
            sign(unun,dzxmin & !
            -dzsnSV(ikl,      isn   ))))!
      i_thin(ikl) =    (1. - OKthin)  * i_thin(ikl) & ! Update   thinest Lay.
            + OKthin   * isn         !                Index
    END DO
  END DO




  ! +   ***************
  call SISVAT_zCr
  ! +   ***************


  ! +--Assign the 2 Layers to agregate
  ! +  -------------------------------

    DO ikl=1,knonv
      isn         =    i_thin(ikl)
      if(LIndsv(ikl)>0) isn=min(nsno-1,isn) ! cXF
      isagr1(ikl) =    istoSV(ikl,isn)
      isagr2(ikl) =    istoSV(ikl,isn+LIndsv(ikl))
      dzagr1(ikl) =    dzsnSV(ikl,isn)
      dzagr2(ikl) =    dzsnSV(ikl,isn+LIndsv(ikl))
      T_agr1(ikl) =    TsisSV(ikl,isn)
      T_agr2(ikl) =    TsisSV(ikl,isn+LIndsv(ikl))
      roagr1(ikl) =    ro__SV(ikl,isn)
      roagr2(ikl) =    ro__SV(ikl,isn+LIndsv(ikl))
      etagr1(ikl) =    eta_SV(ikl,isn)
      etagr2(ikl) =    eta_SV(ikl,isn+LIndsv(ikl))
      G1agr1(ikl) =    G1snSV(ikl,isn)
      G1agr2(ikl) =    G1snSV(ikl,isn+LIndsv(ikl))
      G2agr1(ikl) =    G2snSV(ikl,isn)
      G2agr2(ikl) =    G2snSV(ikl,isn+LIndsv(ikl))
      agagr1(ikl) =    agsnSV(ikl,isn)
      agagr2(ikl) =    agsnSV(ikl,isn+LIndsv(ikl))
      LstLay      = min(1,max(  0,isnoSV(ikl) -1))  ! 0  if single Layer
      isnoSV(ikl) =               isnoSV(ikl) & ! decrement   isnoSV
            -(1-LstLay)* max(zero, & ! if downmost  Layer
            sign(unun,eps_21 & ! <  1.e-21 m
            -dzsnSV(ikl,1)))   !
      isnoSV(ikl) = max(   0,     isnoSV(ikl)   )   !
      Agrege(ikl) = max(zero, & !
            sign(unun,dz_min & ! No Agregation
            -dzagr1(ikl)  )) & ! if too thick Layer
            *LstLay & ! if  a single Layer
            * min( max(0   ,isnoSV(ikl)+1 & ! if Agregation
            -i_thin(ikl) & !    with    a Layer
            -LIndsv(ikl)  ),1) !    above the Pack

      WEagre(ikl) = 0.
    END DO


    DO ikl=1,knonv
    DO   isn=1,isnoSV(ikl)
      WEagre(ikl) = WEagre(ikl) + ro__SV(ikl,isn)*dzsnSV(ikl,isn) &
            *min(1,max(0,i_thin(ikl)+1-isn))
    ENDDO
    ENDDO


  ! +--Agregates
  ! +  ---------

  ! +     ***************
    call SISVAT_zAg &
          (isagr1,isagr2,WEagre &
          ,dzagr1,dzagr2,T_agr1,T_agr2 &
          ,roagr1,roagr2,etagr1,etagr2 &
          ,G1agr1,G1agr2,G2agr1,G2agr2 &
          ,agagr1,agagr2,Agrege &
          )
  ! +     ***************


  ! +--Rearranges the Layers
  ! +  ---------------------

  ! +--New (agregated) Snow layer
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
    DO ikl=1,knonv
      isn     =             i_thin(ikl)
      isn     = min(isn,isn+LIndsv(ikl))
      isnoSV(ikl) =  max(0.,isnoSV(ikl) -Agrege(ikl))
      iiceSV(ikl) =         iiceSV(ikl) &
            -max(0,sign(1,iiceSV(ikl) -isn +icemix)) &
            *Agrege(ikl) &
            *max(0,sign(1,iiceSV(ikl) -1          ))
      istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn) &
            +   Agrege(ikl) *isagr1(ikl)
      dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn) &
            +   Agrege(ikl) *dzagr1(ikl)
      TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn) &
            +   Agrege(ikl) *T_agr1(ikl)
      ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn) &
            +   Agrege(ikl) *roagr1(ikl)
      eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn) &
            +   Agrege(ikl) *etagr1(ikl)
      G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn) &
            +   Agrege(ikl) *G1agr1(ikl)
      G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn) &
            +   Agrege(ikl) *G2agr1(ikl)
      agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn) &
            +   Agrege(ikl) *agagr1(ikl)
    END DO

  ! +--Above
  ! +  ^^^^^
    DO ikl=1,knonv
      isn1(ikl)=max(i_thin(ikl),i_thin(ikl)+LIndsv(ikl))
    END DO
    DO i=  1,nsno-1
    DO ikl=1,knonv
        staggr        =  min(1,max(0,i +1 -isn1(ikl)   ))
        istoSV(ikl,i) = (1.-staggr     )*istoSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*istoSV(ikl,i  ) &
              +   Agrege(ikl) *istoSV(ikl,i+1))
        dzsnSV(ikl,i) = (1.-staggr     )*dzsnSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*dzsnSV(ikl,i  ) &
              +   Agrege(ikl) *dzsnSV(ikl,i+1))
        TsisSV(ikl,i) = (1.-staggr     )*TsisSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*TsisSV(ikl,i  ) &
              +   Agrege(ikl) *TsisSV(ikl,i+1))
        ro__SV(ikl,i) = (1.-staggr     )*ro__SV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*ro__SV(ikl,i  ) &
              +   Agrege(ikl) *ro__SV(ikl,i+1))
        eta_SV(ikl,i) = (1.-staggr     )*eta_SV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*eta_SV(ikl,i  ) &
              +   Agrege(ikl) *eta_SV(ikl,i+1))
        G1snSV(ikl,i) = (1.-staggr     )*G1snSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*G1snSV(ikl,i  ) &
              +   Agrege(ikl) *G1snSV(ikl,i+1))
        G2snSV(ikl,i) = (1.-staggr     )*G2snSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*G2snSV(ikl,i  ) &
              +   Agrege(ikl) *G2snSV(ikl,i+1))
        agsnSV(ikl,i) = (1.-staggr     )*agsnSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*agsnSV(ikl,i  ) &
              +   Agrege(ikl) *agsnSV(ikl,i+1))
    END DO
    END DO

    DO ikl=1,knonv
      isn             = min(isnoSV(ikl) +1,nsno)
      istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn)
      dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn)
      TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn)
      ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn)
      eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn)
      G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn)
      G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn)
      agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn)
    END DO




  ! +--Constrains Splitting          of too thick Layers
  ! +  =================================================


  ! +--Search the thickest non-zero Layer
  ! +  ----------------------------------

    DO ikl=1,knonv
      dzthin(ikl) =   0.                            ! Arbitrary unrealistic
    END DO
  DO ikl=1,knonv
    DO   isn=1,isnoSV(ikl)
      isno_n    =             isnoSV(ikl)-isn+1     ! Snow Normal.Profile
      iice_n    =             iiceSV(ikl)-isn       ! Ice  Normal.Profile
      iiceOK    = min(1,max(0,iice_n         +1))   ! Ice         Switch
      dz_dif    =(      dzsnSV(ikl,isn) & ! Actual      Profile
            - dz_max *((1-iiceOK)*isno_n*isno_n & ! Theoretical Profile
            +    iiceOK *2.   **iice_n)  ) & !
            /max(dzsnSV(ikl,isn),epsi)       !
      OKthin      = max(zero, & !
            sign(unun, & !
            dz_dif-dzthin(ikl))) & ! 1.=>New thickest Lay.
            * max(0, & ! 1 =>.le. isnoSV
            min(1, & !
            isnoSV(ikl)-isn +1 ))       !
      i_thin(ikl) =    (1. - OKthin)  * i_thin(ikl) & !  Update thickest Lay.
            + OKthin   * isn         !                Index
      dzthin(ikl) =    (1. - OKthin)  * dzthin(ikl) & !
            + OKthin   * dz_dif      !
    END DO

    isn=max(1,isnoSV(ikl)-3)
    if(dzsnSV(ikl,isn)>0.30) then   ! surface layer > 30cm
     i_thin(ikl) = isn              ! XF 04/07/2019
     dzthin(ikl) = dzsnSV(ikl,isn)
    endif

    isn=max(1,isnoSV(ikl)-2)
    if(dzsnSV(ikl,isn)>0.10) then   ! surface layer > 10cm
     i_thin(ikl) = isn              ! XF 04/07/2019
     dzthin(ikl) = dzsnSV(ikl,isn)
    endif

    isn=max(1,isnoSV(ikl)-1)
    if(dzsnSV(ikl,isn)>0.05) then   ! surface layer > 5cm
     i_thin(ikl) = isn              ! XF 04/07/2019
     dzthin(ikl) = dzsnSV(ikl,isn)
    endif

    isn=max(1,isnoSV(ikl))
    if(dzsnSV(ikl,isn)>0.02) then   ! surface layer > 2cm
     i_thin(ikl) = isn              ! XF 04/07/2019
     dzthin(ikl) = dzsnSV(ikl,isn)
    endif

  END DO

  DO   ikl=1,knonv
      ThickL      = max(zero, & ! 1. => a too   thick
            sign(unun,dzthin(ikl) & !         Layer exists
            -epsi       )) & !
            * max(0,1-max(0   ,     isnoSV(ikl) & ! No spliting allowed
            -nsno+1     ))     ! if isno > nsno - 1
      Agrege(ikl) =               ThickL & ! 1. => effective split
            * max(0,1-max(0   ,     NLaysv(ikl) & !
            +isnoSV(ikl) & !
            -nsno+1     ))     !
      NLay_s(ikl) =               ThickL & ! Agregation
            * max(0,1-max(0   ,     NLaysv(ikl) & ! to allow  Splitting
            +isnoSV(ikl) & !   at next Time Step
            -nsno       )) & !
            -Agrege(ikl)       !
      NLay_s(ikl) = max(0   ,     NLay_s(ikl))      ! Agregation effective
  END DO


  ! +--Rearranges the Layers
  ! +  ---------------------

  DO isn=nsno,2,-1
  DO ikl=1,knonv
    IF (Agrege(ikl).gt.0..AND.i_thin(ikl).lt.isnoSV(ikl))       THEN
      staggr          =  min(1,max(0,isn-i_thin(ikl)    -1)) &
            *  min(1,max(0,    isnoSV(ikl)-isn+2))
      istoSV(ikl,isn) =        staggr  * istoSV(ikl ,isn-1) &
            + (1. -  staggr) * istoSV(ikl ,isn  )
      dzsnSV(ikl,isn) =        staggr  * dzsnSV(ikl ,isn-1) &
            + (1. -  staggr) * dzsnSV(ikl ,isn  )
      TsisSV(ikl,isn) =        staggr  * TsisSV(ikl ,isn-1) &
            + (1. -  staggr) * TsisSV(ikl ,isn  )
      ro__SV(ikl,isn) =        staggr  * ro__SV(ikl ,isn-1) &
            + (1. -  staggr) * ro__SV(ikl ,isn  )
      eta_SV(ikl,isn) =        staggr  * eta_SV(ikl ,isn-1) &
            + (1. -  staggr) * eta_SV(ikl ,isn  )
      G1snSV(ikl,isn) =        staggr  * G1snSV(ikl ,isn-1) &
            + (1. -  staggr) * G1snSV(ikl ,isn  )
      G2snSV(ikl,isn) =        staggr  * G2snSV(ikl ,isn-1) &
            + (1. -  staggr) * G2snSV(ikl ,isn  )
      agsnSV(ikl,isn) =        staggr  * agsnSV(ikl ,isn-1) &
            + (1. -  staggr) * agsnSV(ikl ,isn  )
    END IF
  END DO
  END DO

  DO  ikl=1,knonv
      isn             =     i_thin(ikl)
      dzsnSV(ikl,isn) = 0.5*Agrege(ikl) *dzsnSV(ikl,isn) &
            + (1.-Agrege(ikl))*dzsnSV(ikl,isn)

      isn             = min(i_thin(ikl) +1,nsno)
      istoSV(ikl,isn) =     Agrege(ikl) *istoSV(ikl,isn-1) &
            + (1.-Agrege(ikl))*istoSV(ikl,isn)
      dzsnSV(ikl,isn) =     Agrege(ikl) *dzsnSV(ikl,isn-1) &
            + (1.-Agrege(ikl))*dzsnSV(ikl,isn)
      TsisSV(ikl,isn) =     Agrege(ikl) *TsisSV(ikl,isn-1) &
            + (1.-Agrege(ikl))*TsisSV(ikl,isn)
      ro__SV(ikl,isn) =     Agrege(ikl) *ro__SV(ikl,isn-1) &
            + (1.-Agrege(ikl))*ro__SV(ikl,isn)
      eta_SV(ikl,isn) =     Agrege(ikl) *eta_SV(ikl,isn-1) &
            + (1.-Agrege(ikl))*eta_SV(ikl,isn)
      G1snSV(ikl,isn) =     Agrege(ikl) *G1snSV(ikl,isn-1) &
            + (1.-Agrege(ikl))*G1snSV(ikl,isn)
      G2snSV(ikl,isn) =     Agrege(ikl) *G2snSV(ikl,isn-1) &
            + (1.-Agrege(ikl))*G2snSV(ikl,isn)
      agsnSV(ikl,isn) =     Agrege(ikl) *agsnSV(ikl,isn-1) &
            + (1.-Agrege(ikl))*agsnSV(ikl,isn)
      isnoSV(ikl)     = min(Agrege(ikl) +isnoSV(ikl),real(nsno))
      iiceSV(ikl)     =                  iiceSV(ikl) &
            +     Agrege(ikl) *max(0,sign(1,iiceSV(ikl) &
            -isn +icemix)) &
            *max(0,sign(1,iiceSV(ikl) &
            -1          ))
  END DO


  ! +--Constrains Agregation in case of too much  Layers
  ! +  =================================================

  ! +--Search the thinest   non-zero Layer
  ! +  -----------------------------------



    DO ikl=1,knonv
      dzthin(ikl) =   0.                            ! Arbitrary unrealistic
    END DO                                          !       Layer Thickness
  DO ikl=1,knonv
    DO isn=1,isnoSV(ikl)-3 ! no agregation of 3 first snowlayers
                           ! ! XF 04/07/2019

      isno_n    =             isnoSV(ikl)-isn+1     ! Snow Normal.Profile
      iice_n    =             iiceSV(ikl)-isn       ! Ice  Normal.Profile
      iiceOK    = min(1,max(0,iice_n         +1))   ! Ice         Switch
  ! #vz     dz_ref(isn) =                                 !
  ! #vz.          dz_min *((1-iiceOK)*isno_n*isno_n       ! Theoretical Profile
  ! #vz.                 +    iiceOK *    2**iice_n)      !
  ! #vz.               /max(1,isnoSV(ikl))                !
      dz_dif      =     dz_min & ! Actual      Profile
            - dzsnSV(ikl    ,isn) & !
            /max(epsi,((1-iiceOK)*isno_n*isno_n & ! Theoretical Profile
            +    iiceOK *2.   **iice_n))     !
  ! #vz     dzwdif(isn) =     dz_dif                      !
      OKthin      = max(zero, & !
            sign(unun, & !
            dz_dif  - dzthin(ikl))) & ! 1.=> New thinest Lay.
            * max(0, & ! 1 => .le. isnoSV
            min(1, & !
            isnoSV(ikl)-isn +1 ))       !
      i_thin(ikl) =    (1. - OKthin) * i_thin(ikl) & ! Update   thinest Lay.
            + OKthin  * isn          !                Index
      dzthin(ikl) =    (1. - OKthin) * dzthin(ikl) & !
            + OKthin  * dz_dif       !


    END DO
  END DO





  ! +--Index of the contiguous Layer to agregate
  ! +  -----------------------------------------

  ! +   ***************
  call SISVAT_zCr
  ! +   ***************


  ! +--Assign the 2 Layers to agregate
  ! +  -------------------------------

    DO ikl=1,knonv
      isn         =    i_thin(ikl)
      if(LIndsv(ikl)>0) isn=min(isn, nsno-1) !cXF
      isagr1(ikl) =    istoSV(ikl,isn)
      isagr2(ikl) =    istoSV(ikl,isn+LIndsv(ikl))
      dzagr1(ikl) =    dzsnSV(ikl,isn)
      dzagr2(ikl) =    dzsnSV(ikl,isn+LIndsv(ikl))
      T_agr1(ikl) =    TsisSV(ikl,isn)
      T_agr2(ikl) =    TsisSV(ikl,isn+LIndsv(ikl))
      roagr1(ikl) =    ro__SV(ikl,isn)
      roagr2(ikl) =    ro__SV(ikl,isn+LIndsv(ikl))
      etagr1(ikl) =    eta_SV(ikl,isn)
      etagr2(ikl) =    eta_SV(ikl,isn+LIndsv(ikl))
      G1agr1(ikl) =    G1snSV(ikl,isn)
      G1agr2(ikl) =    G1snSV(ikl,isn+LIndsv(ikl))
      G2agr1(ikl) =    G2snSV(ikl,isn)
      G2agr2(ikl) =    G2snSV(ikl,isn+LIndsv(ikl))
      agagr1(ikl) =    agsnSV(ikl,isn)
      agagr2(ikl) =    agsnSV(ikl,isn+LIndsv(ikl))
      LstLay      = min(1,max(  0,    isnoSV(ikl)-1   ))
      Agrege(ikl) = min(1, &
            max(0, &
            NLaysv(ikl)   +isnoSV(ikl)-nsno &
            +NLay_s(ikl)                    ) &
            *LstLay           )

  ! + minimum uppermost layer thickness to guarantee a correct reproduction of the snow
  ! + atmosphere coupling
    if(dzsnSV(ikl,max(1,isnoSV(ikl)-0))>0.02 .or. & ! surface layers> 2-5-10
          dzsnSV(ikl,max(1,isnoSV(ikl)-1))>0.05 .or. & ! XF 04/07/2019
          dzsnSV(ikl,max(1,isnoSV(ikl)-2))>0.10 .or. &
          dzsnSV(ikl,max(1,isnoSV(ikl)-3))>0.30 )then
      Agrege(ikl) = min(1, &
            max(0, &
            NLaysv(ikl)   +isnoSV(ikl)+1-nsno & ! nsno-1 layers ma
            +NLay_s(ikl)                    ) &
            *LstLay           )
    endif

      isnoSV(ikl) =                    isnoSV(ikl) &
            -(1-LstLay)*max(zero, &
            sign(unun,      eps_21 &
            -dzsnSV(ikl,1)   ))
      isnoSV(ikl) =max(   0,           isnoSV(ikl)      )

      WEagre(ikl) = 0.
    END DO

    DO isn=1,nsno
    DO ikl=1,knonv
      WEagre(ikl) = WEagre(ikl) + ro__SV(ikl,isn)*dzsnSV(ikl,isn) &
            *min(1,max(0,i_thin(ikl)+1-isn))
    ENDDO
    ENDDO

  ! +--Agregates
  ! +  ---------

  ! +     ***************
    call SISVAT_zAg &
          (isagr1,isagr2,WEagre &
          ,dzagr1,dzagr2,T_agr1,T_agr2 &
          ,roagr1,roagr2,etagr1,etagr2 &
          ,G1agr1,G1agr2,G2agr1,G2agr2 &
          ,agagr1,agagr2,Agrege &
          )
  ! +     ***************


  ! +--Rearranges the Layers
  ! +  ---------------------

  ! +--New (agregated) Snow layer
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
    DO ikl=1,knonv
      isn     =             i_thin(ikl)
      isn     = min(isn,isn+LIndsv(ikl))
      isnoSV(ikl) =  max(0.,isnoSV(ikl) -Agrege(ikl))
      iiceSV(ikl) =         iiceSV(ikl) &
            -max(0,sign(1,iiceSV(ikl) -isn +icemix)) &
            *Agrege(ikl) &
            *max(0,sign(1,iiceSV(ikl) -1          ))
      istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn) &
            +   Agrege(ikl) *isagr1(ikl)
      dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn) &
            +   Agrege(ikl) *dzagr1(ikl)
      TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn) &
            +   Agrege(ikl) *T_agr1(ikl)
      ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn) &
            +   Agrege(ikl) *roagr1(ikl)
      eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn) &
            +   Agrege(ikl) *etagr1(ikl)
      G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn) &
            +   Agrege(ikl) *G1agr1(ikl)
      G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn) &
            +   Agrege(ikl) *G2agr1(ikl)
      agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn) &
            +   Agrege(ikl) *agagr1(ikl)
    END DO

  ! +--Above
  ! +  ^^^^^
    DO ikl=1,knonv
      isn1(ikl)=max(i_thin(ikl),i_thin(ikl)+LIndsv(ikl))
    END DO
    DO i=  1,nsno-1
    DO ikl=1,knonv
        staggr        =  min(1,max(0,i +1 -isn1(ikl)   ))
        istoSV(ikl,i) = (1.-staggr     )*istoSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*istoSV(ikl,i  ) &
              +   Agrege(ikl) *istoSV(ikl,i+1))
        dzsnSV(ikl,i) = (1.-staggr     )*dzsnSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*dzsnSV(ikl,i  ) &
              +   Agrege(ikl) *dzsnSV(ikl,i+1))
        TsisSV(ikl,i) = (1.-staggr     )*TsisSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*TsisSV(ikl,i  ) &
              +   Agrege(ikl) *TsisSV(ikl,i+1))
        ro__SV(ikl,i) = (1.-staggr     )*ro__SV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*ro__SV(ikl,i  ) &
              +   Agrege(ikl) *ro__SV(ikl,i+1))
        eta_SV(ikl,i) = (1.-staggr     )*eta_SV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*eta_SV(ikl,i  ) &
              +   Agrege(ikl) *eta_SV(ikl,i+1))
        G1snSV(ikl,i) = (1.-staggr     )*G1snSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*G1snSV(ikl,i  ) &
              +   Agrege(ikl) *G1snSV(ikl,i+1))
        G2snSV(ikl,i) = (1.-staggr     )*G2snSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*G2snSV(ikl,i  ) &
              +   Agrege(ikl) *G2snSV(ikl,i+1))
        agsnSV(ikl,i) = (1.-staggr     )*agsnSV(ikl,i  ) &
              + staggr*((1.-Agrege(ikl))*agsnSV(ikl,i  ) &
              +   Agrege(ikl) *agsnSV(ikl,i+1))
    END DO
    END DO

    DO ikl=1,knonv
      isn             = min(isnoSV(ikl) +1,nsno)
      istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn)
      dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn)
      TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn)
      ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn)
      eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn)
      G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn)
      G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn)
      agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn)
    END DO


  ! +--Search new Ice/Snow Interface (option II in MAR)
  ! +  ===============================================

    IF (ok_zsn_ii) THEN

    DO ikl=1,knonv
      iiceSV(ikl) =  0
    END DO

    DO ikl=1,knonv
    DO   isn=1,isnoSV(ikl)
      OK_ICE      = max(zero,sign(unun,ro__SV(ikl,isn)-ro_ice+20.)) &
            * max(zero,sign(unun,dzsnSV(ikl,isn)-epsi))
      iiceSV(ikl) = (1.-OK_ICE)       *iiceSV(ikl) &
            +     OK_ICE        *isn
    END DO
    END DO

    END IF

  return
end subroutine sisvat_zsn
