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

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

      IMPLICIT NONE

 
C +--Internal Variables
C +  ==================
 
      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
C +                                           ! 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

 
C +--DATA
C +  ====
 
      data      icemix /   0    /             ! 0 ====> Agregated Snow+Ice=Snow
      data      dzepsi /   0.0020/            ! Min single Layer Thickness
      data      dzxmin /   0.0025/            ! Min accept.Layer Thickness
c #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
C +   CAUTION:  dz_max > dz_min*2 is required ! Otherwise re-agregation is
C +                                           ! activated  after splitting
 

 

 
C +--Constrains Agregation         of too thin  Layers
C +  =================================================
 
C +--Search the thinest  non-zero Layer
C +  ----------------------------------
 
        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
cXF
      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))        !
C +                                                     !
     .             +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
C +                                                     !
     .                 +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
 


 
C +   ***************
      CALL SISVAT_zCr
C +   ***************
 
 
C +--Assign the 2 Layers to agregate
C +  -------------------------------
 
        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
 
 
C +--Agregates
C +  ---------
 
C +     ***************
        CALL SISVAT_zAg
     .                 (isagr1,isagr2,WEagre
     .                 ,dzagr1,dzagr2,T_agr1,T_agr2
     .                 ,roagr1,roagr2,etagr1,etagr2
     .                 ,G1agr1,G1agr2,G2agr1,G2agr2
     .                 ,agagr1,agagr2,Agrege
     .                 )
C +     ***************
 
 
C +--Rearranges the Layers
C +  ---------------------
 
C +--New (agregated) Snow layer
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
        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
 
C +--Above
C +  ^^^^^
        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
 

 
 
C +--Constrains Splitting          of too thick Layers
C +  =================================================
 
 
C +--Search the thickest non-zero Layer
C +  ----------------------------------
 
        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
 

C +--Rearranges the Layers
C +  ---------------------
 
      DO isn=nsno,2,-1
      DO ikl=1,knonv
        IF (Agrege(ikl)>0..AND.i_thin(ikl)<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
 
 
C +--Constrains Agregation in case of too much  Layers
C +  =================================================
 
C +--Search the thinest   non-zero Layer
C +  -----------------------------------


 
        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
 
 



C +--Index of the contiguous Layer to agregate
C +  -----------------------------------------
 
C +   ***************
      CALL SISVAT_zCr
C +   ***************
 
 
C +--Assign the 2 Layers to agregate
C +  -------------------------------
 
        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           )
 
C + minimum uppermost layer thickness to guarantee a correct reproduction of the snow
C + 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

C +--Agregates
C +  ---------
 
C +     ***************
        CALL SISVAT_zAg
     .                 (isagr1,isagr2,WEagre
     .                 ,dzagr1,dzagr2,T_agr1,T_agr2
     .                 ,roagr1,roagr2,etagr1,etagr2
     .                 ,G1agr1,G1agr2,G2agr1,G2agr2
     .                 ,agagr1,agagr2,Agrege
     .                 )
C +     ***************
 
 
C +--Rearranges the Layers
C +  ---------------------
 
C +--New (agregated) Snow layer
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
        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
 
C +--Above
C +  ^^^^^
        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


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

        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
