SUBROUTINE SISVAT_zCr
  ! +
  ! +------------------------------------------------------------------------+
  ! | MAR          SISVAT_zCr                                12-12-2002  MAR |
  ! |   SubRoutine SISVAT_zCr determines criteria for Layers Agregation      |
  ! |                                                                        |
  ! +------------------------------------------------------------------------+
  ! |                                                                        |
  ! |   PARAMETERS:  klonv: Total Number of columns =                        |
  ! |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
  ! |                     X       Number of Mosaic Cell per grid box         |
  ! |                                                                        |
  ! |   INPUT /  isnoSV   = total Nb of Ice/Snow Layers                      |
  ! |   OUTPUT:  iiceSV   = total Nb of Ice      Layers                      |
  ! |   ^^^^^^   ispiSV   = 0,...,nsno: Uppermost Superimposed Ice Layer     |
  ! |            istoSV   = 0,...,5 :   Snow     History (see istdSV data)   |
  ! |                                                                        |
  ! |   INPUT /  ro__SV   : Soil/Snow Volumic Mass                   [kg/m3] |
  ! |   OUTPUT:           & Snow     Temperatures (layers  1,2,...,nsno) [K] |
  ! |   ^^^^^^   G1snSV   : Dendricity (<0) or Sphericity (>0) of Snow Layer |
  ! |            G2snSV   : Sphericity (>0) or Size            of Snow Layer |
  ! |            agsnSV   : Snow       Age                             [day] |
  ! |                                                                        |
  ! |   OUTPUT:  LIndsv   : Relative Index of a contiguous Layer to agregate |
  ! |   ^^^^^^                                                               |
  ! +------------------------------------------------------------------------+
  ! +
  ! +
  ! +
  ! +
  ! +--Global Variables
  ! +  ================
  !
  use VARphy
  use VAR_SV
  use VARdSV
  use VAR0SV
  use VARxSV
  use VARySV
  use VARtSV

  IMPLICIT NONE

  ! +
  ! +
  ! +--Internal Variables
  ! +  ==================
  ! +
  integer :: ikl   ,isn   ,is0   ,is1
  integer :: isno_1                        ! Switch:  ! Snow Layer over Ice
  real :: Dtyp_0,Dtyp_1                 ! Snow Grains Difference Measure
  real :: DenSph                        ! 1. when contiguous spheric
  ! +                                           !     and dendritic  Grains
  real :: DendOK                        ! 1. when dendritic  Grains
  real :: dTypMx                        ! Grain Type Differ.
  real :: dTypSp                        ! Sphericity Weight
  real :: dTypRo                        ! Density    Weight
  real :: dTypDi                        ! Grain Diam.Weight
  real :: dTypHi                        ! History    Weight


  ! +--DATA
  ! +  ====

  data      dTypMx / 200.0  /             ! Grain Type Weight
  data      dTypSp /   0.5  /             ! Sphericity Weight
  data      dTypRo /   0.5  /             ! Density    Weight
  data      dTypDi /  10.0  /             ! Grain Diam.Weight
  data      dTypHi / 100.0  /             ! History    Weight


  ! +--Agregation Criteria
  ! +  ===================
  ! +
  DO  ikl=1,knonv
      i_thin(ikl) = min(i_thin(ikl),isnoSV(ikl))
      isn         = max(1          ,i_thin(ikl))
  ! +
  ! +
  ! +--Comparison with the downward Layer
  ! +  ----------------------------------
  ! +

      is0    = max(1,        i_thin(ikl)-1 )        ! Downward Layer Index
      DenSph = max(zero, & ! isn/is1
            sign(unun, & ! Dendricity/Sphericity
            epsi-G1snSV(ikl,isn) & !            Switch
            *G1snSV(ikl,is0)))      !
      DendOK = max(zero, & ! Dendricity Switch
            sign(unun, & !
            epsi-G1snSV(ikl,isn)))      !
  ! +
      Dtyp_0 = &
            DenSph *      dTypMx &
            +(1.-DenSph) &
            *    DendOK *((abs(G1snSV(ikl,isn) & ! Dendricity
            -G1snSV(ikl,is0)) & !     Contribution
            +abs(G2snSV(ikl,isn) & ! Sphericity
            -G2snSV(ikl,is0))) *dTypSp & !     Contribution
            +abs(ro__SV(ikl,isn) & ! Density
            -ro__SV(ikl,is0))  *dTypRo) & !     Contribution
            +(1.-DenSph) & !
            *(1.-DendOK)*((abs(G1snSV(ikl,isn) & ! Sphericity
            -G1snSV(ikl,is0)) & !     Contribution
            +abs(G2snSV(ikl,isn) & ! Size
            -G2snSV(ikl,is0))) *dTypDi & !     Contribution
            +abs(ro__SV(ikl,isn) & ! Density
            -ro__SV(ikl,is0))  *dTypRo) !     Contribution
      Dtyp_0 = & !
            min(dTypMx, & !
            Dtyp_0 & !
            +abs(istoSV(ikl,isn) & ! History
            -istoSV(ikl,is0))  *dTypHi) & !     Contribution
            +             (1 -abs(isn-is0))  * 1.e+6 & !"Same Layer"Score
            +  max(0,1-abs(iiceSV(ikl) & !"Ice /Snow
            -is0))  * 1.e+6  ! Interface" Score
  ! +
  ! +
  ! +--Comparison with the   upward Layer
  ! +  ----------------------------------
  ! +
      is1    = min(          i_thin(ikl)+1, & ! Upward   Layer Index
            max(1,    isnoSV(ikl)  ))        !
      DenSph = max(zero, & ! isn/is1
            sign(unun, & ! Dendricity/Sphericity
            epsi-G1snSV(ikl,isn) & !            Switch
            *G1snSV(ikl,is1)))      !
      DendOK = max(zero, & ! Dendricity Switch
            sign(unun, & !
            epsi-G1snSV(ikl,isn)))      !
  ! +
      Dtyp_1 = &
            DenSph *      dTypMx &
            +(1.-DenSph) &
            *    DendOK *((abs(G1snSV(ikl,isn) & ! Dendricity
            -G1snSV(ikl,is1)) & !     Contribution
            +abs(G2snSV(ikl,isn) & ! Sphericity
            -G2snSV(ikl,is1))) *dTypSp & !     Contribution
            +abs(ro__SV(ikl,isn) & ! Density
            -ro__SV(ikl,is1))  *dTypRo) & !     Contribution
            +(1.-DenSph) & !
            *(1.-DendOK)*((abs(G1snSV(ikl,isn) & ! Sphericity
            -G1snSV(ikl,is1)) & !     Contribution
            +abs(G2snSV(ikl,isn) & ! Size
            -G2snSV(ikl,is1))) *dTypDi & !     Contribution
            +abs(ro__SV(ikl,isn) & ! Density
            -ro__SV(ikl,is1))  *dTypRo) !     Contribution
      Dtyp_1 = & !
            min(dTypMx, & !
            Dtyp_1 & !
            +abs(istoSV(ikl,isn) & ! History
            -istoSV(ikl,is1))  *dTypHi) & !     Contribution
            +             (1 -abs(isn-is1))  * 1.e+6 & !"Same Layer"Score
            +  max(0,1-abs(iiceSV(ikl) & !"Ice /Snow
            -isn))  * 1.e+6  ! Interface" Score
  ! +
  ! +
  ! +--Index of the Layer to agregate
  ! +  ==============================
  ! +
      LIndsv(ikl) = sign(unun,Dtyp_0 &
            -Dtyp_1)
      isno_1      = (1 -min (abs(isnoSV(ikl) & ! Switch = 1
            -iiceSV(ikl)-1),1)) & !   if isno = iice +1
            * (1 -min (abs(isnoSV(ikl) & ! Switch = 1
            -i_thin(ikl)  ),1)) !   if isno = i_ithin
      LIndsv(ikl) = (1 -isno_1) *LIndsv(ikl) & ! Contiguous Layer is
            -isno_1                      ! downward for top L.
      i_thin(ikl) =  max(1,   i_thin(ikl)   )
  END DO
  ! +

end subroutine sisvat_zcr