MODULE module_scientific
! Module of the scientific function/subroutines
  
!!!!!!! Functions & subroutines
! all_polygons_properties: Subroutine to determine the properties of all polygons in a 2D field:
! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.)
! coincidence_all_polys: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
!   to a map of integer polygons
! coincidence_all_polys_area: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
!   to a map of integer polygons filtrered by a minimal area
! coincidence_poly: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
!   to a map of integer polygons
! coincidence_poly_area: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
!   to a map of integer polygons filtrered by a minimal area
! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon 
!   following ray casting algorithm
! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders 
!   (limits of a region)
! paths_border: Subroutine to search the paths of a border field.
! path_properties: Subroutine to determine the properties of a path
! polygon_properties: Subroutine to determine the properties of a polygon (as .TRUE. matrix)
! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1!
! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
! poly_overlap_tracks: Subroutine to determine tracks of a series of consecutive 2D field with polygons using 
!   maximum overlaping/coincidence! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
! poly_overlap_tracks_area: Subroutine to determine tracks of a series of consecutive 2D field with polygons using 
!   maximum overlaping/coincidence filtrered by a minimal area
! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k'
! rand_sample: Subroutine to randomly sample a range of indices
! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
! StatsR_K: Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a 
!   series of r_k numbers
! SwapR_K*: Subroutine swaps the values of its two formal arguments.

!!! *Functions/Subroutines to sort values adpated. The method used is usually referred to as "selection" method.
! from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90

  USE module_definitions
  USE module_generic

  CONTAINS


  SUBROUTINE poly_overlap_tracks_area(dbg, dx, dy, dt, minarea, inNallpolys, allpolys, ctrpolys,      &
    areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks)
! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum 
!   overlaping/coincidence filtrered by a minimal area

    IMPLICIT NONE

    LOGICAL, INTENT(in)                                  :: dbg
    INTEGER, INTENT(in)                                  :: dx, dy, dt, Nmaxpoly, Nmaxtracks
    INTEGER, DIMENSION(dt), INTENT(in)                   :: inNallpolys
    INTEGER, DIMENSION(dx,dy,dt), INTENT(in)             :: allpolys
    REAL(r_k), INTENT(in)                                :: minarea
    REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in)      :: ctrpolys
    REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in)        :: areapolys
    REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,dt),                                                   &
      INTENT(out)                                        :: tracks
    REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out)   :: finaltracks

! Local
    INTEGER                                              :: i, j, ip, it, iip, itt
    INTEGER                                              :: ierr
    REAL(r_k), DIMENSION(2,Nmaxpoly,dt)                  :: Actrpolys
    REAL(r_k), DIMENSION(Nmaxpoly,dt)                    :: Aareapolys
    INTEGER, DIMENSION(dt)                               :: Nallpolys
    INTEGER, DIMENSION(dx,dy)                            :: meetpolys, searchpolys
    INTEGER, DIMENSION(Nmaxpoly)                         :: coincidencies
    INTEGER, DIMENSION(Nmaxpoly)                         :: prevID, currID
    INTEGER, DIMENSION(:), ALLOCATABLE                   :: coins
    INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: coinsNpts
    INTEGER, DIMENSION(Nmaxpoly,dt)                      :: polycoincidencies
    INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt)             :: coincidenciesNpts
    INTEGER                                              :: Nmeet, Nsearch, Nindep
    INTEGER, DIMENSION(dt)                               :: Nindeppolys
    CHARACTER(len=5)                                     :: NcoinS
    INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt)             :: polysIndep
    INTEGER, DIMENSION(Nmaxpoly,dt)                      :: NpolysIndep
    INTEGER, DIMENSION(Nmaxpoly,dt)                      :: SpolysIndep
    INTEGER                                              :: iindep, iiprev
    INTEGER                                              :: Nprev, NNprev, Ntprev
    LOGICAL                                              :: Indeppolychained
    INTEGER                                              :: itrack, ictrack
    INTEGER                                              :: ixp, iyp, ttrack
    INTEGER, DIMENSION(dt)                               :: Ntracks
    INTEGER                                              :: idtrack, maxtrack

!!!!!!! Variables
! dx,dy,dt: space/time dimensions
! Nallpolys: Vector with the number of polygons at each time-step
! allpolys: Series of 2D field with the polygons
! minarea: minimal area (in same units as areapolys) to perform the tracking 
! ctrpolys: center of the polygons
! areapolys: area of the polygons
! Nmaxpoly: Maximum possible number of polygons
! Nmaxtracks: maximum number of tracks
! tracks: series of consecutive polygons
! trackperiod: period of the track in time-steps

    fname = 'poly_overlap_tracks_area'

    IF (dbg) PRINT *,TRIM(fname)

    polycoincidencies = fillvalI
    coincidenciesNpts = fillvalI
    ! Number of independent polygons by time step
    Nindeppolys = 0
    ! Number of polygons attached to each independent polygons by time step
    NpolysIndep = 0
    ! ID of searching polygon attached to each independent polygons by time step
    SpolysIndep = 0
    ! ID of polygons attached to each independent polygons by time step
    polysIndep = 0
    ! ID of polygons from previous time-step
    prevID = 0
    ! ID of polygons from current time-step
    currID = 0

    ! First time-step all are independent polygons
    it = 1
    Nmeet = inNallpolys(it)
    Nindeppolys(it) = Nmeet
    ip = 0
    meetpolys = allpolys(:,:,it)
    DO i=1, Nmeet
      IF (areapolys(i,it) >= minarea) THEN
        ip = ip + 1
        SpolysIndep(ip,it) = i
        currID(ip) = i
        NpolysIndep(ip,it) = 1
        polysIndep(ip,1,it) = i
        Actrpolys(1,ip,it) = ctrpolys(1,i,it)
        Actrpolys(2,ip,it) = ctrpolys(2,i,it)
        Aareapolys(ip,it) = areapolys(i,it)
      ELSE
        WHERE (meetpolys == i)
          meetpolys = 0
        END WHERE
      END IF
    END DO
    Nallpolys(1) = ip
    Nindeppolys(1) = ip

    ! Starting step
    it = 0
    IF (dbg) THEN
      PRINT *,'  time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',0
      PRINT *,'    number of independent polygons:', Nindeppolys(it+1)
      PRINT *,'    indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
      DO i=1, Nindeppolys(it+1)
        PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':',                                     &
          polysIndep(i,1:NpolysIndep(i,it+1),it+1)
      END DO
    END IF

    ! Looking for the coincidencies at each time step
    DO it=1, dt-1
      ! Number of times that a polygon has a coincidence
      coincidencies = 0

      Nmeet = inNallpolys(it+1)
      searchpolys = meetpolys
      meetpolys = allpolys(:,:,it+1)
      prevID = 0
      prevID = currID
      ip = 0

      DO i=1, Nmeet
        IF (areapolys(i,it+1) >= minarea) THEN
          ip = ip + 1
          currID(ip) = i
          Actrpolys(1,ip,it+1) = ctrpolys(1,i,it+1)
          Actrpolys(2,ip,it+1) = ctrpolys(2,i,it+1)
          Aareapolys(ip,it+1) = areapolys(i,it+1)
        ELSE
          WHERE (meetpolys == i)
            meetpolys = 0
          END WHERE
        END IF
      END DO
      Nallpolys(it+1) = ip
      Nindeppolys(it+1) = ip

      Nmeet = Nallpolys(it+1)
      ! Looking throughout the independent polygons
      Nsearch = Nindeppolys(it)

      IF (ALLOCATED(coins)) DEALLOCATE(coins)
      ALLOCATE(coins(Nmeet), STAT=ierr)
      msg="Problems allocating 'coins'"
      CALL ErrMsg(msg,fname,ierr)

      IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts)
      ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr)
      msg="Problems allocating 'coinsNpts'"
      CALL ErrMsg(msg,fname,ierr)

      CALL coincidence_all_polys_area(dbg, dx,dy, Nmeet, currID, meetpolys, Actrpolys(:,1:Nmeet,it+1),&
        Nsearch, prevID, searchpolys, Actrpolys(:,1:Nsearch,it), Aareapolys(1:Nsearch,it), coins,     &
        coinsNpts)

      polycoincidencies(1:Nmeet,it+1) = coins
      coincidenciesNpts(1:Nmeet,1:Nsearch,it+1) = coinsNpts

      ! Counting the number of times that a polygon has a coincidency
      IF (dbg) THEN
        PRINT *,'  Coincidencies for the given time-step:', it+1,' _______'
        DO i=1, Nmeet
          PRINT *,currID(i), coins(i),' N search pts:', coinsNpts(i,:)
        END DO
      END IF

      ! Looking for the same equivalencies
      Nindep = 0
      DO i=1, Nmeet
        IF (coins(i) == -1) THEN
          Nindep = Nindep + 1
          SpolysIndep(Nindep,it+1) = -1
          NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1
          polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = currID(i)
        ELSE IF (coins(i) == -9) THEN
          WRITE(NcoinS,'(I5)')coins(i)
          msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " //      &
            "coincidence of polygon"
          CALL ErrMsg(msg, fname, -1)
        ELSE
          ! Looking for coincidences with previous independent polygons
          DO ip=1, Nsearch
            ! Looking into the polygons associated
            NNprev = NpolysIndep(ip,it)
            DO j=1, NNprev
              IF (coins(i) == polysIndep(ip,j,it)) THEN
                IF (coincidencies(ip) == 0) THEN
                  Nindep = Nindep + 1
                  SpolysIndep(Nindep,it+1) = coins(i)
                END IF
                ! Which index corresponds to this coincidence?
                iindep = Index1DArrayI(SpolysIndep(1:Nmeet,it+1), Nindep, coins(i))
                IF (iindep < 0) THEN
                  msg = 'Wrong index! There must be an index here'
                  CALL ErrMsg(msg,fname,iindep)
                END IF
                coincidencies(ip) = coincidencies(ip) + 1
                NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1
                polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = currID(i)
                EXIT
              END IF
            END DO
          END DO
        END IF
      END DO
      Nindeppolys(it+1) = Nindep

      IF (dbg) THEN
        PRINT *,'  time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch
        PRINT *,'    number of independent polygons:', Nindeppolys(it+1)
        PRINT *,'    indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
        DO i=1, Nindeppolys(it+1)
          PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':',                                   &
            polysIndep(i,1:NpolysIndep(i,it+1),it+1)
        END DO
      END IF
    END DO

    IF (dbg) THEN
      PRINT *,  'Coincidencies to connect _______'
      DO it=1, dt
        PRINT *,'  it:', it, ' Nindep:', Nindeppolys(it)
        PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs'
        DO ip=1, Nindeppolys(it)
          PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':',                         &
            polysIndep(ip,1:NpolysIndep(ip,it),it)
        END DO
      END DO

    END IF

    ! Trajectories
    ! It should be done following the number of 'independent' polygons
    ! One would concatenate that independent polygons which share IDs from one step to another

    ! First time-step. Take all polygons
    itrack = 0
    tracks = 0.
    Ntracks = 0
    DO ip=1, Nindeppolys(1)
      itrack = itrack + 1
      tracks(1,1,itrack,1) = itrack*1.
      tracks(2,1,itrack,1) = SpolysIndep(ip,1)
      tracks(3,1,itrack,1) = ctrpolys(1,ip,1)
      tracks(4,1,itrack,1) = ctrpolys(2,ip,1)
      tracks(5,1,itrack,1) = 1
      Ntracks(1) = Ntracks(1) + 1
    END DO

    ! Looping allover already assigned tracks
    timesteps: DO it=2, dt
      IF (dbg) PRINT *,'track-timestep:', it, 'N indep polys:', Nindeppolys(it)
      ! Indep polygons current time-step
      current_poly: DO i=1, Nindeppolys(it)
        IF (dbg) PRINT *,'  curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:',   &
          NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it)
        Indeppolychained = .FALSE.

        ! Number of tracks previous time-step
        ! Looping overall
        it1_tracks: DO itt=1, Ntracks(it-1)
          itrack = tracks(1,1,itt,it-1)
          ! Number polygons ID assigned
          Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0)
          IF (dbg) PRINT *,itt,'  track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1)

          ! Looking for coincidencies
          DO iip=1, Ntprev
            IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN
              Indeppolychained = .TRUE.
              IF (dbg) PRINT *,'    coincidence found by polygon:', tracks(2,iip,itt,it-1)
              EXIT
            END IF
          END DO
          IF (Indeppolychained) THEN
            Ntracks(it) = Ntracks(it) + 1
            ictrack = Ntracks(it)
            ! Assigning all the IDs to the next step of the track
            DO iip=1, NpolysIndep(i,it)
              iiprev = polysIndep(i,iip,it)
              tracks(1,iip,ictrack,it) = itrack
              tracks(2,iip,ictrack,it) = iiprev
              ixp = ctrpolys(1,iiprev,it)
              iyp = ctrpolys(2,iiprev,it)
              tracks(3,iip,ictrack,it) = ixp
              tracks(4,iip,ictrack,it) = iyp
              tracks(5,iip,ictrack,it) = it
            END DO
            EXIT
          END IF
          IF (Indeppolychained) EXIT
        END DO it1_tracks

        ! Creation of a new track
        IF (.NOT.Indeppolychained) THEN
          Ntracks(it) = Ntracks(it) + 1
          ictrack = Ntracks(it)
          ! ID of new track
          maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.))
          IF (dbg) PRINT *,'  New track!', maxtrack+1

          ! Assigning all the IDs to the next step of the track
          DO j=1, NpolysIndep(i,it)
            iiprev = polysIndep(i,j,it)
            tracks(1,j,ictrack,it) = maxtrack+1
            tracks(2,j,ictrack,it) = iiprev
            ixp = ctrpolys(1,iiprev,it)
            iyp = ctrpolys(2,iiprev,it)
            tracks(3,j,ictrack,it) = ixp
            tracks(4,j,ictrack,it) = iyp
            tracks(5,j,ictrack,it) = it
          END DO
        END IF

      END DO current_poly

      IF (dbg) THEN
        PRINT *,'  At this time-step:', it, ' N trajectories:', Ntracks(it)
        DO i=1, Ntracks(it)
          Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0)
          PRINT *,i ,'ID tracks:', tracks(1,1,i,it), 'ID polygon:', tracks(2,1:Nprev,i,it)
        END DO
      END IF

    END DO timesteps

    ! Summarizing trajectories
    ! When multiple polygons are available, the mean of their central positions determines the position

    finaltracks = 0.
    maxtrack = MAXVAL(tracks(1,:,:,:))

    DO it=1, dt
      DO itt=1, Ntracks(it)
        itrack = INT(tracks(1,1,itt,it))
        Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0)
        PRINT *,'it:', it,'itrack:', itrack, 'Nprev:', Nprev
        PRINT *,'  coordinates ix iy _______'
        DO i=1,Nprev
          PRINT *,'    ',tracks(3,i,itt,it),tracks(4,i,itt,it)
        END DO
        finaltracks(1,itrack,it) = itrack*1.
        finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev
        finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev
        finaltracks(4,itrack,it) = it*1.
        PRINT *,'  finaltrack:', finaltracks(:,itrack,it)
      END DO
    END DO

    RETURN

  END SUBROUTINE poly_overlap_tracks_area

  SUBROUTINE coincidence_all_polys_area(dbg, dx, dy, Nallpoly, IDallpoly, allpoly, icpolys, Npoly, &
    IDpolys, polys, cpolys, apolys, polycoins, coinNptss)
! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
!   In case of multiple coincidencies, the closest and then the biggest is taken filtrered by a minimal area
!   Here the difference is that the index does not coincide with its ID

    IMPLICIT NONE

    LOGICAL, INTENT(in)                                  :: dbg
    INTEGER, INTENT(in)                                  :: dx, dy, Nallpoly, Npoly
    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: allpoly, polys
    INTEGER, DIMENSION(Nallpoly), INTENT(in)             :: IDallpoly
    INTEGER, DIMENSION(Npoly), INTENT(in)                :: IDpolys
    REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in)         :: icpolys
    REAL(r_k), DIMENSION(2,Npoly), INTENT(in)            :: cpolys
    REAL(r_k), DIMENSION(Npoly), INTENT(in)              :: apolys
    INTEGER, DIMENSION(Nallpoly), INTENT(out)            :: polycoins
    INTEGER, DIMENSION(Nallpoly,Npoly), INTENT(out)      :: coinNptss

! Local
    INTEGER                                              :: i, j, ip
    INTEGER                                              :: maxcorr
    INTEGER                                              :: Nmaxcorr
    LOGICAL, DIMENSION(dx,dy)                            :: boolpoly
    INTEGER                                              :: maxcoin
    REAL                                                 :: dist, maxcoindist, maxcoinarea
    INTEGER, DIMENSION(Npoly)                            :: IDcoins

!!!!!!! Variables
! dx,dy: dimension of the space
! Nallpoly: Number of polygons to find coincidence
! allpoly: space with the polygons to meet
! icpolys: center of the polygons to look for the coincidence
! Npoly: number of polygons on the 2D space
! polys: 2D field of polygons identified by their integer number (0 for no polygon)
! cpolys: center of the polygons
! apolys: area of the polygons
! polycoins: coincident polyogn
!          -1: no-coincidence
!   1 < Npoly: single coincidence with a given polygon
!          -9: coincidence with more than one polygon
! coinNptss: number of points coincident with each polygon

    fname = 'coincidence_all_polys_area'
    IF (dbg) PRINT *,TRIM(fname)

    DO ip=1, Nallpoly
      boolpoly = allpoly == IDallpoly(ip)
      CALL coincidence_poly_area(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), IDcoins,         &
        coinNptss(ip,:))
      IF (dbg) PRINT *,'  polygon', IDallpoly(ip), ' coincidence with:', polycoins(ip), 'IDpolys:', IDpolys(1:Npoly)

      ! Coincidence with more than one polygon
      IF (polycoins(ip) == -9) THEN
        maxcoindist = -10.
        maxcoinarea = -10.
        maxcoin = MAXVAL(coinNptss(ip,:))
        DO j=1, Npoly
          IF (coinNptss(ip,j) == maxcoin) THEN
            dist = SQRT( (icpolys(1,ip)*1.-cpolys(1,j)*1.)**2 + (icpolys(2,ip)*1.-cpolys(2,j)*1.)**2 )
            IF ( dist > maxcoindist) THEN
              maxcoindist = dist
              maxcoinarea = apolys(j)
              polycoins(ip) = IDcoins(j)
            ELSE IF ( dist == maxcoindist) THEN
              IF (apolys(j) > maxcoinarea) THEN
                polycoins(ip) = IDcoins(j)
                maxcoinarea = apolys(j)
              END IF
            END IF 
          END IF
        END DO
      END IF
    END DO

    RETURN

  END SUBROUTINE coincidence_all_polys_area

  SUBROUTINE coincidence_poly_area(dbg, dx, dy, poly, Npoly, polys, polycoin, IDpoly, coinNpts)
! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
!   Here the difference is that the index does not coincide with its ID

    IMPLICIT NONE

    LOGICAL, INTENT(in)                                  :: dbg
    INTEGER, INTENT(in)                                  :: dx, dy, Npoly
    LOGICAL, DIMENSION(dx,dy), INTENT(in)                :: poly
    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: polys
    INTEGER, INTENT(out)                                 :: polycoin
    INTEGER, DIMENSION(Npoly), INTENT(out)               :: IDpoly, coinNpts

! Local
    INTEGER                                              :: i, j, ip
    INTEGER                                              :: maxcorr
    INTEGER                                              :: Nmaxcorr

!!!!!!! Variables
! dx,dy: dimension of the space
! poly: bolean polygon to meet
! Npoly: number of polygons on the 2D space
! polys: 2D field of polygons identified by their integer number (0 for no polygon)
! polycoin: coincident polyogn
!          -1: no-coincidence
!   1 < Npoly: single coincidence with a given polygon
!          -9: coincidence with more than one polygon
! IDpoly: ID of the found polygon
! coinNpts: number of points coincident with each polygon

    fname = 'coincidence_poly_area'
    IF (dbg) PRINT *,TRIM(fname)

    IF (dbg) THEN
      PRINT *,'  Boolean polygon to search coincidences ...'
      DO i=1,dx
        PRINT *,poly(i,:)
      END DO

      PRINT *,'  2D polygons space ...'
      DO i=1,dx
        PRINT '(1000(I3,1x))',polys(i,:)
      END DO
    END IF

    ! Looking for coincient points for the polygon
    coinNpts = 0
    IDpoly = 0
    ip = 0
    DO i=1,dx
      DO j=1,dy
        IF (poly(i,j) .AND. polys(i,j) .NE. 0) THEN
          IF (.NOT.ANY(IDpoly == polys(i,j))) THEN
            ip = ip + 1
            IDpoly(ip) = polys(i,j)
          ELSE
            ip = Index1DarrayI(IDpoly, Npoly, polys(i,j))
          END IF
          coinNpts(ip) = coinNpts(ip) + 1
        END IF
      END DO
    END DO

    maxcorr = 0
    maxcorr = INT(MAXVAL(coinNpts*1.))

    IF (dbg) PRINT *,'  Maximum coincidence:', maxcorr
    IF (maxcorr == 0) THEN
      polycoin = -1
    ELSE
      Nmaxcorr = 0
      DO ip=1, Npoly
        IF (coinNpts(ip) == maxcorr) THEN
          Nmaxcorr = Nmaxcorr+1
          polycoin = IDpoly(ip)
        END IF
      END DO
      IF (Nmaxcorr > 1) polycoin = -9
    END IF

    IF (dbg) THEN
      PRINT *,'  Coincidences for each polygon _______', Npoly
      DO ip=1, Npoly
        PRINT *,'  ',ip, ' ID:', IDpoly(ip),': ', coinNpts(ip)
      END DO
    END IF

    RETURN

END SUBROUTINE coincidence_poly_area

  SUBROUTINE poly_overlap_tracks(dbg, dx, dy, dt, minarea, Nallpolys, allpolys, ctrpolys,        &
    areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks)
! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum overlaping/coincidence

    IMPLICIT NONE

    LOGICAL, INTENT(in)                                  :: dbg
    INTEGER, INTENT(in)                                  :: dx, dy, dt, Nmaxpoly, Nmaxtracks
    INTEGER, DIMENSION(dt), INTENT(in)                   :: Nallpolys
    INTEGER, DIMENSION(dx,dy,dt), INTENT(in)             :: allpolys
    REAL(r_k), INTENT(in)                                :: minarea
    REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in)      :: ctrpolys
    REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in)        :: areapolys
    REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,dt),                                                   &
      INTENT(out)                                        :: tracks
    REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out)   :: finaltracks

! Local
    INTEGER                                              :: i, j, ip, it, iip, itt
    INTEGER                                              :: ierr
    INTEGER, DIMENSION(Nmaxpoly,dt)                      :: coincidencies, NOcoincidencies
    INTEGER, DIMENSION(:), ALLOCATABLE                   :: coins
    INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: coinsNpts
    INTEGER, DIMENSION(Nmaxpoly,dt)                      :: polycoincidencies
    INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt)             :: coincidenciesNpts
    INTEGER                                              :: Nmeet, Nsearch, Nindep
    INTEGER, DIMENSION(dt)                               :: Nindeppolys
    CHARACTER(len=5)                                     :: NcoinS
    INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt)             :: polysIndep
    INTEGER, DIMENSION(Nmaxpoly,dt)                      :: NpolysIndep
    INTEGER, DIMENSION(Nmaxpoly,dt)                      :: SpolysIndep
    INTEGER                                              :: iindep, iiprev
    INTEGER                                              :: Nprev, NNprev, Ntprev
    LOGICAL                                              :: Indeppolychained
    INTEGER                                              :: itrack, ictrack
    INTEGER                                              :: ixp, iyp, ttrack
    INTEGER, DIMENSION(dt)                               :: Ntracks
    INTEGER                                              :: idtrack, maxtrack

!!!!!!! Variables
! dx,dy,dt: space/time dimensions
! Nallpolys: Vector with the number of polygons at each time-step
! allpolys: Series of 2D field with the polygons
! minarea: minimal area (in same units as areapolys) to perform the tracking 
! ctrpolys: center of the polygons
! areapolys: area of the polygons
! Nmaxpoly: Maximum possible number of polygons
! Nmaxtracks: maximum number of tracks
! tracks: series of consecutive polygons
! trackperiod: period of the track in time-steps

    fname = 'poly_overlap_tracks'

    IF (dbg) PRINT *,TRIM(fname)

    polycoincidencies = fillvalI
    coincidenciesNpts = fillvalI
    ! Number of times that a polygon has a coincidence
    coincidencies = 0
    ! Polygons without a coincidence
    NOcoincidencies = 0
    ! Number of independent polygons by time step
    Nindeppolys = 0
    ! Number of polygons attached to each independent polygons by time step
    NpolysIndep = 0
    ! ID of searching polygon attached to each independent polygons by time step
    SpolysIndep = 0
    ! ID of polygons attached to each independent polygons by time step
    polysIndep = 0

    ! First time-step all are independent polygons
    it = 1
    Nmeet = Nallpolys(it)
    Nindeppolys(it) = Nmeet
    DO i=1, Nmeet
      SpolysIndep(i,it) = i
      NpolysIndep(1:Nmeet,it) = 1
      polysIndep(1,i,it) = i
    END DO

    ! Looking for the coincidencies at each time step
    DO it=1, dt-1
      Nmeet = Nallpolys(it+1)
      Nsearch = Nallpolys(it)

      IF (ALLOCATED(coins)) DEALLOCATE(coins)
      ALLOCATE(coins(Nmeet), STAT=ierr)
      msg="Problems allocating 'coins'"
      CALL ErrMsg(msg,fname,ierr)

      IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts)
      ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr)
      msg="Problems allocating 'coinsNpts'"
      CALL ErrMsg(msg,fname,ierr)

      CALL coincidence_all_polys(dbg, dx, dy, Nmeet, allpolys(:,:,it+1), ctrpolys(:,1:Nmeet,it+1),    &
        Nsearch, allpolys(:,:,it), ctrpolys(:,1:Nsearch,it), areapolys(1:Nsearch,it), coins, coinsNpts)

      polycoincidencies(1:Nmeet,it+1) = coins
      coincidenciesNpts(1:Nmeet,1:Nsearch,it+1) = coinsNpts

      ! Counting the number of times that a polygon has a coincidency
      IF (dbg) THEN
        PRINT *,'  Coincidencies for the given time-step:', it+1,' _______'
        DO i=1, Nmeet
          PRINT *,coins(i),' N search pts:', coinsNpts(i,:)
        END DO
      END IF

      Nindep = 0
      DO i=1, Nmeet
        IF (coins(i) == -1) THEN
          Nindep = Nindep + 1
          NOcoincidencies(i,it+1) = 1
          SpolysIndep(Nindep,it+1) = -1
          NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1
          polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = i
        ELSE IF (coins(i) == -9) THEN
          WRITE(NcoinS,'(I5)')coins(i)
          msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " //      &
            "coincidence of polygon"
          CALL ErrMsg(msg, fname, -1)
        ELSE
          DO ip=1, Nsearch
            IF (coins(i) == ip) THEN
              IF (coincidencies(ip,it+1) == 0) THEN
                Nindep = Nindep + 1
                SpolysIndep(Nindep,it+1) = ip
              END IF
              coincidencies(ip,it+1) = coincidencies(ip,it+1) + 1
              DO iindep=1, Nindep
                IF (SpolysIndep(iindep,it+1) == coins(i)) THEN
                  NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1
                  polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = i
                END IF
              END DO
            END IF
          END DO
        END IF
      END DO
      Nindeppolys(it+1) = Nindep

      IF (dbg) THEN
        PRINT *,'  time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch
        PRINT *,'    number of independent polygons:', Nindeppolys(it+1)
        PRINT *,'    indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
        DO i=1, Nindeppolys(it+1)
          PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':',                                   &
            polysIndep(i,1:NpolysIndep(i,it+1),it+1)
        END DO
      END IF
    END DO

    IF (dbg) THEN
      PRINT *,  'Coincidencies to connect _______'
      DO it=1, dt
        PRINT *,'  it:', it, ' Nindep:', Nindeppolys(it)
        PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs'
        DO ip=1, Nindeppolys(it)
          PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':',                         &
            polysIndep(ip,1:NpolysIndep(ip,it),it)
        END DO
      END DO

    END IF

    ! Trajectories
    ! It should be done following the number of 'independent' polygons
    ! One would concatenate that independent polygons which share IDs from one step to another

    ! First time-step. Take all polygons
    itrack = 0
    tracks = 0.
    Ntracks = 0
    DO ip=1, Nindeppolys(1)
      itrack = itrack + 1
      tracks(1,1,itrack,1) = itrack*1.
      tracks(2,1,itrack,1) = SpolysIndep(ip,1)
      tracks(3,1,itrack,1) = ctrpolys(1,ip,1)
      tracks(4,1,itrack,1) = ctrpolys(2,ip,1)
      tracks(5,1,itrack,1) = 1
      Ntracks(1) = Ntracks(1) + 1
    END DO

    ! Looping allover already assigned tracks
    timesteps: DO it=2, dt
      IF (dbg) PRINT *,'timestep:', it, 'N indep polys:', Nindeppolys(it)
      ! Indep polygons current time-step
      current_poly: DO i=1, Nindeppolys(it)
        IF (dbg) PRINT *,'  curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:',   &
          NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it)
        Indeppolychained = .FALSE.

        ! Number of tracks previous time-step
        ! Looping overall
        it1_tracks: DO itt=1, Ntracks(it-1)
          itrack = tracks(1,1,itt,it-1)
          ! Number polygons ID assigned
          Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0)
          IF (dbg) PRINT *,itt,'  track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1)

          ! Looking for coincidencies
          DO iip=1, Ntprev
            IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN
              Indeppolychained = .TRUE.
              IF (dbg) PRINT *,'    coincidence found by polygon:', tracks(2,iip,itt,it-1)
              EXIT
            END IF
          END DO
          IF (Indeppolychained) THEN
            Ntracks(it) = Ntracks(it) + 1
            ictrack = Ntracks(it)
            ! Assigning all the IDs to the next step of the track
            DO iip=1, NpolysIndep(i,it)
              iiprev = polysIndep(i,iip,it)
              tracks(1,iip,ictrack,it) = itrack
              tracks(2,iip,ictrack,it) = iiprev
              ixp = ctrpolys(1,iiprev,it)
              iyp = ctrpolys(2,iiprev,it)
              tracks(3,iip,ictrack,it) = ixp
              tracks(4,iip,ictrack,it) = iyp
              tracks(5,iip,ictrack,it) = it
            END DO
            EXIT
          END IF
        END DO it1_tracks

        ! Creation of a new track
        IF (.NOT.Indeppolychained) THEN
          Ntracks(it) = Ntracks(it) + 1
          ictrack = Ntracks(it)
          ! ID of new track
          maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.))
          IF (dbg) PRINT *,'  New track!', maxtrack+1

          ! Assigning all the IDs to the next step of the track
          DO j=1, NpolysIndep(i,it)
            iiprev = polysIndep(i,j,it)
            tracks(1,j,ictrack,it) = maxtrack+1
            tracks(2,j,ictrack,it) = iiprev
            ixp = ctrpolys(1,iiprev,it)
            iyp = ctrpolys(2,iiprev,it)
            tracks(3,j,ictrack,it) = ixp
            tracks(4,j,ictrack,it) = iyp
            tracks(5,j,ictrack,it) = it
          END DO
        END IF

      END DO current_poly

      IF (dbg) THEN
        PRINT *,'  At this time-step:', it, ' N trajectories:', Ntracks(it)
        DO i=1, Ntracks(it)
          Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0)
          PRINT *,i,'tracks:', tracks(2,1:Nprev,i,it)
        END DO
      END IF

    END DO timesteps

    ! Summarizing trajectories
    ! When multiple polygons are available, the mean of their central positions determines the position

    finaltracks = 0.
    maxtrack = MAXVAL(tracks(1,:,:,:))

    DO it=1, dt
      DO itt=1, Ntracks(it)
        itrack = INT(tracks(1,1,itt,it))
        Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0)
        PRINT *,'it:', it,'itrack:', itrack, 'Nprev:', Nprev
        finaltracks(1,itrack,it) = itrack*1.
        finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev
        finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev
        finaltracks(4,itrack,it) = it*1.
        PRINT *,'  finaltrack:', finaltracks(:,itrack,it)
      END DO
    END DO

    RETURN

  END SUBROUTINE poly_overlap_tracks

  SUBROUTINE coincidence_all_polys(dbg, dx, dy, Nallpoly, allpoly, icpolys, Npoly, polys, cpolys,     &
    apolys, polycoins, coinNptss)
! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
!   In case of multiple coincidencies, the closest and then the biggest is taken

    IMPLICIT NONE

    LOGICAL, INTENT(in)                                  :: dbg
    INTEGER, INTENT(in)                                  :: dx, dy, Nallpoly, Npoly
    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: allpoly, polys
    REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in)         :: icpolys
    REAL(r_k), DIMENSION(2,Npoly), INTENT(in)            :: cpolys
    REAL(r_k), DIMENSION(Npoly), INTENT(in)              :: apolys
    INTEGER, DIMENSION(Nallpoly), INTENT(out)            :: polycoins
    INTEGER, DIMENSION(Nallpoly,Npoly), INTENT(out)      :: coinNptss

! Local
    INTEGER                                              :: i, j, ip
    INTEGER                                              :: maxcorr
    INTEGER                                              :: Nmaxcorr
    LOGICAL, DIMENSION(dx,dy)                            :: boolpoly
    INTEGER                                              :: maxcoin
    REAL                                                 :: dist, maxcoindist, maxcoinarea

!!!!!!! Variables
! dx,dy: dimension of the space
! Nallpoly: Number of polygons to find coincidence
! allpoly: space with the polygons to meet
! icpolys: center of the polygons to look for the coincidence
! Npoly: number of polygons on the 2D space
! polys: 2D field of polygons identified by their integer number (0 for no polygon)
! cpolys: center of the polygons
! apolys: area of the polygons
! polycoins: coincident polyogn
!          -1: no-coincidence
!   1 < Npoly: single coincidence with a given polygon
!          -9: coincidence with more than one polygon
! coinNptss: number of points coincident with each polygon

    fname = 'coincidence_all_polys'
    IF (dbg) PRINT *,TRIM(fname)

    DO ip=1, Nallpoly
      boolpoly = allpoly == ip
      CALL coincidence_poly(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), coinNptss(ip,:))
      IF (dbg) PRINT *,'  polygon', ip, ' coincidence with:', polycoins(ip)

      ! Coincidence with more than one polygon
      IF (polycoins(ip) == -9) THEN
        maxcoindist = -10.
        maxcoinarea = -10.
        maxcoin = MAXVAL(coinNptss(ip,:))
        DO j=1, Npoly
          IF (coinNptss(ip,j) == maxcoin) THEN
            dist = SQRT( (icpolys(1,ip)*1.-cpolys(1,j)*1.)**2 + (icpolys(2,ip)*1.-cpolys(2,j)*1.)**2 )
            IF ( dist > maxcoindist) THEN
              maxcoindist = dist
              maxcoinarea = apolys(j)
              polycoins(ip) = j
            ELSE IF ( dist == maxcoindist) THEN
              IF (apolys(j) > maxcoinarea) THEN
                polycoins(ip) = j
                maxcoinarea = apolys(j)
              END IF
            END IF 
          END IF
        END DO
      END IF
    END DO

    RETURN

  END SUBROUTINE coincidence_all_polys

  SUBROUTINE coincidence_poly(dbg, dx, dy, poly, Npoly, polys, polycoin, coinNpts)
! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons

    IMPLICIT NONE

    LOGICAL, INTENT(in)                                  :: dbg
    INTEGER, INTENT(in)                                  :: dx, dy, Npoly
    LOGICAL, DIMENSION(dx,dy), INTENT(in)                :: poly
    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: polys
    INTEGER, INTENT(out)                                 :: polycoin
    INTEGER, DIMENSION(Npoly), INTENT(out)               :: coinNpts

! Local
    INTEGER                                              :: i, j, ip
    INTEGER                                              :: maxcorr
    INTEGER                                              :: Nmaxcorr

!!!!!!! Variables
! dx,dy: dimension of the space
! poly: bolean polygon to meet
! Npoly: number of polygons on the 2D space
! polys: 2D field of polygons identified by their integer number (0 for no polygon)
! polycoin: coincident polyogn
!          -1: no-coincidence
!   1 < Npoly: single coincidence with a given polygon
!          -9: coincidence with more than one polygon
! coinNpts: number of points coincident with each polygon

    fname = 'coincidence_poly'
    IF (dbg) PRINT *,TRIM(fname)

    IF (dbg) THEN
      PRINT *,'  Boolean polygon to search coincidences ...'
      DO i=1,dx
        PRINT *,poly(i,:)
      END DO

      PRINT *,'  2D polygons space ...'
      DO i=1,dx
        PRINT '(1000(I3,1x))',polys(i,:)
      END DO
    END IF

    ! Looking for coincient points for the polygon
    coinNpts = 0
    DO i=1,dx
      DO j=1,dy
        IF (poly(i,j) .AND. polys(i,j) .NE. 0) coinNpts(polys(i,j)) = coinNpts(polys(i,j)) + 1
      END DO
    END DO

    maxcorr = 0
    maxcorr = INT(MAXVAL(coinNpts*1.))

    IF (dbg) PRINT *,'  Maximum coincidence:', maxcorr
    IF (maxcorr == 0) THEN
      polycoin = -1
    ELSE
      Nmaxcorr = 0
      DO ip=1, Npoly
        IF (coinNpts(ip) == maxcorr) THEN
          Nmaxcorr=Nmaxcorr+1
          polycoin = ip
        END IF
      END DO
      IF (Nmaxcorr > 1) polycoin = -9
    END IF

    IF (dbg) THEN
      PRINT *,'  Coincidences for each polygon _______'
      DO ip=1, Npoly
        PRINT *,'  ', ip,': ', coinNpts(ip)
      END DO
    END IF

    RETURN

END SUBROUTINE coincidence_poly

  SUBROUTINE all_polygons_properties(dbg, dx, dy, Npoly, polys, lon, lat, values, xres, yres, projN,  &
    Npolyptss, xxtrms, yxtrms, meanctrs, meanwctrs, areas, nvals, xvals, mvals, m2vals, stdvals,      &
    Nquant, quants, nvcoords, xvcoords, meanvnctrs, meanvxctrs)
! Subroutine to determine the properties of all polygons in a 2D field:
!   Number of grid points
!   grid-point coordinates of the minimum and maximum of the path along x,y axes
!   grid coordinates of center from the mean of the coordinates of the poygon locations
!   lon, lat center from the area weighted mean of the coordinates of the polygon locations
!   area of the polygon (km2)
!   minimum and maximum of the values within the polygon
!   mean of the values within the polygon
!   quadratic mean of the values within the polygon
!   standard deviation of the values within the polygon
!   number of quantiles
!   quantiles of the values within the polygon
!   grid coordinates of the minimum, maximum value within the polygon
!   lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon

  IMPLICIT NONE

  LOGICAL, INTENT(in)                                    :: dbg
  INTEGER, INTENT(in)                                    :: dx, dy, Npoly, Nquant
  INTEGER, DIMENSION(dx,dy), INTENT(in)                  :: polys
  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: lon, lat, values
  REAL(r_k), INTENT(in)                                  :: xres, yres
  CHARACTER(len=1000), INTENT(in)                        :: projN
  INTEGER, DIMENSION(Npoly), INTENT(out)                 :: Npolyptss
  INTEGER, DIMENSION(Npoly,2), INTENT(out)               :: xxtrms, yxtrms, meanctrs
  REAL(r_k), DIMENSION(Npoly), INTENT(out)               :: areas
  REAL(r_k), DIMENSION(Npoly), INTENT(out)               :: nvals, xvals, mvals, m2vals, stdvals
  REAL(r_k), DIMENSION(Npoly, Nquant), INTENT(out)       :: quants
  INTEGER, DIMENSION(Npoly,2), INTENT(out)               :: nvcoords, xvcoords
  REAL(r_k), DIMENSION(Npoly,2), INTENT(out)             :: meanwctrs, meanvnctrs, meanvxctrs

! Local
  INTEGER                                                :: ip
  LOGICAL, DIMENSION(dx,dy)                              :: boolpoly

!!!!!!! Variables
! dx,dy: size of the space
! Npoly: number of polygons
! polys: polygon matrix with all polygons (as integer number per polygon)
! lon, lat: geographical coordinates of the grid-points of the matrix
! values: values of the 2D field to use
! [x/y]res resolution along the x and y axis
! projN: name of the projection
!   'lon/lat': for regular longitude-latitude
!   'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
! Npolyptss: number of points of the polygons
! [x/y]xtrms: grid-point coordinates of minimum and maximum coordinates of the polygons
! meanctrs: center from the mean of the coordinates of the polygons
! meanwctrs: lon, lat coordinates of the center from the spatial-weighted mean of the polygons
! areas: area of the polygons [km]
! [n/x]vals: minimum and maximum of the values within the polygons
! mvals: mean of the values within the polygons
! m2vals: quadratic mean of the values within the polygons
! stdvals: standard deviation of the values within the polygons
! Nquant: number of quantiles
! quants: quantiles of the values within the polygons
! [n/x]vcoords: grid coordinates of the minimum/maximum of the values within the polygons
! meanv[n/x]ctrs: lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygons

  fname = 'all_polygons_properties'

  ! Initializing matrices
  Npolyptss = -1
  xxtrms = fillval64
  yxtrms = fillval64
  meanctrs = fillval64
  meanwctrs = fillval64
  areas = fillval64
  nvals = fillvalI
  xvals = fillval64
  mvals = fillval64
  m2vals = fillval64
  stdvals = fillval64
  quants = fillval64
  nvcoords = fillvalI
  xvcoords = fillvalI
  meanvnctrs = fillval64
  meanvxctrs = fillval64

  DO ip=1, Npoly
    boolpoly = polys == ip
    CALL polygon_properties(dbg, dx, dy, boolpoly, lon, lat, values, xres, yres, projN, Npolyptss(ip),&
      xxtrms(ip,:), yxtrms(ip,:), meanctrs(ip,:), meanwctrs(ip,:), areas(ip), nvals(ip), xvals(ip),   &
      mvals(ip), m2vals(ip), stdvals(ip), Nquant, quants(ip,:), nvcoords(ip,:), xvcoords(ip,:),       &
      meanvnctrs(ip,:), meanvxctrs(ip,:))
  END DO

  RETURN

  END SUBROUTINE all_polygons_properties

  SUBROUTINE polygon_properties(dbg, dx, dy, poly, lon, lat, values, xres, yres, projN, Npolypts,     &
    xxtrm, yxtrm, meanctr, meanwctr, area, nval, xval, mval, m2val, stdval, Nquant, quant, nvcoord,   &
    xvcoord, meanvnctr, meanvxctr)
! Subroutine to determine the properties of a polygon (as .TRUE. matrix)
!   Number of grid points
!   grid-point coordinates of the minimum and maximum of the path along x,y axes
!   grid coordinates of center from the mean of the coordinates of the poygon locations
!   lon, lat center from the area weighted mean of the coordinates of the polygon locations
!   area of the polygon (km2)
!   minimum and maximum of the values within the polygon
!   mean of the values within the polygon
!   quadratic mean of the values within the polygon
!   standard deviation of the values within the polygon
!   number of quantiles
!   quantiles of the values within the polygon
!   grid coordinates of the minimum, maximum value within the polygon
!   lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon

  IMPLICIT NONE

  LOGICAL, INTENT(in)                                    :: dbg
  INTEGER, INTENT(in)                                    :: dx, dy, Nquant
  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: poly
  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: lon, lat, values
  REAL(r_k), INTENT(in)                                  :: xres, yres
  CHARACTER(len=1000), INTENT(in)                        :: projN
  INTEGER, INTENT(out)                                   :: Npolypts
  INTEGER, DIMENSION(2), INTENT(out)                     :: xxtrm, yxtrm, meanctr
  INTEGER, DIMENSION(2), INTENT(out)                     :: nvcoord, xvcoord
  REAL(r_k), DIMENSION(2), INTENT(out)                   :: meanwctr, meanvnctr, meanvxctr
  REAL(r_k), INTENT(out)                                 :: area
  REAL(r_k), INTENT(out)                                 :: nval, xval, mval, m2val, stdval
  REAL(r_k), DIMENSION(Nquant), INTENT(out)              :: quant

! Local
  INTEGER                                                :: i, j, ip
  INTEGER                                                :: ierr
  INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: polypts
  REAL(r_k), DIMENSION(:), ALLOCATABLE                   :: polyvals, distvn, distvx
  REAL(r_k)                                              :: lonNADIR, latNADIR
  REAL(r_k)                                              :: sumRESx, sumRESy
  REAL(r_k), DIMENSION(dx,dy)                            :: xcorr, ycorr
  CHARACTER(len=200), DIMENSION(3)                       :: satSvals
  CHARACTER(len=50)                                      :: projS
  REAL(r_k)                                              :: sumDISTnlon, sumDISTnlat, sumDISTxlon,   &
    sumDISTxlat

!!!!!!! Variables
! dx,dy: size of the space
! poly: polygon matrix (boolean)
! lon, lat: geographical coordinates of the grid-points of the matrix
! values: values of the 2D field to use
! [x/y]res resolution along the x and y axis
! projN: name of the projection
!   'ctsarea': Constant Area
!   'lon/lat': for regular longitude-latitude
!   'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
! Npolypts: number of points of the polygon
! [x/y]xtrm: grid-point coordinates of minimum and maximum coordinates of the polygon
! meanctr: center from the mean of the coordinates of the polygon
! meanwctr: lon, lat coordinates of the center from the spatial-weighted mean of the polygon
! area: area of the polygon [km]
! [n/x]val: minimum and maximum of the values within the polygon
! mval: mean of the values within the polygon
! m2val: quadratic mean of the values within the polygon
! stdval: standard deviation of the values within the polygon
! Nquant: number of quantiles
! quant: quantiles of the values within the polygon 
! [n/x]vcoord: grid coordinates of the minimum/maximum of the values within the polygon
! meanv[n/x]ctr: lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon

  fname = 'polygon_properties'

  IF (dbg) PRINT *,"  '" // TRIM(fname) // "' ..."

  ! Getting grid-point coordinates of the polygon
  Npolypts = COUNT(poly)

  IF (ALLOCATED(polypts)) DEALLOCATE(polypts)
  ALLOCATE(polypts(Npolypts,2), STAT=ierr)
  msg = "Problems allocating 'polypts'"
  CALL ErrMsg(msg, fname, ierr)

  IF (ALLOCATED(polyvals)) DEALLOCATE(polyvals)
  ALLOCATE(polyvals(Npolypts), STAT=ierr)
  msg = "Problems allocating 'polyvals'"
  CALL ErrMsg(msg, fname, ierr)

  IF (ALLOCATED(distvn)) DEALLOCATE(distvn)
  ALLOCATE(distvn(Npolypts), STAT=ierr)
  msg = "Problems allocating 'distvn'"
  CALL ErrMsg(msg, fname, ierr)

  IF (ALLOCATED(distvx)) DEALLOCATE(distvx)
  ALLOCATE(distvx(Npolypts), STAT=ierr)
  msg = "Problems allocating 'distvx'"
  CALL ErrMsg(msg, fname, ierr)

  IF (projN(1:7) == 'lon/lat') THEN
    projS = projN
  ELSE IF (projN(1:7) == 'ctsarea') THEN
    projS = projN
  ELSE IF (projN(1:9) == 'nadir-sat') THEN
    projS = 'nadir-sat'
    CALL split(projN, ',', 3, satSvals)
    READ(satSvals(2),'(F200.100)')lonNadir
    READ(satSvals(3),'(F200.100)')latNadir
    IF (dbg) PRINT *,"  'nadir-geostationary-satellite' based projection of data with nadir " //      &
      "location at:", lonNadir, latNadir
  ELSE
    msg = "Projection '" // TRIM(projN) // "' not ready" // CHAR(10) // " available ones: " //        &
       "'ctsarea', 'lon/lat', 'nadir-sat'"
    CALL ErrMsg(msg,fname,-1)
  END IF

  area = 0.
  sumRESx = 0.
  sumRESy = 0.
  meanwctr = 0.
  meanvnctr = 0.
  meanvxctr = 0.

  nval = fillval64
  xval = -fillval64

  ip = 1
  DO i=1,dx
    DO j=1,dy
      IF (poly(i,j)) THEN
        polypts(ip,1) = i
        polypts(ip,2) = j
        polyvals(ip) = values(i,j)
        SELECT CASE (TRIM(projS))
          CASE ('ctsarea')
            ! Constant Area
            xcorr(i,j) = 1.
            ycorr(i,j) = 1.
          CASE ('lon/lat')
            ! Area as fixed yres and sinus-lat varying for xres
            xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
            ycorr(i,j) = 1.
          CASE ('nadir-sat')
            ! Area from nadir resolution and degrading as we get far from satellite's nadir
            ! GOES-E: 0 N, 75 W
            xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
            ycorr(i,j) = 1.
        END SELECT
        area = area + xres*xcorr(i,j)*yres*ycorr(i,j)
        meanwctr(1) = meanwctr(1) + lon(i,j)*xres*xcorr(i,j)
        meanwctr(2) = meanwctr(2) + lat(i,j)*yres*ycorr(i,j)
        IF (nval > values(i,j)) THEN
          nvcoord(1) = i
          nvcoord(2) = j
          nval = values(i,j)
        END IF
        IF (xval < values(i,j)) THEN
          xvcoord(1) = i
          xvcoord(2) = j
          xval = values(i,j)
        END IF
        ip = ip + 1
      END IF
    END DO
  END DO

  IF (dbg) THEN
    PRINT *,'  grdi_coord lon lat value _______ '
    DO ip=1, Npolypts
      PRINT *, polypts(ip,:), ';', lon(polypts(ip,1),polypts(ip,2)), lat(polypts(ip,1),polypts(ip,2)),&
         ':', polyvals(ip)
    END DO
  END IF

  sumRESx = xres*SUM(xcorr)
  sumRESy = yres*SUM(ycorr)

  xxtrm = (/ MINVAL(polypts(:,1)), MAXVAL(polypts(:,1)) /)
  yxtrm = (/ MINVAL(polypts(:,2)), MAXVAL(polypts(:,2)) /)
  meanctr = (/ SUM(polypts(:,1))/Npolypts, SUM(polypts(:,2))/Npolypts /)
  meanwctr = (/ meanwctr(1)/sumRESx, meanwctr(2)/sumRESy /)

  IF (dbg) THEN
    PRINT *,'  mean grid center: ', meanctr, ' weighted mean center: ', meanwctr
  END IF
  
  ! Statistics of the values within the polygon
  CALL StatsR_K(Npolypts, polyvals, nval, xval, mval, m2val, stdval)

  IF (dbg) THEN
    PRINT *,'    minimum value: ', nval, ' maximum:', xval, ' mean:', mval
    PRINT *,'    coor. minimum: ', nvcoord
    PRINT *,'    coor. maximum: ', xvcoord
  END IF

  ! Mean center weighted to minimum and maximum values
  distvn = DABS(polyvals - nval)
  distvx = DABS(polyvals - xval)

  meanvnctr = 0.
  meanvxctr = 0.
  sumDISTnlon = 0.
  sumDISTnlat = 0.
  sumDISTxlon = 0.
  sumDISTxlat = 0.

  ip = 1
  DO i=1,dx
    DO j=1,dy
      IF (poly(i,j)) THEN
        meanvnctr(1) = meanvnctr(1)+lon(i,j)*distvn(ip)*xres*xcorr(i,j)
        meanvnctr(2) = meanvnctr(2)+lat(i,j)*distvn(ip)*yres*ycorr(i,j)

        meanvxctr(1) = meanvxctr(1)+lon(i,j)*distvx(ip)*xres*xcorr(i,j)      
        meanvxctr(2) = meanvxctr(2)+lat(i,j)*distvx(ip)*yres*ycorr(i,j)

        sumDISTnlon = sumDISTnlon + distvn(ip)*xres*xcorr(i,j)
        sumDISTnlat = sumDISTnlat + distvn(ip)*yres*ycorr(i,j)
        sumDISTxlon = sumDISTxlon + distvx(ip)*xres*xcorr(i,j)
        sumDISTxlat = sumDISTxlat + distvx(ip)*yres*ycorr(i,j)

        ip = ip + 1
      END IF
    END DO
  END DO

  meanvnctr = (/ meanvnctr(1)/sumDISTnlon, meanvnctr(2)/sumDISTnlat /)
  meanvxctr = (/ meanvxctr(1)/sumDISTxlon, meanvxctr(2)/sumDISTxlat /)

  IF (dbg) THEN
    PRINT *,'  mean center to minimum: ', meanvnctr, ' to maximum: ', meanvxctr
  END IF

  ! Quantiles of the values within the polygon
  quant = -9999.d0
  IF (Npolypts > Nquant) THEN
    CALL quantilesR_K(Npolypts, polyvals, Nquant, quant)
  ELSE
    CALL SortR_K(polyvals, Npolypts)
  END IF

  DEALLOCATE (polypts)
  DEALLOCATE (polyvals)

  RETURN

  END SUBROUTINE polygon_properties

SUBROUTINE polygons_t(dbg, dx, dy, dt, boolmatt, polys, Npoly)
! Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: dx, dy, dt
  LOGICAL, DIMENSION(dx,dy,dt), INTENT(in)               :: boolmatt
  LOGICAL, INTENT(in)                                    :: dbg
  INTEGER, DIMENSION(dt), INTENT(out)                    :: Npoly
  INTEGER, DIMENSION(dx,dy,dt), INTENT(out)              :: polys

! Local
  INTEGER                                                :: i,it

!!!!!!! Variables
! dx,dy: spatial dimensions of the space
! boolmatt: boolean matrix tolook for the polygons (.TRUE. based)
! polys: found polygons
! Npoly: number of polygons found

  fname = 'polygons'

  IF (dbg) PRINT *,TRIM(fname)

  polys = -1
  Npoly = 0

  DO it=1,dt
    IF (ANY(boolmatt(:,:,it))) THEN
      IF (dbg) THEN
        PRINT *,'  it:', it, ' num. TRUE:', COUNT(boolmatt(:,:,it)), 'bool _______'
        DO i=1,dx
          PRINT *,boolmatt(i,:,it)
        END DO
      END IF
      CALL polygons(dbg, dx, dy, boolmatt(:,:,it), polys(:,:,it), Npoly(it))
    ELSE
      IF (dbg) THEN
        PRINT *,'  it:', it, " without '.TRUE.' values skipiing it!!"
      END IF
    END IF
  END DO

END SUBROUTINE polygons_t

SUBROUTINE polygons(dbg, dx, dy, boolmat, polys, Npoly)
! Subroutine to search the polygons of a boolean field. FORTRAN based. 1st = 1!

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: dx, dy
  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: boolmat
  LOGICAL, INTENT(in)                                    :: dbg
  INTEGER, INTENT(out)                                   :: Npoly
  INTEGER, DIMENSION(dx,dy), INTENT(out)                 :: polys

! Local
  INTEGER                                                :: i, j, ip, ipp, Nppt
  INTEGER                                                :: ierr
  INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: borders
  LOGICAL, DIMENSION(dx,dy)                              :: isborder, isbordery
  INTEGER, DIMENSION(:,:,:), ALLOCATABLE                 :: paths
  INTEGER                                                :: Npath
  INTEGER, DIMENSION(:), ALLOCATABLE                     :: Nptpaths
  INTEGER, DIMENSION(2)                                  :: xtrx, xtry, meanpth
  INTEGER                                                :: Nvertx, Npts
  INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: vertxs, points
  LOGICAL, DIMENSION(:), ALLOCATABLE                     :: isin

!!!!!!! Variables
! dx,dy: spatial dimensions of the space
! boolmat: boolean matrix tolook for the polygons (.TRUE. based)
! polys: found polygons
! Npoly: number of polygons found

  fname = 'polygons'

  polys = -1

  ! The mathematical maximum woiuld be dx*dy/4, but let's be optimistic... (sorry Jero)
  Nppt = dx*dy/10

  IF (ALLOCATED(borders)) DEALLOCATE(borders)
  ALLOCATE(borders(Nppt,2), STAT=ierr)
  msg = "Problems allocating matrix 'borders'"
  CALL ErrMsg(msg, fname, ierr)

  IF (ALLOCATED(paths)) DEALLOCATE(paths)
  ALLOCATE(paths(Nppt,Nppt,2), STAT=ierr)
  msg = "Problems allocating matrix 'paths'"
  CALL ErrMsg(msg, fname, ierr)

  IF (ALLOCATED(Nptpaths)) DEALLOCATE(Nptpaths)
  ALLOCATE(Nptpaths(Nppt), STAT=ierr)
  msg = "Problems allocating matrix 'Nptpaths'"
  CALL ErrMsg(msg, fname, ierr)

  ! Filling with the points of all the space with .TRUE.
  Npts = COUNT(boolmat)

  IF (ALLOCATED(points)) DEALLOCATE(points)
  ALLOCATE(points(Npts,2), STAT=ierr)
  msg = "Problems allocating matrix 'points'"
  CALL ErrMsg(msg, fname, ierr)

  ! We only want to localize that points 'inside'
  ip = 1
  DO i=1, dx
    DO j=1, dy
      IF (boolmat(i,j)) THEN
        points(ip,1) = i
        points(ip,2) = j
        ip = ip + 1
      END IF
    END DO
  END DO

  CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery)
  CALL paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)

  Npoly = Npath

  DO ip=1, Npath
    IF (ALLOCATED(vertxs)) DEALLOCATE(vertxs)
    ALLOCATE(vertxs(Nptpaths(ip),2))
    msg = "Problems allocating matrix 'vertxs'"
    CALL ErrMsg(msg, fname, ierr)

    IF (ALLOCATED(isin)) DEALLOCATE(isin)
    ALLOCATE(isin(Npts), STAT=ierr)
    msg = "Problems allocating matrix 'isin'"
    CALL ErrMsg(msg, fname, ierr)

    isin = .FALSE.

    IF (dbg) PRINT *, '  path:', ip, ' N pts:', Nptpaths(ip)

    CALL path_properties(dx, dy, boolmat, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), xtrx, xtry,       &
      meanpth, 'y', Nvertx, vertxs)

    IF (dbg) THEN
      PRINT *, '    properties  _______'
      PRINT *, '    x-extremes:', xtrx
      PRINT *, '    y-extremes:', xtry
      PRINT *, '    center mean:', meanpth
      PRINT *, '    y-vertexs:', Nvertx,' ________'
      DO i=1, Nvertx
        PRINT *,'      ',i,':',vertxs(i,:)
      END DO
    END IF
 
    CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,&
      xtrx, xtry, vertxs, Npts, points, isin)

    ! Filling polygons
    DO ipp=1, Npts
      IF (isin(ipp)) polys(points(ipp,1),points(ipp,2)) = ip
    END DO

    IF (dbg) THEN
      PRINT *,'  boolmat isborder isbordery polygon (',xtrx(1),',',xtry(1),')x(',xtrx(2),',',xtry(2), &
        ') _______' 
      DO i=xtrx(1), xtrx(2)
        PRINT *,i,':',boolmat(i,xtry(1):xtry(2)), ' border ', isborder(i,xtry(1):xtry(2)),            &
          ' isbordery ', isbordery(i,xtry(1):xtry(2)), ' polygon ', polys(i,xtry(1):xtry(2))
      END DO
    END IF

  END DO

  ! Cleaning polygons matrix of not-used and paths around holes
  CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg)

  DEALLOCATE (borders)  
  DEALLOCATE (Nptpaths)
  DEALLOCATE (paths)
  DEALLOCATE (vertxs)
  DEALLOCATE (points)
  DEALLOCATE (isin)

  RETURN

END SUBROUTINE polygons

SUBROUTINE clean_polygons(dx, dy, Lmat, pols, Npols, dbg)
! Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: dx, dy
  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: Lmat
  INTEGER, INTENT(inout)                                 :: Npols
  INTEGER, DIMENSION(dx,dy), INTENT(inout)               :: pols
  LOGICAL, INTENT(in)                                    :: dbg

! Local
  INTEGER                                                :: i,j,ip,iprm
  INTEGER, DIMENSION(Npols)                              :: origPol, NotPol, neigPol
  INTEGER                                                :: ispol, NnotPol
  CHARACTER(len=4)                                       :: ISa

!!!!!!! Variables
! dx, dy: size of the space
! Lmat: original bolean matrix from which the polygons come from
! Npols: original number of polygons
! pols: polygons space

  fname = 'clean_polygons'
  IF (dbg) PRINT *,"  At '" // TRIM(fname) // "' ..."

  origPol = -1

  ! Looking for polygons already in space
  NnotPol = 0
  DO ip=1, Npols
    ispol = COUNT(pols-ip == 0)
    IF (ispol > 0) THEN
      origPol(ip) = ip
    ELSE
      NnotPol = NnotPol + 1
      NotPol(NnotPol) = ip
      neigPol(NnotPol) = -1
    END IF
  END DO

  IF (dbg) THEN
    PRINT *,'  It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol
    PRINT *,'  Polygons to remove:', NotPol(1:NnotPol)
  END IF
  
  ! Looking for the hole border of a polygon. This is identify as such polygon point which along 
  !   y-axis has NpolygonA, Npolygon, .FALSE.
  DO i=1,dx
    DO j=2,dy-1
      IF  ( (pols(i,j-1) /= pols(i,j) .AND. pols(i,j+1) == -1) .AND. (COUNT(NotPol-pols(i,j)==0)==0)  &
        .AND. (pols(i,j) /= -1) .AND. (pols(i,j-1) /= -1)) THEN
        IF (dbg) PRINT *,'  Polygon:', pols(i,j), ' to be removed at point (',i,',',j,'); j-1:',      &
          pols(i,j-1), ' j:', pols(i,j), ' j+1:', pols(i,j+1)
        NnotPol = NnotPol + 1
        NotPol(NnotPol) = pols(i,j)
        neigPol(NnotPol) = pols(i,j-1)
      END IF
    END DO
  END DO

  IF (dbg) THEN
    PRINT *,'  It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol
    PRINT *,'  Polygons to remove after looking for fake border-of-hole polygons _______'
    DO i=1, NnotPol
      PRINT *, '      Polygon:', NotPol(i), ' to be replaced by:', neigPol(i)
    END DO
  END IF

  ! Removing polygons
  DO iprm=1, NnotPol
    IF (neigPol(iprm) == -1) THEN
      WHERE (pols == NotPol(iprm))
        pols = -1
      END WHERE
      IF (dbg) THEN
        PRINT *,'    removing polygon:', NotPol(iprm)
      END IF
    ELSE
      WHERE (pols == NotPol(iprm))
        pols = neigPol(iprm)
      END WHERE
      IF (dbg) THEN
        PRINT *,'       replacing polygon:', NotPol(iprm), ' by:', neigPol(iprm)
      END IF
    END IF
  END DO

  ! Re-numbering (descending values)
  DO i = 1, NnotPol
    iprm = MAXVAL(NotPol(1:NnotPol))
    WHERE(pols > iprm)
      pols = pols - 1
    END WHERE
    j = Index1DArrayI(NotPol, NnotPol, iprm)
    NotPol(j) = -9
  END DO

  Npols = Npols - NnotPol

  RETURN

END SUBROUTINE clean_polygons

  SUBROUTINE path_properties(dx, dy, Lmat, Nptspth, pth, xxtrm, yxtrm, meanctr, axs, Nvrtx, vrtxs)
! Subroutine to determine the properties of a path:
!   extremes: minimum and maximum of the path along x,y axes
!   meancenter: center from the mean of the coordinates of the paths locations
!   vertexs: path point, without neighbours along a given axis

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: dx, dy, Nptspth
  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: Lmat
  INTEGER, DIMENSION(Nptspth,2), INTENT(in)              :: pth
  CHARACTER, INTENT(in)                                  :: axs
  INTEGER, DIMENSION(2), INTENT(out)                     :: meanctr, xxtrm, yxtrm
  INTEGER, INTENT(out)                                   :: Nvrtx
  INTEGER, DIMENSION(Nptspth,2), INTENT(out)             :: vrtxs

! Local
  INTEGER                                                :: i, ip, jp
  INTEGER                                                :: neig1, neig2

!!!!!!! Variables
! dx,dy: size of the space
! Lmat: original matrix of logical values for the path
! Nptspth: number of points of the path
! pth: path coordinates (clockwise)
! axs: axis of finding the vertex
! [x/y]xtrm: minimum and maximum coordinates of the path
! meanctr: center from the mean of the coordinates of the path
! Nvrtx: Number of vertexs of the path
! vrtxs: coordinates of the vertexs

  fname = 'path_properties'

  vrtxs = -1
  Nvrtx = 0

  xxtrm = (/ MINVAL(pth(:,1)), MAXVAL(pth(:,1)) /)
  yxtrm = (/ MINVAL(pth(:,2)), MAXVAL(pth(:,2)) /)
  meanctr = (/ SUM(pth(:,1))/Nptspth, SUM(pth(:,2))/Nptspth /)

  IF (axs == 'x' .OR. axs == 'X') THEN
    ! Looking vertexs along x-axis
    DO i=1, Nptspth
      ip = pth(i,1)
      jp = pth(i,2)
      neig1 = 0
      neig2 = 0
      ! W-point 
      IF (ip == 1) THEN
        neig1 = -1
      ELSE
        IF (.NOT.Lmat(ip-1,jp)) neig1 = -1
      END IF
      ! E-point 
      IF (ip == dx) THEN
        neig2 = -1
      ELSE
        IF (.NOT.Lmat(ip+1,jp)) neig2 = -1
      END IF
    
      IF (neig1 == -1 .AND. neig2 == -1) THEN
        Nvrtx = Nvrtx + 1
        vrtxs(Nvrtx,:) = (/ip,jp/)
      END IF
    END DO
  ELSE IF (axs == 'y' .OR. axs == 'Y') THEN
    ! Looking vertexs along x-axis
    DO i=1, Nptspth
      ip = pth(i,1)
      jp = pth(i,2)

      neig1 = 0
      neig2 = 0
      ! S-point 
      IF (jp == 1) THEN
        neig1 = -1
      ELSE
        IF (.NOT.Lmat(ip,jp-1)) neig1 = -1
      END IF
      ! N-point 
      IF (jp == dy) THEN
        neig2 = -1
      ELSE
        IF (.NOT.Lmat(ip,jp+1)) neig2 = -1
      END IF

      IF (neig1 == -1 .AND. neig2 == -1) THEN
        Nvrtx = Nvrtx + 1
        vrtxs(Nvrtx,:) = (/ ip, jp /)
      END IF
    END DO
  ELSE
    msg = "Axis '" // axs // "' not available" // CHAR(10) // "  Available ones: 'x', 'X', 'y, 'Y'"
    CALL ErrMsg(msg, fname, -1)
  END IF

  RETURN

  END SUBROUTINE path_properties

  SUBROUTINE gridpoints_InsidePolygon(dx, dy, isbrdr, Npath, path, Nvrtx, xpathxtrm, ypathxtrm,       &
    vrtxs, Npts, pts, inside)
! Subroutine to determine if a series of grid points are inside a polygon following ray casting algorithm
! FROM: https://en.wikipedia.org/wiki/Point_in_polygon

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: dx,dy,Npath,Nvrtx,Npts
  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: isbrdr
  INTEGER, DIMENSION(Npath,2), INTENT(in)                :: path
  INTEGER, DIMENSION(2), INTENT(in)                      :: xpathxtrm, ypathxtrm
  INTEGER, DIMENSION(Npath,2)                            :: vrtxs
  INTEGER, DIMENSION(Npts,2), INTENT(in)                 :: pts
  LOGICAL, DIMENSION(Npts), INTENT(out)                  :: inside

! Local
  INTEGER                                                :: i,j,ip,ix,iy
  INTEGER                                                :: Nintersecs, isvertex, ispath
  INTEGER                                                :: ierr
  LOGICAL, DIMENSION(:,:), ALLOCATABLE                   :: halo_brdr
  INTEGER                                                :: Nbrbrdr

!!!!!!! Variables
! dx,dy: space size
! Npath: number of points of the path of the polygon
! path: path of the polygon
! isbrdr: boolean matrix of the space wqith .T. on polygon border
! Nvrtx: number of vertexs of the path
! [x/y]pathxtrm extremes of the path
! vrtxs: vertexs of the path along y-axis
! Npts: number of points
! pts: points to look for
! inside: vector wether point is inside or not (coincident to a border is inside)

  fname = 'gridpoints_InsidePolygon'

  ! Creation of a 1-grid point larger matrix to deal with points reaching the limits
  IF (ALLOCATED(halo_brdr)) DEALLOCATE(halo_brdr)
  ALLOCATE(halo_brdr(dx+2,dy+2), STAT=ierr)
  msg = "Problems allocating matrix 'halo_brdr'"
  CALL ErrMsg(msg, fname, ierr)
  halo_brdr = .FALSE.

  DO i=1,dx
    halo_brdr(i+1,2:dy+1) = isbrdr(i,:)
  END DO

  inside = .FALSE.

  DO ip=1,Npts
    Nintersecs = 0
    ix = pts(ip,1)
    iy = pts(ip,2)
    ! Point might be outside path range...
    IF (ix >= xpathxtrm(1) .AND. ix <= xpathxtrm(2) .AND. iy >= ypathxtrm(1) .AND.                    &
      iy <= ypathxtrm(2)) THEN

      ! It is a border point?
      ispath = index_list_coordsI(Npath, path, (/ix,iy/))
      IF (isbrdr(ix,iy) .AND. (ispath /= -1)) THEN
        inside(ip) = .TRUE.
        CYCLE
      END IF

      ! Looking along y-axis
      ! Accounting for consecutives borders
      Nbrbrdr = 0
      DO j=MAX(1,ypathxtrm(1)-1),iy-1
        ! Only counting that borders that are not vertexs
        ispath = index_list_coordsI(Npath, path, (/ix,j/))
        isvertex = index_list_coordsI(Npath, vrtxs, (/ix,j/))

        IF (halo_brdr(ix+1,j+1) .AND. (ispath /= -1) .AND. (isvertex == -1) ) Nintersecs = Nintersecs + 1 
        IF (halo_brdr(ix+1,j+1) .AND. (ispath /= -1) .AND. (halo_brdr(ix+1,j+1) .EQV. halo_brdr(ix+1,j+2))) THEN
          Nbrbrdr = Nbrbrdr + 1
        ELSE
          ! Will remove that consecutive borders above 2
          IF (Nbrbrdr /= 0) THEN
            Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0)
            Nbrbrdr = 0
          END IF
        END IF
      END DO
      IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE.
    END IF

  END DO

  RETURN

END SUBROUTINE gridpoints_InsidePolygon

SUBROUTINE look_clockwise_borders(dx,dy,Nbrdrs,brdrs,gbrdr,isbrdr,ix,iy,dbg,xf,yf,iff)
! Subroutine to look clock-wise for a next point within a collection of borders (limits of a region)

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: dx, dy, Nbrdrs, ix, iy
  INTEGER, DIMENSION(Nbrdrs,2), INTENT(in)               :: brdrs
  LOGICAL, DIMENSION(Nbrdrs), INTENT(in)                 :: gbrdr
  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: isbrdr
  LOGICAL, INTENT(in)                                    :: dbg
  INTEGER, INTENT(out)                                   :: xf, yf, iff

! Local
  INTEGER                                                :: isch
  CHARACTER(len=2), DIMENSION(8)                         :: Lclock
  INTEGER, DIMENSION(8,2)                                :: spt
  INTEGER                                                :: iif, jjf

!!!!!!! Variables
! dx, dy: 2D shape ot the space
! Nbrdrs: number of brdrs found in this 2D space
! brdrs: list of coordinates of the borders
! gbrdr: accounts for the use if the given border point
! isbrdr: accounts for the matrix of the point is a border or not
! ix,iy: coordinates of the point to start to find for
! xf,yf: coordinates of the found point
! iff: position of the border found within the list of borders

  fname = 'look_clockwise_borders'

  ! Looking clock-wise assuming that one starts from the westernmost point

  ! Label of the search
  lclock = (/ 'W ', 'NW', 'N ', 'NE', 'E ', 'SE', 'S ', 'SW' /)
  ! Transformation to apply
  !spt = (/ (/-1,0/), (/-1,1/), (/0,1/), (/1,1/), (/1,0/), (/1,-1/), (/0,-1/), (/-1,-1/) /)
  spt(:,1) = (/ -1, -1, 0, 1, 1, 1, 0, -1 /)
  spt(:,2) = (/ 0, 1, 1, 1, 0, -1, -1, -1 /)

  xf = -1
  yf = -1
  DO isch=1, 8
    ! clock-wise search
    IF (spt(isch,1) >= 0) THEN
      iif = MIN(dx,ix+spt(isch,1))
    ELSE
      iif = MAX(1,ix+spt(isch,1))
    END IF
    IF (spt(isch,2) >= 0) THEN
      jjf = MIN(dy,iy+spt(isch,2))
    ELSE
      jjf = MAX(1,iy+spt(isch,2))
    END IF
    iff = index_list_coordsI(Nbrdrs, brdrs,(/iif,jjf/))
    IF (iff > 0) THEN
      IF (dbg) PRINT *,'    ' // lclock(isch) // '-point:', iif,jjf, ':', iff, 'is',isbrdr(iif,jjf),  &
        'got',gbrdr(iff)
      IF (isbrdr(iif,jjf) .AND. .NOT.gbrdr(iff)) THEN
        xf = iif
        yf = jjf
        EXIT
      END IF
    END IF
  END DO

  RETURN

END SUBROUTINE look_clockwise_borders

SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
! Subroutine to provide the borders of a logical array (interested in .TRUE.)

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: dx,dy,dxy
  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: Lmat
  INTEGER, DIMENSION(dxy,2), INTENT(out)                 :: brdrs
  LOGICAL, DIMENSION(dx,dy), INTENT(out)                 :: isbrdr, isbrdry

! Local
  INTEGER                                                :: i,j,ib

!!!!!!! Variables
! dx,dy: size of the space
! dxy: maximum number of border points
! Lmat: Matrix to look for the borders
! brdrs: list of coordinates of the borders
! isbrdr: matrix with .T./.F. wether the given matrix point is a border or not
! isbrdry: matrix with .T./.F. wether the given matrix point is a border or not only along y-axis

  fname = 'borders_matrixL'

  isbrdr = .FALSE.
  brdrs = -1
  ib = 1

  ! Starting with the borders. If a given point is TRUE it is a path-vertex
  ! Along y-axis
  DO i=1, dx
    IF (Lmat(i,1) .AND. .NOT.isbrdr(i,1)) THEN 
      brdrs(ib,1) = i
      brdrs(ib,2) = 1
      isbrdr(i,1) = .TRUE.
      ib=ib+1
    END IF
    IF (Lmat(i,dy) .AND. .NOT.isbrdr(i,dy)) THEN 
      brdrs(ib,1) = i
      brdrs(ib,2) = dy
      isbrdr(i,dy) = .TRUE.
      ib=ib+1
    END IF
  END DO
  ! Along x-axis
  DO j=1, dy
    IF (Lmat(1,j) .AND. .NOT.isbrdr(1,j)) THEN 
      brdrs(ib,1) = 1
      brdrs(ib,2) = j
      isbrdr(1,j) = .TRUE.
      ib=ib+1
     END IF
    IF (Lmat(dx,j) .AND. .NOT.isbrdr(dx,j)) THEN 
      brdrs(ib,1) = dx
      brdrs(ib,2) = j
      isbrdr(dx,j) = .TRUE.
      ib=ib+1
    END IF
  END DO

  isbrdry = isbrdr

  ! Border as that when looking on x-axis points with Lmat(i) /= Lmat(i+1)
  DO i=1, dx-1
    DO j=1, dy-1
      IF ( Lmat(i,j) .NEQV. Lmat(i+1,j) ) THEN 
        IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
          brdrs(ib,1) = i
          brdrs(ib,2) = j
          isbrdr(i,j) = .TRUE.
          ib=ib+1
        ELSE IF (Lmat(i+1,j) .AND. .NOT.isbrdr(i+1,j)) THEN
          brdrs(ib,1) = i+1
          brdrs(ib,2) = j
          isbrdr(i+1,j) = .TRUE.
          ib=ib+1
        END IF
      END IF
      ! y-axis
      IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN 
        IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
          brdrs(ib,1) = i
          brdrs(ib,2) = j
          isbrdr(i,j) = .TRUE.
          isbrdry(i,j) = .TRUE.
          ib=ib+1
        ELSE IF (Lmat(i,j+1) .AND. .NOT.isbrdr(i,j+1)) THEN
          brdrs(ib,1) = i
          brdrs(ib,2) = j+1
          isbrdr(i,j+1) = .TRUE.
          isbrdry(i,j+1) = .TRUE.
          ib=ib+1
        END IF
      END IF
    END DO        
  END DO

  DO i=1, dx-1
    DO j=1, dy-1
      ! y-axis
      IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN 
        IF (Lmat(i,j)) THEN
          isbrdry(i,j) = .TRUE.
        ELSE IF (Lmat(i,j+1)) THEN
          isbrdry(i,j+1) = .TRUE.
        END IF
      END IF
    END DO        
  END DO
  ! only y-axis adding bands of 2 grid points
  DO i=1, dx-1
    DO j=2, dy-2
      IF ( (Lmat(i,j) .EQV. Lmat(i,j+1)) .AND. (Lmat(i,j).NEQV.Lmat(i,j-1)) .AND. (Lmat(i,j).NEQV.Lmat(i,j+2)) ) THEN 
        IF (Lmat(i,j)) THEN
          isbrdry(i,j) = .TRUE.
          isbrdry(i,j+1) = .TRUE.
        END IF
      END IF
    END DO        
  END DO

  RETURN

END SUBROUTINE borders_matrixL

SUBROUTINE paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
! Subroutine to search the paths of a border field.

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: dx, dy, Nppt
  LOGICAL, INTENT(in)                                    :: dbg
  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: isborder
  INTEGER, DIMENSION(Nppt,2), INTENT(in)                 :: borders
  INTEGER, DIMENSION(Nppt,Nppt,2), INTENT(out)           :: paths
  INTEGER, INTENT(out)                                   :: Npath
  INTEGER, DIMENSION(Nppt), INTENT(out)                  :: Nptpaths

! Local
  INTEGER                                                :: i,j,k,ib
  INTEGER                                                :: ierr
  INTEGER                                                :: Nbrdr
  LOGICAL, DIMENSION(:), ALLOCATABLE                     :: gotbrdr, emptygotbrdr
  INTEGER                                                :: iipth, ipath, ip, Nptspath
  INTEGER                                                :: iib, jjb, iip, ijp, iif, jjf, iff
  LOGICAL                                                :: found, finishedstep

!!!!!!! Variables
! dx,dy: spatial dimensions of the space
! Nppt: possible number of paths and points that the paths can have
! isborder: boolean matrix which provide the borders of the polygon
! borders: coordinates of the borders of the polygon
! paths: coordinates of each found path
! Npath: number of paths found
! Nptpaths: number of points per path

  fname = 'paths_border'

  ! Sarting matrix
  paths = -1
  Npath = 0
  Nptspath = 0
  Nptpaths = -1

  ib=1
  finishedstep = .FALSE.

  ! Number of border points
  DO ib=1, Nppt
    IF (borders(ib,1) == -1 ) EXIT
  END DO
  Nbrdr = ib-1
    
  IF (dbg) THEN
    PRINT *,'    borders _______'
    DO i=1,Nbrdr
      PRINT *,'    ',i,':',borders(i,:)
    END DO
  END IF

  ! Matrix which keeps track if a border point has been located
  IF (ALLOCATED(gotbrdr)) DEALLOCATE(gotbrdr)
  ALLOCATE(gotbrdr(Nbrdr), STAT=ierr)
  msg = "Problems allocating matrix 'gotbrdr'"
  CALL ErrMsg(msg, fname, ierr)
  IF (ALLOCATED(emptygotbrdr)) DEALLOCATE(emptygotbrdr)
  ALLOCATE(emptygotbrdr(Nbrdr), STAT=ierr)
  msg = "Problems allocating matrix 'emptygotbrdr'"
  CALL ErrMsg(msg, fname, ierr)

  gotbrdr = .FALSE.
  emptygotbrdr = .FALSE.

  ! Starting the fun...
    
  ! Looking along the lines and when a border is found, starting from there in a clock-wise way
  iipth = 1
  ipath = 1    
  DO ib=1,Nbrdr
    iib = borders(iipth,1)
    jjb = borders(iipth,2)
    ! Starting new path
    newpath: IF (.NOT.gotbrdr(iipth)) THEN
      ip = 1
      Nptspath = 1
      paths(ipath,ip,:) = borders(iipth,:)
      gotbrdr(iipth) = .TRUE.
      ! Looking for following clock-wise search
      ! Not looking for W, because search starts from the W
      iip = iib
      ijp = jjb
      DO k=1,Nbrdr
        IF (dbg) PRINT *,ipath,'iip jip:', iip, ijp
        found = .FALSE.
        CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE.,iif,jjf,iff)
        IF (iif /= -1) THEN
          ip=ip+1
          paths(ipath,ip,:) = (/ iif,jjf /)
          found = .TRUE.
          gotbrdr(iff) = .TRUE.
          iip = iif
          ijp = jjf
          Nptspath = Nptspath + 1          
        END IF

        IF (dbg) THEN
          PRINT *,iib,jjb,'    end of this round path:', ipath, '_____', gotbrdr
          DO i=1, Nptspath
            PRINT *,'      ',i,':',paths(ipath,i,:)
          END DO
        END IF
        ! If it is not found a next point, might be because it is a non-polygon related value
        IF (.NOT.found) THEN
          IF (dbg) PRINT *,'NOT FOUND !!!', gotbrdr
          ! Are still there available borders?  
          IF (ALL(gotbrdr) .EQV. .TRUE.) THEN
            finishedstep = .TRUE.
            Npath = ipath
            Nptpaths(ipath) = Nptspath
            EXIT
          ELSE
            Nptpaths(ipath) = Nptspath
            ! Let's have a look if the previous points in the path have already some 'non-located' neighbourgs
            DO i=Nptspath,1,-1
              iip = paths(ipath,i,1)
              ijp = paths(ipath,i,2)
              CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE., iif,   &
                jjf,iff)
              IF (iif /= -1 .AND. iff /= -1) THEN
                IF (dbg) PRINT *,'    re-take path from point:', iif,',',jjf,' n-path:', iff
                found = .TRUE.
                iipth = index_list_coordsI(Nppt, borders, (/iip,ijp/))
                EXIT
              END IF
            END DO
            IF (.NOT.found) THEN
              ! Looking for the next available border point for the new path
              DO i=1,Nbrdr
                IF (.NOT.gotbrdr(i)) THEN
                  iipth = i
                  EXIT
                END IF
              END DO
              IF (dbg) PRINT *,'  Looking for next path starting at:', iipth, ' point:',              &
                borders(iipth,:)
              ipath=ipath+1
              EXIT
            END IF
          END IF
        ELSE
          IF (dbg) PRINT *,'  looking for next point...'
        END IF
        IF (finishedstep) EXIT
      END DO
    END IF newpath
  END DO
  Npath = ipath
  Nptpaths(ipath) = Nptspath
    
  DEALLOCATE (gotbrdr)
  DEALLOCATE (emptygotbrdr)

  RETURN

END SUBROUTINE paths_border

SUBROUTINE rand_sample(Nvals, Nsample, sample)
! Subroutine to randomly sample a range of indices

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: Nvals, Nsample
  INTEGER, DIMENSION(Nsample), INTENT(out)               :: sample

! Local
  INTEGER                                                :: i, ind, jmax
  REAL, DIMENSION(Nsample)                               :: randv
  CHARACTER(len=50)                                      :: fname
  LOGICAL                                                :: found
  LOGICAL, DIMENSION(Nvals)                              :: issampled
  CHARACTER(len=256)                                     :: msg
  CHARACTER(len=10)                                      :: IS1, IS2

!!!!!!! Variables
! Nvals: number of values
! Nsamples: number of samples
! sample: samnple 
  fname = 'rand_sample'

  IF (Nsample > Nvals) THEN
    WRITE(IS1,'(I10)')Nvals
    WRITE(IS2,'(I10)')Nsample
    msg = 'Sampling of ' // TRIM(IS1) // ' is too big for ' // TRIM(IS1) // 'values'
    CALL ErrMsg(msg, fname, -1)
  END IF

  ! Generation of random numbers always the same series during the whole program!
  CALL RANDOM_NUMBER(randv)

  ! Making sure that we do not repeat any value
  issampled = .FALSE.

  DO i=1, Nsample
    ! Generation of the index from the random numbers
    ind = MAX(INT(randv(i)*Nvals), 1)

    IF (.NOT.issampled(ind)) THEN
      sample(i) = ind
      issampled(ind) = .TRUE.
    ELSE
      ! Looking around the given index
      !PRINT *,' Index :', ind, ' already sampled!', issampled(ind)
      found = .FALSE.
      DO jmax=1, Nvals
        ind = MIN(ind+jmax, Nvals)
        IF (.NOT.issampled(ind)) THEN
          sample(i) = ind
          issampled(ind) = .TRUE.
          found = .TRUE.
          EXIT
        END IF
        ind = MAX(1, ind-jmax)
        IF (.NOT.issampled(ind)) THEN
          sample(i) = ind
          issampled(ind) = .TRUE.
          found = .TRUE.
          EXIT
        END IF
      END DO
      IF (.NOT.found) THEN
        msg = 'sampling could not be finished due to absence of available value!!'
        CALL ErrMsg(msg, fname, -1)
      END IF
    END IF

  END DO

  RETURN

END SUBROUTINE rand_sample

SUBROUTINE PrintQuantilesR_K(Nvals, vals, Nquants, qtvs, bspc)
! Subroutine to print the quantiles of values REAL(r_k)

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: Nvals, Nquants
  REAL(r_k), DIMENSION(Nvals), INTENT(in)                :: vals
  REAL(r_k), DIMENSION(Nquants), INTENT(in)              :: qtvs
  CHARACTER(len=1000), OPTIONAL                          :: bspc

! Local
  INTEGER                                                :: iq
  LOGICAL, DIMENSION(Nvals)                              :: search1, search2, search
  CHARACTER(len=6)                                       :: RS1
  CHARACTER(len=50)                                      :: fname
  CHARACTER(len=1000)                                    :: bspcS

!!!!!!! Variables
! vals: series of values
! qtvs: values of the quantiles
! bspc: base quantity of spaces

  fname = 'PrintQuantilesR_K'

  IF (PRESENT(bspc)) THEN
    bspcS = bspc
  ELSE
    bspcS = '      '
  END IF

  DO iq=1, Nquants-1

    WHERE (vals >= qtvs(iq)) 
      search1 = .TRUE.
    ELSEWHERE
      search1 = .FALSE.
    END WHERE

    WHERE (vals < qtvs(iq+1)) 
      search2 = .TRUE.
    ELSEWHERE
      search2 = .FALSE.
    END WHERE

    WHERE (search1 .AND. search2) 
      search = .TRUE.
    ELSEWHERE
      search = .FALSE.
    END WHERE

    WRITE(RS1, '(F6.2)')(iq)*100./(Nquants-1)
    PRINT *, TRIM(bspcS) // '[',iq,']', TRIM(RS1) // ' %:', qtvs(iq), 'N:', COUNT(search)

  END DO

  RETURN

END SUBROUTINE PrintQuantilesR_K

   INTEGER FUNCTION FindMinimumR_K(x, dsize, Startv, Endv)
! Function returns the location of the minimum in the section between Start and End.

      IMPLICIT NONE

      INTEGER, INTENT(in)                                :: dsize
      REAL(r_k), DIMENSION(dsize), INTENT(in)            :: x
      INTEGER, INTENT(in)                                :: Startv, Endv

! Local
      REAL(r_k)                                          :: Minimum
      INTEGER                                            :: Location
      INTEGER                                            :: i

      Minimum  = x(Startv)                               ! assume the first is the min
      Location = Startv                                  ! record its position
      DO i = Startv+1, Endv                              ! start with next elements
         IF (x(i) < Minimum) THEN                        !   if x(i) less than the min?
            Minimum  = x(i)                              !      Yes, a new minimum found
            Location = i                                 !      record its position
         END IF
      END DO

      FindMinimumR_K = Location                          ! return the position

   END FUNCTION  FindMinimumR_K

   SUBROUTINE SwapR_K(a, b)
! Subroutine swaps the values of its two formal arguments.

      IMPLICIT  NONE

      REAL(r_k), INTENT(INOUT)                           :: a, b
! Local
      REAL(r_k)                                          :: Temp

      Temp = a
      a    = b
      b    = Temp

   END SUBROUTINE  SwapR_K

   SUBROUTINE  SortR_K(x, Nx)
! Subroutine receives an array x() r_K and sorts it into ascending order.

      IMPLICIT NONE

      INTEGER, INTENT(IN)                                :: Nx
      REAL(r_k), DIMENSION(Nx), INTENT(INOUT)            :: x

! Local
      INTEGER                                            :: i
      INTEGER                                            :: Location

      DO i = 1, Nx-1                                     ! except for the last
         Location = FindMinimumR_K(x, Nx-i+1, i, Nx)     ! find min from this to last
         CALL  SwapR_K(x(i), x(Location))                ! swap this and the minimum
      END DO

   END SUBROUTINE  SortR_K

SUBROUTINE quantilesR_K(Nvals, vals, Nquants, quants)
! Subroutine to provide the quantiles of a given set of values of type real 'r_k'

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: Nvals, Nquants
  REAL(r_k), DIMENSION(Nvals), INTENT(in)                :: vals
  REAL(r_k), DIMENSION(Nquants), INTENT(out)             :: quants

! Local
  INTEGER                                                :: i
  REAL(r_k)                                              :: minv, maxv
  REAL(r_k), DIMENSION(Nvals)                            :: sortedvals

!!!!!!! Variables
! Nvals: number of values
! Rk: kind of real of the values
! vals: values
! Nquants: number of quants
! quants: values at which the quantile start

  minv = MINVAL(vals)
  maxv = MAXVAL(vals)

  sortedvals = vals
  ! Using from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
  CALL SortR_K(sortedvals, Nvals)

  quants(1) = minv
  DO i=2, Nquants
    quants(i) = sortedvals(INT((i-1)*Nvals/Nquants))
  END DO

END SUBROUTINE quantilesR_K


SUBROUTINE StatsR_K(Nvals, vals, minv, maxv, mean, mean2, stdev)
! Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a 
!   series of r_k numbers

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: Nvals
  REAL(r_k), DIMENSION(Nvals), INTENT(in)                :: vals
  REAL(r_k), INTENT(out)                                 :: minv, maxv, mean, mean2, stdev

!!!!!!! Variables
! Nvals: number of values
! vals: values
! minv: minimum value of values
! maxv: maximum value of values
! mean: mean value of values
! mean2: quadratic mean value of values
! stdev: standard deviation of values

  minv = MINVAL(vals)
  maxv = MAXVAL(vals)

  mean=SUM(vals)
  mean2=SUM(vals*vals)

  mean=mean/Nvals
  mean2=mean2/Nvals

  stdev=SQRT(mean2 - mean*mean)

  RETURN

END SUBROUTINE StatsR_k

END MODULE module_scientific