source: lmdz_wrf/trunk/tools/module_scientific.f90 @ 1958

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.
! fill3DI_2Dvec: Subroutine to fill a 3D integer matrix from a series of indices from a 2D matrix
! fill3DR_2Dvec: Subroutine to fill a 3D float matrix from a series of indices from a 2D matrix
! 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)
! NcountR: Subroutine to count real values
! 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
! poly_overlap_tracks_area_ascii: Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum 
!   overlaping/coincidence filtrered by a minimal area writting theoutput on an ASCII file (memory limitations)
! 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
! read_finaltrack_ascii: Subroutine to read the final trajectories from an ASCII file
! read_overlap_single_track_ascii: Subroutine to read the values for a given trajectory
! read_overlap_polys_ascii: Subroutine to read from an ASCII file the associated polygons at a given time-step
! read_overlap_tracks_ascii: Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step
! reconstruct_matrix: Subroutine to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x 
!   and y coordinates
! 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.
! write_finaltrack_ascii: Subroutine to read the final trajectories into an ASCII file
! write_overlap_polys_ascii: Subroutine to write to an ASCII file the associated polygons at a given time-step
! write_overlap_tracks_ascii: Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step

!!! *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 fill3DI_2Dvec(matind, inmat, id1, id2, od1, od2, od3, outmat)
! Subroutine to fill a 3D integer matrix from a series of indices from a 2D matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: id1, id2, od1, od2, od3
    INTEGER, DIMENSION(id1,id2), INTENT(in)              :: inmat
    INTEGER, DIMENSION(od1,od2), INTENT(in)              :: matind
    INTEGER, DIMENSION(od1,od2,od3), INTENT(out)         :: outmat

! Local
    INTEGER                                              :: i, j

!!!!!!! Variables
! matind: equivalence on the 2D space of the location in inmat
! inmat: values of the input matrix (1Dvec, time)
! id1/2: dimensions of the input matrix
! od1/3: dimensions of the output matrix
! outmat: output matrix
! NOTE: id2 == od3

    fname = 'fill3DR_2Dvec'

    DO i=1, od1 
      DO j=1, od2
        IF (matind(i,j) /= -1) THEN
          outmat(i,j,:) = inmat(matind(i,j),:)
        ELSE
          outmat(i,j,:) = fillvalI
        END IF
      END DO
    END DO

  END SUBROUTINE fill3DI_2Dvec

  SUBROUTINE fill3DR_2Dvec(matind, inmat, id1, id2, od1, od2, od3, outmat)
! Subroutine to fill a 3D float matrix from a series of indices from a 2D matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: id1, id2, od1, od2, od3
    REAL(r_k), DIMENSION(id1,id2), INTENT(in)            :: inmat
    INTEGER, DIMENSION(od1,od2), INTENT(in)              :: matind
    REAL(r_k), DIMENSION(od1,od2,od3), INTENT(out)       :: outmat

! Local
    INTEGER                                              :: i, j

!!!!!!! Variables
! matind: equivalence on the 2D space of the location in inmat
! inmat: values of the input matrix (1Dvec, time)
! id1/2: dimensions of the input matrix
! od1/3: dimensions of the output matrix
! outmat: output matrix
! NOTE: id2 == od3

    fname = 'fill3DR_2Dvec'

    DO i=1, od1 
      DO j=1, od2
        IF (matind(i,j) /= -1) THEN
          outmat(i,j,:) = inmat(matind(i,j),:)
        ELSE
          outmat(i,j,:) = fillval64
        END IF
      END DO
    END DO

  END SUBROUTINE fill3DR_2Dvec

  SUBROUTINE reconstruct_matrix(vectorXpos, vectorYpos, dvec, Xmin, Xmax, Ymin, Ymax, dmatx, dmaty,   &
    matProj, maxdiff, matind, matX, matY, matdiff)
! Subroutine to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x 
!   and y coordinates

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dvec, dmatx, dmaty
    REAL(r_k), DIMENSION(dvec), INTENT(in)               :: vectorXpos, vectorYpos
    REAL(r_k), INTENT(in)                                :: Xmin, Xmax, Ymin, Ymax, maxdiff
    CHARACTER(len=50), INTENT(in)                        :: matPRoj
    INTEGER, DIMENSION(dmatx, dmaty), INTENT(out)        :: matind
    REAL(r_k), DIMENSION(dmatx, dmaty), INTENT(out)      :: matX, matY, matdiff

! Local
    INTEGER                                              :: i,j,iv, idiff, jdiff
    REAL(r_k)                                            :: ddx, ddy, Xpos, Ypos, mindiff
    REAL(r_k), DIMENSION(dmatx,dmaty)                    :: diff
    REAL(r_k)                                            :: nXpos, xXpos, nYpos, xYpos
    INTEGER, DIMENSION(2)                                :: mindiffloc
    CHARACTER(LEN=50)                                    :: RSa, RSb

!!!!!!! Variables
! vectorXpos, vectorYpos: syncronized vectors with the x,y coordinates from which the matrix will be reconstructed
! dvec: quantitiy of coordinates to use
! Xmin,Xmax,Ymin,Ymax: Location range of the matrix to be reconstructed
! maxdiff: Authorized maximum distance between matrix location and vector location
! dmatx, dmaty: Size in grid points of the matrix to be reconstructed
! matProj: Assumed projection of the values of the matrix
!   'latlon': regular lat/lon projection
! matind: matrix with the correspondant indiced from the vector of positions
! matX, matY: matrix with the coordinates of the elements of the matrix
! matdiff: matrix with the differences at each grid point

    fname = 'reconstruct_matrix'

    matind = -oneRK
    matdiff = -oneRK

    nXpos = MINVAL(vectorXpos)
    xXpos = MAXVAL(vectorXpos)
    nYpos = MINVAL(vectorYpos)
    xYpos = MAXVAL(vectorYpos)
    PRINT *, 'Limits of the positions to localize in a matrix: longitudes:', nXpos,  &
     ' ,',xXpos, ' latitudes:', nYpos, ' ,', xYpos

    Projection: SELECT CASE(TRIM(matProj))

    CASE('latlon')
      ! Resolution along matrix axes
      ddx = (Xmax - Xmin)/(dmatx-1)
      ddy = (Ymax - Ymin)/(dmaty-1)

      ! Filling matrix positions
      DO i=1,dmatx
        Xpos = Xmin + ddx*(i-1)
        DO j=1,dmaty
          Ypos = Ymin + ddy*(j-1)
          matX(i,j) = Xpos
          matY(i,j) = Ypos
        END DO
      END DO

    CASE DEFAULT
      msg = "Projection of matrix '" // TRIM(matProj) // "' not ready " // CHAR(10) //                &
        "  Available ones: 'latlon'"
      CALL ErrMsg(msg, fname, -1)
    END SELECT Projection

    ! Let's do it properly...
    DO iv=1,dvec
      diff = SQRT((matX - vectorXpos(iv))**2 + (matY - vectorYpos(iv))**2)
      mindiffloc = MINLOC(diff)
      mindiff = MINVAL(diff)
      IF (mindiff > maxdiff) THEN
        PRINT *,'  vectorXpos:', vectorXpos(iv), 'vectorYpos:', vectorYpos(iv)
        !PRINT *,'  Xpos:', Xpos, 'Ypos:', Ypos
        WRITE(RSa, '(f20.10)')mindiff
        WRITE(RSb, '(f20.10)')maxdiff
        msg = 'Difference: ' // TRIM(RSa) // ' overpasses the maximum authorized distance: ' //   &
          TRIM(RSb)
        CALL ErrMsg(msg, fname, -1)
      ELSE
        i = mindiffloc(1)
        j = mindiffloc(2)
        matind(i,j) = iv
        matdiff(i,j) = mindiff
      END IF
    END DO

    RETURN

  END SUBROUTINE reconstruct_matrix

  SUBROUTINE read_finaltrack_ascii(funit, dt, itrack, ftrack)
! Subroutine to read the final trajectories from an ASCII file

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: funit, dt, itrack
    REAL(r_k), DIMENSION(5,dt), INTENT(out)              :: ftrack

! Local
    INTEGER                                              :: i, j, it
    LOGICAL                                              :: found

!!!!!!! Variables
! funit: unit where to write the trajectory
! dt: number of time-steps
! itrack: trajectory to read the values
! ftrack: values of the trajectory

    fname = 'read_finaltrack_ascii'

    ftrack = 0.
    
    REWIND(funit)

    it = 1
    DO WHILE (.NOT.found)

      READ(funit,10)ftrack(1,1), Str1, ((ftrack(i,j),Str1,i=2,5),Str1,j=1,dt)
      IF (INT(ftrack(1,1)) == itrack) THEN
        ftrack(1,2:dt) = ftrack(1,1)
        found = .TRUE.
      END IF

      ! Just in case
      IF (it >= dt) found = .TRUE.

    END DO

    RETURN

 10 FORMAT(I10000000,1x,A1,1x,10000000(4(F20.10,A1),A1))

  END SUBROUTINE read_finaltrack_ascii

  SUBROUTINE write_finaltrack_ascii(funit, dt, ftrack)
! Subroutine to write the final trajectories into an ASCII file

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: funit, dt
    REAL(r_k), DIMENSION(5,dt), INTENT(in)               :: ftrack

! Local
    INTEGER                                              :: i, j

!!!!!!! Variables
! funit: unit where to write the trajectory
! dt: number of time-steps
! ftrack: values of the trajectory

    fname = 'write_finaltrack_ascii'
    WRITE(funit,10)INT(ftrack(1,1)), ';', ((ftrack(i,j), ',', i=2,5), ':', j=1,dt)

    RETURN

 10 FORMAT(I10,1x,A1,1x,10000000(4(F20.10,A1),A1))

  END SUBROUTINE write_finaltrack_ascii

  SUBROUTINE read_overlap_single_track_ascii(funit, dt, Nxp, Nxtr, itrack, strack)
! Subroutine to read the values for a given trajectory

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: funit, dt, Nxp, Nxtr, itrack
    REAL(r_k), DIMENSION(5,Nxp,dt), INTENT(out)          :: strack

! Local
    INTEGER                                              :: i,j,k,l
    INTEGER                                              :: read_it, itt, it, Ntrcks
    INTEGER, DIMENSION(Nxp)                              :: Npindep
    LOGICAL                                              :: looking
    REAL(r_k), DIMENSION(5,Nxp,Nxtr)                     :: trcks

!!!!!!! Variables
! funit: unit from where retrieve the values of the trajectory
! dt: time-dimension
! Nxp: maximum allowed number of polygons per time-step
! Nxp: maximum allowed number of trajectories
! itrack: trajectory Id to look for
! strack: Values for the given trajectory

    fname = 'read_overlap_single_track_ascii'

    strack = 0.

    REWIND(funit)

    looking = .TRUE.
    itt = 0
    it = 1
    DO WHILE (looking)
      READ(funit,5,END=100)Str10, read_it

      READ(funit,*)Ntrcks
      DO i=1, Ntrcks
        READ(funit,10)l, Str1, Npindep(i), Str1, ((trcks(k,j,i),Str1,k=1,5),Str1,j=1,Npindep(i))
      END DO

      ! There is the desired trajectory at this time-step?
      IF (ANY(INT(trcks(1,1,:)) == itrack)) THEN
        itt = itt + 1
        DO i=1, Ntrcks
          IF (INT(trcks(1,1,i)) == itrack) THEN
            DO j=1, Npindep(i)
              strack(:,j,it) = trcks(:,j,i)
            END DO
          END IF
        END DO
      ELSE
        ! It trajectory has already been initialized this is the end
        IF (itt > 0) looking = .FALSE.
      END IF

      ! Just in case... ;)
      IF (read_it >= dt) looking = .FALSE.
      it = it + 1

      IF (it > dt) looking = .FALSE.

    END DO

 100 CONTINUE

    RETURN

  5 FORMAT(A10,1x,I4)
 10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1))

  END SUBROUTINE read_overlap_single_track_ascii

  SUBROUTINE read_overlap_tracks_ascii(funit, tstep, Nxp, Ntrcks, trcks)
! Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: funit, tstep, Nxp
    INTEGER, INTENT(out)                                 :: Ntrcks
    REAL(r_k), DIMENSION(5,Nxp,Nxp), INTENT(out)         :: trcks

! Local
    INTEGER                                              :: i, j, k, l, Npindep
    INTEGER                                              :: read_it

!!!!!!! Variables
! funit: unit where to write the file
! tstep: time-step to write the trajectories
! Nxp: maximum number of polygons per time-step
! Nrtcks: Number of trajectories of the given time-step
! trcks: trajectories

    fname = 'read_overlap_tracks_ascii'

    Ntrcks = 0
    trcks = 0

    READ(funit,5)Str10, read_it

    IF (read_it /= tstep) THEN
      WRITE(numSa,'(I4)')read_it
      WRITE(numSb,'(I4)')tstep
      msg = 'File time-step;' // TRIM(numSa) // ' does not coincide with the one from program:' //    &
        TRIM(numSb)
    END IF

    READ(funit,*)Ntrcks
    DO i=1, Ntrcks
      READ(funit,10)l, Str1, Npindep, Str1, ((trcks(k,j,i),Str1,k=1,5),Str1,j=1,Npindep)
    END DO

    RETURN

  5 FORMAT(A10,1x,I4)
 10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1))

  END SUBROUTINE read_overlap_tracks_ascii

  SUBROUTINE write_overlap_tracks_ascii(funit, tstep, Nxp, Ntrcks, trcks)
! Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: funit, tstep, Nxp, Ntrcks
    REAL(r_k), DIMENSION(5,Nxp,Ntrcks)                   :: trcks

! Local
    INTEGER                                              :: i, j, k, ii, Npindep, Nrealtracks

!!!!!!! Variables
! funit: unit where to write the file
! tstep: time-step to write the trajectories
! Nxp: maximum number of polygons per time-step
! Nrtcks: Number of trajectories of the given time-step
! trcks: trajectories

    fname = 'write_overlap_tracks_ascii'

    WRITE(funit,5)'time-step:', tstep

    ! Looking for the non-zero trajectories 
    Nrealtracks = 0
    DO i=1, Ntrcks
      Npindep = COUNT(trcks(2,:,i) /= zeroRK)
      IF (Npindep /= 0) Nrealtracks = Nrealtracks + 1
    END DO
    WRITE(funit,*)Nrealtracks

    ! Only writting the trajectories with values
    ii = 1
    DO i=1, Ntrcks
      Npindep = COUNT(trcks(2,:,i) /= zeroRK)
      IF (Npindep /= 0) THEN
        WRITE(funit,10)ii,';', Npindep, ';', ((trcks(k,j,i),',',k=1,5),':',j=1,Npindep)
        ii = ii + 1
      END IF
    END DO

    RETURN

  5 FORMAT(A10,1x,I4)
 10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1))

  END SUBROUTINE write_overlap_tracks_ascii

  SUBROUTINE read_overlap_polys_ascii(funit, tstep, Nxp, Nindep, SpIndep, NpIndep, pIndep)
! Subroutine to read from an ASCII file the associated polygons at a given time-step

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: funit, tstep, Nxp
    INTEGER, INTENT(out)                                 :: Nindep
    INTEGER, DIMENSION(Nxp), INTENT(out)                 :: SpIndep, NpIndep
    INTEGER, DIMENSION(Nxp,Nxp), INTENT(out)             :: pIndep

! Local
    INTEGER                                              :: i, j, k
    INTEGER                                              :: read_it

!!!!!!! Variables
! funit: unit associated to the file
! tstep: time-step of the values
! Nxp: allowed maximum numbe of polygons per time-step
! Nindpe: Number of independent polygons at this time-step
! SpIndep: Associated polygon to the independent one from the previous time-step
! NpIndep: Number of associated polygons to the independent time-step
! pIndep: polygons associated to a given independent polygon

    fname = 'read_overlap_polys_ascii'

    Nindep = 0
    SpIndep = 0
    NpIndep = 0

    READ(funit,5)Str10, read_it

    IF (read_it /= tstep) THEN
      WRITE(numSa,'(I4)')read_it
      WRITE(numSb,'(I4)')tstep
      msg = 'File time-step;' // TRIM(numSa) // ' does not coincide with the one from program:' //    &
        TRIM(numSb)
    END IF

    READ(funit,*)Nindep
    DO i=1, Nindep
      READ(funit,10) k, Str1, SpIndep(i), Str1, NpIndep(i), Str1, (pIndep(i,j), Str1, j=1,NpIndep(i))
    END DO

    RETURN

  5 FORMAT(A10,1x,I4)
 10 FORMAT(I4,1x,A1,1x,I4,1x,A1,1x,I4,A1,1x,100000(I4,A1))

  END SUBROUTINE read_overlap_polys_ascii

  SUBROUTINE write_overlap_polys_ascii(funit, tstep, Nxp, Nindep, SpIndep, NpIndep, pIndep)
! Subroutine to write into an ASCII file the associated polygons at a given time-step

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: funit, tstep, Nxp, Nindep
    INTEGER, DIMENSION(Nindep), INTENT(in)               :: SpIndep, NpIndep
    INTEGER, DIMENSION(Nindep,Nxp), INTENT(in)           :: pIndep

! Local
    INTEGER                                              :: i, j

!!!!!!! Variables
! funit: unit associated to the file
! tstep: time-step of the values
! Nxp: allowed maximum numbe of polygons per time-step
! Nindpe: Number of independent polygons at this time-step
! SpIndep: Associated polygon to the independent one from the previous time-step
! NpIndep: Number of associated polygons to the independent time-step
! pIndep: polygons associated to a given independent polygon

    fname = 'write_overlap_polys_ascii'

    WRITE(funit,5)'time-step:', tstep
    WRITE(funit,*)Nindep, ' ! Number of independent polygons'
    DO i=1, Nindep
      WRITE(funit,10) i, ';', SpIndep(i), ';', NpIndep(i), ':', (pIndep(i,j), ',', j=1,NpIndep(i))
    END DO

    RETURN

  5 FORMAT(A10,1x,I4)
 10 FORMAT(I4,1x,A1,1x,I4,1x,A1,1x,I4,A1,1x,100000(I4,A1))

  END SUBROUTINE write_overlap_polys_ascii

  SUBROUTINE poly_overlap_tracks_area_ascii(dbg, compute, dx, dy, dt, minarea, inNallpolys, allpolys, &
    ctrpolys, areapolys, Nmaxpoly, Nmaxtracks, methodmulti)
! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum 
!   overlaping/coincidence filtrered by a minimal area writting theoutput on an ASCII file (memory limitations)

    IMPLICIT NONE

    LOGICAL, INTENT(in)                                  :: dbg
    CHARACTER(LEN=*), INTENT(in)                         :: compute, methodmulti
    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

! Local
    INTEGER                                              :: i, j, ip, it, iip, itt, iit
    INTEGER                                              :: fprevunit, ftrackunit, ftrunit, ierr, ios
    LOGICAL                                              :: file_exist, dooverlap, dotracks, doftracks
    REAL(r_k), DIMENSION(Nmaxpoly)                       :: Aprevpolys, Acurrpolys
    REAL(r_k), DIMENSION(2,Nmaxpoly)                     :: Cprevpolys, Ccurrpolys
    INTEGER, DIMENSION(dx,dy)                            :: meetpolys, searchpolys
    INTEGER, DIMENSION(Nmaxpoly)                         :: coincidencies
    INTEGER, DIMENSION(Nmaxpoly)                         :: prevID, currID
    REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,2)        :: tracks
    REAL(r_k), DIMENSION(5,dt)                           :: finaltracks
    INTEGER, DIMENSION(:), ALLOCATABLE                   :: coins
    INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: coinsNpts
    INTEGER                                              :: Nmeet, Nsearch, Nindep
    INTEGER, DIMENSION(2)                                :: Nindeppolys, Npolystime
    CHARACTER(len=5)                                     :: NcoinS
    INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,2)              :: polysIndep
    INTEGER, DIMENSION(Nmaxpoly,2)                       :: NpolysIndep
    INTEGER, DIMENSION(Nmaxpoly,2)                       :: SpolysIndep
    INTEGER                                              :: iindep, iiprev
    INTEGER                                              :: Nprev, NNprev, Ntprev
    LOGICAL                                              :: Indeppolychained
    INTEGER                                              :: itrack, ictrack
    INTEGER                                              :: ixp, iyp, ttrack
    INTEGER, DIMENSION(2)                                :: Ntracks
    INTEGER                                              :: idtrack, maxtrack
    REAL(r_k), DIMENSION(5,Nmaxpoly,dt)                  :: singletrack
    REAL(r_k)                                            :: totArea, dist, mindist, maxarea, areai

!!!!!!! Variables
! dx,dy,dt: space/time dimensions
! compute: how to copmute
!   'scratch': everything from the beginning
!   'continue': skipt that parts which already have the ascii file written
! inNallpolys: Vector with the original 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
! methodmulti: methodology to follow when multiple polygons are given for the same track
!   'mean': get coordinates from the areal-weighted mean of the centers of the given polygons and their areas
!   'largest': get the coorindates of the largest polygon
!   'closest': get the coordinates of the closest polygon

    fname = 'poly_overlap_tracks_area_ascii'

    IF (dbg) PRINT *,TRIM(fname)

    SELECT CASE (TRIM(compute))
      CASE ('scratch')
        dooverlap = .TRUE.
        dotracks = .TRUE.
        doftracks = .TRUE.
      CASE ('continue')
        INQUIRE(file='polygons_overlap.dat', exist=file_exist)
        IF (.NOT.file_exist) THEN
          dooverlap = .TRUE.
        ELSE
          IF (dbg) THEN
            PRINT *, TRIM(warnmsg)
            PRINT *,"  "//TRIM(fname) // ": File 'polygons_overlap.dat' already exists, skipping it !!"
          END IF
          dooverlap = .FALSE.
        END IF
        INQUIRE(file='trajectories_overlap.dat', exist=file_exist)
        IF (.NOT.file_exist) THEN
          dotracks = .TRUE.
        ELSE
          IF (dbg) THEN
            PRINT *, TRIM(warnmsg)
            PRINT *, "  " // TRIM(fname) // ": File 'trajectories_overlap.dat' already exists, " //   &
              "skipping it !!"
          END IF
          dotracks = .FALSE.
        END IF
        INQUIRE(file='trajectories.dat', exist=file_exist)
        IF (.NOT.file_exist) THEN
          doftracks = .TRUE.
        ELSE
          IF (dbg) THEN
            PRINT *, TRIM(warnmsg)
            PRINT *,"  "//TRIM(fname) // ": File 'trajectories.dat' already exists, skipping it !!"
          END IF
          doftracks = .FALSE.
        END IF
    CASE DEFAULT
      msg = "compute case: '" // TRIM(compute) // "' not ready !!"
      CALL ErrMsg(msg, fname, -1)
    END SELECT

    ! Checking multi-polygon methodology
    IF ( (TRIM(methodmulti) /= 'mean') .AND. (TRIM(methodmulti) /= 'largest') .AND.                   &
      (TRIM(methodmulti) /= 'closest')) THEN
      msg= "methodology for multiple-polygons: '"//TRIM(methodmulti)//"' not ready" // NEW_LINE('a')//&
        " available ones: 'mean', 'largest', 'closest'"
      CALL ErrMsg(msg, fname, -1)
    END IF

    IF (dooverlap) THEN
      ! ASCII file for all the polygons and their previous associated one
      fprevunit = freeunit()
      OPEN(fprevunit, file='polygons_overlap.dat', status='new', form='formatted', iostat=ios)
      msg = "Problems opening file: 'polygons_overlap.dat'"
      IF (ios == 17) PRINT *,"  Careful: 'polygons_overlap.dat' already exists!!"
      CALL ErrMsg(msg, fname, ios)

      ! 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
          Acurrpolys(ip) = areapolys(i,it)
          Ccurrpolys(1,ip) = ctrpolys(1,i,it)
          Ccurrpolys(2,ip) = ctrpolys(2,i,it)
          NpolysIndep(ip,it) = 1
          polysIndep(ip,1,it) = i
        ELSE
          WHERE (meetpolys == i)
            meetpolys = 0
          END WHERE
        END IF
      END DO
      Nindeppolys(1) = ip
      Npolystime(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
      ! Writting to the ASCII file Independent polygons and their associated
      CALL write_overlap_polys_ascii(fprevunit,it+1, Nmaxpoly, Nindeppolys(it+1),                       &
        SpolysIndep(1:Nindeppolys(it+1),it+1), NpolysIndep(1:Nindeppolys(it+1),it+1),                   &
        polysIndep(1:Nindeppolys(it+1),:,it+1))

      it = 1
      ! Looking for the coincidencies at each time step
      DO iit=1, dt-1
        ! Number of times that a polygon has a coincidence
        coincidencies = 0
  
        ! Preparing for next time-step
        searchpolys = meetpolys
        prevID = 0
        prevID = currID
        Aprevpolys = Acurrpolys
        Cprevpolys = Ccurrpolys

        Nmeet = inNallpolys(iit+1)
        meetpolys = allpolys(:,:,iit+1)
        ip = 0
        DO i=1, Nmeet
          IF (areapolys(i,iit+1) >= minarea) THEN
            ip = ip + 1
            currID(ip) = i
            Acurrpolys(ip) = areapolys(i,iit+1)
            Acurrpolys(ip) = areapolys(i,iit+1)
            Ccurrpolys(1,ip) = ctrpolys(1,i,iit+1)
            Ccurrpolys(2,ip) = ctrpolys(2,i,iit+1)
          ELSE
            WHERE (meetpolys == i)
              meetpolys = 0
            END WHERE
          END IF
        END DO
        Nindeppolys(it+1) = ip
        Npolystime(it+1) = ip

        ! Looking throughout the independent polygons
        Nmeet = Nindeppolys(it+1)
        !Nsearch = Nindeppolys(it)
        ! Previous space might have more polygons that their number of independent ones
        Nsearch = Npolystime(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, Ccurrpolys(:,1:Nmeet),   &
          Nsearch, prevID, searchpolys, Cprevpolys(:,1:Nsearch), Aprevpolys(1:Nsearch), coins,          &
          coinsNpts)

        ! Counting the number of times that a polygon has a coincidency
        IF (dbg) THEN
          PRINT *,'  Coincidencies for the given time-step:', iit+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
                  ! Which index corresponds to this coincidence?
                  iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
                  IF (iindep == -1) THEN
                    Nindep = Nindep + 1
                    SpolysIndep(Nindep,it+1) = coins(i)
                  END IF
                  iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
                  IF (iindep < 0) THEN
                    PRINT *,'    Looking for:', coins(i)
                    PRINT *,'    Within:', SpolysIndep(1:Nindep,it+1)
                    PRINT *,'    Might content:', polysIndep(ip,1:NNprev,it)
                    PRINT *,'    From an initial list:', coins(1:Nmeet)
                    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:',iit+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

        ! Writting to the ASCII file Independent polygons and their associated
        CALL write_overlap_polys_ascii(fprevunit, iit+1, Nmaxpoly, Nindeppolys(it+1),                   &
          SpolysIndep(1:Nindeppolys(it+1),it+1), NpolysIndep(1:Nindeppolys(it+1),it+1),                 &
          polysIndep(1:Nindeppolys(it+1),:,it+1))
        ! Preparing for the next time-step
        SpolysIndep(:,it) = 0
        NpolysIndep(:,it) = 0
        polysIndep(:,:,it) = 0
        Nindeppolys(it) = Nindeppolys(it+1)
        SpolysIndep(1:Nindeppolys(it),it) = SpolysIndep(1:Nindeppolys(it+1),it+1)
        NpolysIndep(1:Nindeppolys(it),it) = NpolysIndep(1:Nindeppolys(it+1),it+1)
        Npolystime(it) = Npolystime(it+1)

        DO ip=1, Nindeppolys(it)
          polysIndep(ip,1,it) = polysIndep(ip,1,it+1)
          polysIndep(ip,2,it) = polysIndep(ip,2,it+1)
        END DO
        SpolysIndep(:,it+1) = 0
        NpolysIndep(:,it+1) = 0
        polysIndep(:,:,it+1) = 0

      END DO
      CLOSE(fprevunit)
      IF (dbg) PRINT *,"  Succesful writting of ASCII chain of polygons 'polygons_overlap.dat' !!"
    END IF
    ! ASCII file for all the polygons and their previous associated one
    fprevunit = freeunit()
    OPEN(fprevunit, file='polygons_overlap.dat', status='old', form='formatted', iostat=ios)
    msg = "Problems opening file: 'polygons_overlap.dat'"
    CALL ErrMsg(msg, fname, ios)

    it = 1
    IF (dbg) THEN
      PRINT *,  'Coincidencies to connect _______'
      DO iit=1, dt
        ! Reading from the ASCII file Independent polygons and their associated
        CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it),   &
          NpolysIndep(:,it), polysIndep(:,:,it))
        PRINT *,'  it:', iit, ' 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

    REWIND(fprevunit)

    ! 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
    IF (dotracks) THEN

      ! ASCII file for the trajectories
      ftrackunit = freeunit()
      OPEN(ftrackunit, file='trajectories_overlap.dat', status='new', form='formatted', iostat=ios)
      msg = "Problems opening file: 'trajectories_overlap.dat'" 
      IF (ios == 17) PRINT *,"  Careful: 'trajectories_overlap.dat' already exists!!"
      CALL ErrMsg(msg,fname,ios)

      ! First time-step. Take all polygons
      itrack = 0
      tracks = zeroRK
      Ntracks = 0
      it = 1
      iit = 1
      CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it),       &
        NpolysIndep(:,it), polysIndep(:,:,it))

      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

      ! Writting first time-step trajectories to the intermediate file
      CALL write_overlap_tracks_ascii(ftrackunit,iit,Nmaxpoly, Ntracks(it), tracks(:,:,1:Ntracks(it),it))

      ! Looping allover already assigned tracks
      it = 2
      maxtrack = Ntracks(1)
      timesteps: DO iit=2, dt
        CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it),     &
          NpolysIndep(:,it), polysIndep(:,:,it))
        IF (dbg) PRINT *,'track-timestep:', iit, '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
                tracks(3,iip,ictrack,it) = ctrpolys(1,iiprev,iit)
                tracks(4,iip,ictrack,it) = ctrpolys(2,iiprev,iit)
                tracks(5,iip,ictrack,it) = iit
              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 = maxtrack + 1
            IF (dbg) PRINT *,'  New track!', maxtrack

            ! 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
              tracks(2,j,ictrack,it) = iiprev
              tracks(3,j,ictrack,it) = ctrpolys(1,iiprev,iit)
              tracks(4,j,ictrack,it) = ctrpolys(2,iiprev,iit)
              tracks(5,j,ictrack,it) = iit
            END DO
          END IF

        END DO current_poly

        IF (dbg) THEN
          PRINT *,'  At this time-step:', iit, ' 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

        CALL write_overlap_tracks_ascii(ftrackunit,iit,Nmaxpoly,Ntracks(it),tracks(:,:,1:Ntracks(it),it))
        ! Re-initializing for the next time-step
        tracks(:,:,:,it-1) = zeroRK
        Ntracks(it-1) = Ntracks(it)
        tracks(:,:,1:Ntracks(it-1),it-1) = tracks(:,:,1:Ntracks(it),it)
        Ntracks(it) = 0
        tracks(:,:,:,it) = zeroRK

      END DO timesteps
      CLOSE(ftrackunit)
      IF (dbg) PRINT *,"  Succesful writting of ASCII chain of polygons 'trajectories_overlap.dat' !!"
      CLOSE(fprevunit)
    END IF

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

    IF (doftracks) THEN
      ! ASCII file for the trajectories
      ftrackunit = freeunit()
      OPEN(ftrackunit, file='trajectories_overlap.dat', status='old', form='formatted', iostat=ios)
      msg = "Problems opening file: 'trajectories_overlap.dat'" 
      CALL ErrMsg(msg,fname,ios)

      ! ASCII file for the final trajectories
      ftrunit = freeunit()
      OPEN(ftrunit, file='trajectories.dat', status='new', form='formatted', iostat=ios)
      msg = "Problems opening file: 'trajectories.dat'" 
      IF (ios == 17) PRINT *,"  Careful: 'trajectories.dat' already exists!!"
      CALL ErrMsg(msg,fname,ios)

      finaltracks = zeroRK

      DO itt=1, Nmaxtracks
        CALL read_overlap_single_track_ascii(ftrackunit, dt, Nmaxpoly, Nmaxtracks, itt, singletrack)

        ! It might reach the las trajectory
        IF (ALL(singletrack == zeroRK)) EXIT

        itrack = INT(MAXVAL(singletrack(1,1,:)))
        IF (dbg) THEN
          PRINT *,'  Trajectory:', itt, '_______', itrack
          DO it=1, dt
            IF (singletrack(2,1,it) /= zeroRK) THEN
              j = COUNT(singletrack(2,:,it) /= zeroRK)
              PRINT *,it,':',(singletrack(3,i,it),',',singletrack(4,i,it),' ; ',i=1,j)
            END IF
          END DO
        END IF

        finaltracks = zeroRK
        finaltracks(1,:) = itrack*oneRK
        DO it =1, dt      
          Nprev = COUNT(INT(singletrack(2,:,it)) /= zeroRK)
          IF (Nprev /= 0) THEN
            finaltracks(5,it) = it*oneRK
            IF (TRIM(methodmulti) == 'largest') THEN
              maxarea = -10.*oneRK
              DO ip=1, Nprev
                IF (areapolys(singletrack(2,ip,it),it) > maxarea) THEN
                  maxarea = areapolys(singletrack(2,ip,it),it)
                  i = ip
                END IF
              END DO
              IF (dbg) THEN
                PRINT *,'  Determine the trajectory coordinates to the largest polygon:', i,          &
                  ' area:', maxarea
              END IF
              finaltracks(2,it) = singletrack(2,i,it)*oneRK
              finaltracks(3,it) = singletrack(3,i,it)
              finaltracks(4,it) = singletrack(4,i,it)
            ELSE IF (TRIM(methodmulti) == 'closest') THEN
              IF (it > 1) THEN
                mindist = 10000000.*oneRK
                DO ip=1, Nprev
                  dist = SQRT((singletrack(3,ip,it)-finaltracks(3,it-1))**2 +                         &
                    (singletrack(4,ip,it)-finaltracks(4,it-1))**2 )
                  IF (dist < mindist) THEN
                    mindist = dist
                    i = ip
                  END IF
                END DO
                finaltracks(2,it) = singletrack(3,i,it)*oneRK
                finaltracks(3,it) = singletrack(3,i,it)
                finaltracks(4,it) = singletrack(4,i,it)
                IF (dbg) THEN
                  PRINT *,'  Determine the trajectory coordinates to the closest previous polygon:',i,&
                    ' distance:', mindist
                END IF
              ELSE
                maxarea = -10.*oneRK
                DO ip=1, Nprev
                  IF (areapolys(singletrack(2,ip,it),it) > maxarea) THEN
                    maxarea = areapolys(singletrack(2,ip,it),it)
                    i = ip
                  END IF
                END DO
                IF (dbg) THEN
                  PRINT *, '  Determine the trajectory coordinates to the largest polygon:', i,        &
                    ' area:', maxarea, ' at the first time-step then to the closest'
                END IF
                finaltracks(2,it) = i*oneRK
                finaltracks(3,it) = singletrack(3,i,it)
                finaltracks(4,it) = singletrack(4,i,it)              
              END IF
            ELSE
              totArea = zeroRK
              finaltracks(2,it) = -oneRK
              finaltracks(3,it) = zeroRK
              finaltracks(4,it) = zeroRK
              DO ip=1, Nprev
                areai = areapolys(singletrack(2,ip,it),it)
                totArea = totArea + areai
                finaltracks(3,it) = finaltracks(3,it) + singletrack(3,ip,it)*areai
                finaltracks(4,it) = finaltracks(4,it) + singletrack(4,ip,it)*areai
              END DO
              finaltracks(3,it) = finaltracks(3,it)/totArea
              finaltracks(4,it) = finaltracks(4,it)/totArea
              IF (dbg) THEN
                PRINT *,'  Determine the trajectory coordinates to the area-averaged polygon ' //     &
                  ' total area:', totArea
              END IF

            END IF

          END IF
        END DO
        ! Writting the final track into the ASCII file
        CALL write_finaltrack_ascii(ftrunit, dt, finaltracks)

      END DO
      CLOSE(ftrackunit)
      IF (dbg) PRINT *,"  Succesful writting of ASCII trajectories 'trajectories.dat' !!"
      CLOSE(ftrunit)
    END IF

    IF (ALLOCATED(coins)) DEALLOCATE(coins)
    IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts)

    RETURN

  END SUBROUTINE poly_overlap_tracks_area_ascii

  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(Nmaxpoly)                       :: Aprevpolys, Acurrpolys
    REAL(r_k), DIMENSION(2,Nmaxpoly)                     :: Cprevpolys, Ccurrpolys
    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                                              :: 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
    REAL(r_k)                                            :: ixp, iyp
    INTEGER                                              :: 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)

    ! 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
        Acurrpolys(ip) = areapolys(i,it)
        Ccurrpolys(1,ip) = ctrpolys(1,i,it)
        Ccurrpolys(2,ip) = ctrpolys(2,i,it)
        NpolysIndep(ip,it) = 1
        polysIndep(ip,1,it) = i
      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
      Aprevpolys = Acurrpolys
      Cprevpolys = Ccurrpolys
      ip = 0

      DO i=1, Nmeet
        IF (areapolys(i,it+1) >= minarea) THEN
          ip = ip + 1
          currID(ip) = i
          Acurrpolys(ip) = areapolys(i,it+1)
          Acurrpolys(ip) = areapolys(i,it+1)
          Ccurrpolys(1,ip) = ctrpolys(1,i,it+1)
          Ccurrpolys(2,ip) = ctrpolys(2,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, Acurrpolys(1:Nmeet),      &
        Nsearch, prevID, searchpolys, Cprevpolys(:,1:Nsearch), Aprevpolys(1:Nsearch), coins,          &
        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
                ! Which index corresponds to this coincidence?
                iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
                IF (iindep == -1) THEN
                  Nindep = Nindep + 1
                  SpolysIndep(Nindep,it+1) = coins(i)
                END IF
                iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
                IF (iindep < 0) THEN
                  PRINT *,'    Looking for:', coins(i)
                  PRINT *,'    Within:', SpolysIndep(1:Nindep,it+1)
                  PRINT *,'    Might content:', polysIndep(ip,1:NNprev,it)
                  PRINT *,'    From an initial list:', coins(1:Nmeet)
                  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)
        finaltracks(1,itrack,it) = itrack*1.
        finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev*1.
        finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev*1.
        finaltracks(4,itrack,it) = it*1.
      END DO
    END DO

    DEALLOCATE(coins)
    DEALLOCATE(coinsNpts)

    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 largest 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
! IDallpoly: ID of the polygon to find coincidence
! 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)
! IDpolys: ID of the polygon to search for coincidences
! 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
! Lluis
    INTEGER                                              :: Ndiffvals
    INTEGER, DIMENSION(:), ALLOCATABLE                   :: diffvals

!!!!!!! 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(I7,1x))',polys(i,:)
      END DO
    END IF

    IF (ALLOCATED(diffvals)) DEALLOCATE(diffvals)
    ALLOCATE(diffvals(dx*dy))

    ! Checking for consistency on number of polygons and real content (except 0 value)
    CALL Nvalues_2DArrayI(dx, dy, dx*dy, polys, Ndiffvals, diffvals)
    IF (Ndiffvals -1 /= Npoly) THEN
      PRINT *,TRIM(emsg)
      PRINT *,'    number of different values:', Ndiffvals-1, ' theoretical Npoly:', Npoly
      PRINT *,'    Different values:', diffvals(1:Ndiffvals)
      msg = 'Number of different values and Npoly must coincide'
      CALL ErrMsg(msg, fname, -1)
    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 largest 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(I7,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
!   'ctsarea': Constant Area
!   '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.
  xcorr = 0.
  ycorr = 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
!            IF (KIND(xcorr(i,j)) == KIND(1.d0)) THEN
!              xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
!            ELSE
              xcorr(i,j) = ABS(SIN(lon(i,j)*DegRad))
!            END IF
            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
!            IF (KIND(xcorr(i,j)) == KIND(1.d0)) THEN
!              xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
!            ELSE
              xcorr(i,j) = ABS(SIN(lon(i,j)*DegRad))
!            END IF
            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 *,'  grid_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
!  IF (KIND(polyvals(1)) == KIND(1.d0)) THEN
!    distvn = DABS(polyvals - nval)
!    distvx = DABS(polyvals - xval)
!  ELSE
    distvn = ABS(polyvals - nval)
    distvx = ABS(polyvals - xval)
!  END IF

  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

  SUBROUTINE NcountR(values, d1, Ndiffvals, counts)
! Subroutine to count real values

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    REAL(r_k), DIMENSION(d1), INTENT(in)                 :: values
    INTEGER, INTENT(out)                                 :: Ndiffvals
    REAL(r_k), DIMENSION(d1,2), INTENT(out)              :: counts
! Local
    INTEGER                                              :: i, ival
    REAL(r_k), DIMENSION(d1)                             :: diffv

!!!!!!! Variables
! values: values to count
! counts: counts of time for each value
    
    fname = 'NcountR'

    counts = -1.

    counts(1,1) = values(1)
    counts(1,2) = 1
    Ndiffvals = 1
    DO i=2,d1
      diffv(1:Ndiffvals) = counts(1:Ndiffvals,1) - values(i)
      IF (ANY(diffv(1:Ndiffvals) == 0)) THEN
        ival = Index1DArrayR(counts(1:Ndiffvals,1), Ndiffvals, values(i))
        counts(ival,2) = counts(ival,2) + 1
      ELSE
        Ndiffvals = Ndiffvals + 1
        counts(Ndiffvals,1) = values(i)
        counts(Ndiffvals,2) = 1
      END IF
    END DO

  END SUBROUTINE NcountR

END MODULE module_scientific