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

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
! coincident_gridsin2D: Subroutine to determine which lists of 2D gridsin points of an A list are also found in a B list
! coincident_list_2Dcoords: Subroutine to determine which 2D points of an A list are also found in a B list
! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
! coincident_polygon: Subroutine to provide the intersection polygon between two polygons
! distanceRK: Function to provide the distance between two points
! 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
! grid_within_polygon: Subroutine to determine which grid cells from a matrix lay inside a polygon
! grid_spacepercen: Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another 
!   series of grid-cells (A)
! grid_spacepercen_providing_polys: Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another 
!   series of grid-cells (A) providing coincident polygons
! grid_spacepercen_within_reg: Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed 
!   by a polygon
! grid_spacepercen_within_reg_providing_polys: Subroutine to compute the percentage of grid space of a series of grid 
!   cells which are encompassed by a polygon providing coordinates of the resultant polygons
! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon 
!   following ray casting algorithm
! Heron_area_triangle: Subroutine to compute area of a triangle using Heron's formula
! intersectfaces: Subroutine to provide if two faces of two polygons intersect
! intersection_2Dlines: Subroutine to provide the intersection point between two lines on the plane using Cramer's method
! join_polygon: Subroutine to join two polygons
! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders 
!   (limits of a region)
! multi_index_mat2DI: Subroutine to provide the indices of the different locations of a value inside a 2D integer matrix
! multi_index_mat3DI: Subroutine to provide the indices of the different locations of a value inside a 3D integer matrix
! multi_index_mat4DI: Subroutine to provide the indices of the different locations of a value inside a 4D integer matrix
! multi_index_mat2DRK: Subroutine to provide the indices of the different locations of a value inside a 2D RK matrix
! multi_index_mat3DRK: Subroutine to provide the indices of the different locations of a value inside a 3D RK matrix
! multi_index_mat4DRK: Subroutine to provide the indices of the different locations of a value inside a 4D RK matrix
! multi_spaceweightstats_in1DRKno_slc3v3: Subroutine to compute an spatial statistics value from a 1D RK matrix without
!   running one into a matrix of 3-variables slices of rank 3 using spatial weights
! multi_spaceweightstats_in2DRKno_slc3v3: Subroutine to compute an spatial statistics value from a 2D RK matrix without
!   running one into a matrix of 3-variables slices of rank 3 using spatial weights
! multi_spaceweightstats_in3DRK3_slc3v3: Subroutine to compute an spatial statistics value from a 3D RK matrix using 
!   3rd dimension as running one into a matrix of 3-variables slices of rank 3 using spatial weights
! multi_spaceweightstats_in3DRK3_slc3v4: Subroutine to compute an spatial statistics value from a 3D RK matrix using 
!   3rd dimension as running one into a matrix of 3-variables slices of rank 4 using spatial weights
! 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
! percentiles_R_K*D: Subroutine to compute the percentiles of a *D R_K array along given set of axis
! point_in_face: Function to determine if a given point is on a face of a polygon
! point_inside: Function to determine if a given point is inside a polygon providing its sorted vertices
! poly_has_point: Function to determine if a polygon has already a given point as one of its vertex
! 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
! runmean_F1D: Subroutine fo computing the running mean of a given set of float 1D values
! shoelace_area_polygon: Computing the area of a polygon using sholace formula
! sort_polygon: Subroutine to sort a polygon using its center as average of the coordinates and remove duplicates
! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
! spacepercen: Subroutine to compute the space-percentages of a series of grid cells (B) into another series of grid-cells (A)
! spaceweightstats: Subroutine to compute an spatial statistics value from a matrix B into a matrix A using weights
! 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.
! unique_matrixRK2D: Subroutine to provide the unique values within a 2D RK matrix
! 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

  SUBROUTINE runmean_F1D(d1, values, Nmean, headertail, runmean)
! Subroutine fo computing the running mean of a given set of float 1D values

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: d1, Nmean
  REAL(r_k), DIMENSION(d1), INTENT(in)                   :: values
  CHARACTER(len=*), INTENT(in)                           :: headertail
  REAL(r_k), DIMENSION(d1), INTENT(out)                  :: runmean
  
! Local
  INTEGER                                                :: i, j, Nmean2
  CHARACTER(len=5)                                       :: NmeanS

!!!!!!! Variables
! values: values to compute the running mean
! Nmean: number of odd points to use for the running mean
! headertail: How to proceed for the grid points at the beginning of the values which are not 
!   encompassed by the Nmean
!   'miss': set as missing values (1.d20)
!   'original': keep the original values
!   'progressfill': mean the values as a progressive running filter (e.g. for Nmean=5):
!     runmean[values(1)] = values(1)
!     runmean[values(2)] = MEAN(values(1:3))
!     runmean[values(3)] = MEAN(values(1:5))
!     runmean[values(4)] = MEAN(values(2:6))
!     (...)
!     runmean[values(d1-2)] = MEAN(values(d1-5:d1))
!     runmean[values(d1-1)] = MEAN(values(d1-2:d1))
!     runmean[values(d1)] = MEAN(values(dd1))
!   'zero': set as zero values
! runmean: runnig mean values

  fname = 'runmean_F1D'

  IF (MOD(Nmean,2) == 0) THEN
    WRITE(NmeanS,'(I5)')Nmean
    msg="Nmean has to be odd!! value provided: "// NmeanS
    CALL ErrMsg(msg, fname, -1)
  END IF
  Nmean2 = Nmean/2
 
  SELECT CASE (TRIM(headertail))
    CASE ('missing')
      runmean = fillval64
    CASE ('original')
      runmean = values
    CASE ('progressfill')
      DO i=1, Nmean2
        runmean(i) = SUM(values(1:2*(i-1)+1))/(2*(i-1)+1)
      END DO
      runmean(d1) = values(d1)
      DO i=2, Nmean2
        j = d1-(2*(i-1))
        runmean(d1-(i-1)) = SUM(values(j:d1))/(2*(i-1)+1)
      END DO
    CASE ('zero')
      runmean = zeroRK
    CASE DEFAULT
      msg = "'" // TRIM(headertail) // "' not available !!" //CHAR(44) // "  available ones: " //     &
        "'missing', 'original', 'progressfill', 'zero'"
      CALL ErrMsg(msg, fname, -1)
  END SELECT

  DO i= 1+Nmean2, d1 - Nmean2
    runmean(i) = SUM(values(i-Nmean2:i+Nmean2))/Nmean
  END DO

  END SUBROUTINE runmean_F1D

  SUBROUTINE percentiles_R_K2D(values, axisS, Npercen, d1, d2, percentiles)
  ! Subroutine to compute the percentiles of a 2D R_K array along given set of axis

    IMPLICIT NONE
  
    INTEGER, INTENT(in)                                    :: d1, d2, Npercen
    REAL(r_k), DIMENSION(d1,d2), INTENT(in)                :: values
    CHARACTER(LEN=*), INTENT(in)                           :: axisS
    REAL(r_k), DIMENSION(d1, d2, Npercen), INTENT(out)     :: percentiles

    ! Local
    INTEGER                                                :: i
    INTEGER                                                :: Lstring, LaxisS, iichar
    CHARACTER(LEN=1000)                                    :: splitaxis
    INTEGER, DIMENSION(1)                                  :: axis1
    CHARACTER(LEN=200), DIMENSION(2)                       :: axis2S
    INTEGER, DIMENSION(2)                                  :: axis2
    CHARACTER(LEN=1)                                       :: Naxs

!!!!!!! Variables
! d1,d2: length of the 2D dimensions
! values: values to use to compute the percentiles
! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes)
! Npercen: number of percentiles
! percentiles: percentiles of the daata

    fname = 'percentiles_R_K2D'

    LaxisS = LEN_TRIM(axisS)
    iichar = numberTimes(axisS(1:LaxisS), ':')

    splitaxis = ''
    splitaxis(1:LaxisS) = axisS(1:LaxisS)
    percentiles = 0.

    IF (iichar == 0) THEN
      READ(axisS,'(I1)')axis1(1)
    ELSE IF (iichar == 1) THEN
      CALL split(splitaxis, ':', 2, axis2S)
    ELSE
      WRITE(Naxs,'(A1)')iichar
      msg = "' rank 2 values can not compute percentiles using " // Naxs // "' number of axis !!"
      CALL ErrMsg(msg, fname, -1)
    END IF

    IF (TRIM(axisS) == 'all') iichar = 2
  
    IF (iichar == 0) THEN
      ! Might be a better way, but today I can't think it !!
      IF (axis1(1) == 1) THEN
        DO i=1, d2
          CALL quantilesR_K(d1, values(:,i), Npercen, percentiles(1,i,:))
        END DO
      ELSE IF (axis1(1) == 2) THEN
        DO i=1, d1
          CALL quantilesR_K(d2, values(i,:), Npercen, percentiles(i,1,:))
        END DO
      ELSE
        WRITE(Naxs,'(A1)')axis1(1)
        msg = "' rank 2 values can not compute percentiles using axis " // Naxs // "' !!"
        CALL ErrMsg(msg, fname, -1)
      END IF
    ELSE
      CALL quantilesR_K(d1*d2, RESHAPE(values, (/d1*d2/)), Npercen, percentiles(1,1,:))
    END IF

  END SUBROUTINE percentiles_R_K2D

  SUBROUTINE percentiles_R_K3D(values, axisS, Npercen, d1, d2, d3, percentiles)
  ! Subroutine to compute the percentiles of a 3D R_K array along given set of axis

    IMPLICIT NONE
  
    INTEGER, INTENT(in)                                    :: d1, d2, d3, Npercen
    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)             :: values
    CHARACTER(LEN=*), INTENT(in)                           :: axisS
    REAL(r_k), DIMENSION(d1, d2, d3, Npercen), INTENT(out) :: percentiles

    ! Local
    INTEGER                                                :: i, j
    INTEGER                                                :: Lstring, LaxisS, iichar
    CHARACTER(LEN=1000)                                    :: splitaxis
    INTEGER, DIMENSION(1)                                  :: axis1
    CHARACTER(LEN=200), DIMENSION(2)                       :: axis2S
    INTEGER, DIMENSION(2)                                  :: axis2
    CHARACTER(LEN=200), DIMENSION(3)                       :: axis3S
    INTEGER, DIMENSION(3)                                  :: axis3
    CHARACTER(LEN=1)                                       :: Naxs1, Naxs2

!!!!!!! Variables
! d1,d2: length of the 2D dimensions
! values: values to use to compute the percentiles
! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes)
! Npercen: number of percentiles
! percentiles: percentiles of the daata

    fname = 'percentiles_R_K3D'

    LaxisS = LEN_TRIM(axisS)
    iichar = numberTimes(axisS(1:LaxisS), ':')

    splitaxis = ''
    splitaxis(1:LaxisS) = axisS(1:LaxisS)

    percentiles = 0.

    IF (iichar == 0) THEN
      READ(axisS,'(I1)')axis1(1)
    ELSE IF (iichar == 1) THEN
      CALL split(splitaxis, ':', 2, axis2S)
      DO i=1,2
        READ(axis2S(i), '(I1)')axis2(i)
      END DO
    ELSE IF (iichar == 2) THEN
      CALL split(splitaxis, ':', 3, axis3S)
    ELSE
      READ(Naxs1,'(A1)')iichar
      msg = "' rank 3 values can not compute percentiles using " // Naxs1 // "' number of axis !!"
      CALL ErrMsg(msg, fname, -1)
    END IF

    IF (TRIM(axisS) == 'all') iichar = 3
  
    IF (iichar == 0) THEN
      ! Might be a better way, but today I can't think it !!
      IF (axis1(1) == 1) THEN
        DO i=1, d2
          DO j=1, d3
            CALL quantilesR_K(d1, values(:,i,j), Npercen, percentiles(1,i,j,:))
          END DO
        END DO
      ELSE IF (axis1(1) == 2) THEN
        DO i=1, d1
          DO j=1, d3
            CALL quantilesR_K(d2, values(i,:,j), Npercen, percentiles(i,1,j,:))
          END DO
        END DO
      ELSE IF (axis1(1) == 3) THEN
        DO i=1, d1
          DO j=1, d2
            CALL quantilesR_K(d3, values(i,j,:), Npercen, percentiles(i,j,1,:))
          END DO
        END DO
      ELSE
        WRITE(Naxs1,'(A1)')axis1(1)
        msg = "' rank 3 values can not compute percentiles using axis " // Naxs1 // "' !!"
        CALL ErrMsg(msg, fname, -1)
      END IF
    ELSE IF (iichar == 1) THEN
      ! Might be a better way, but today I can't think it !!
      IF (axis2(1) == 1 .AND. axis2(2) == 2) THEN
        DO i=1, d3
          CALL quantilesR_K(d1*d2, RESHAPE(values(:,:,i), (/d1*d2/)), Npercen, percentiles(1,1,i,:))
        END DO
      ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3) THEN
        DO i=1, d2
          CALL quantilesR_K(d1*d3, RESHAPE(values(:,i,:), (/d1*d3/)), Npercen, percentiles(1,i,1,:))
        END DO
      ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3) THEN
        DO i=1, d1
          CALL quantilesR_K(d2*d3, RESHAPE(values(i,:,:), (/d2*d3/)), Npercen, percentiles(i,1,1,:))
        END DO
      ELSE
        WRITE(Naxs1,'(A1)')axis2(1)
        WRITE(Naxs2,'(A1)')axis2(2)
        msg="' rank 3 values can not compute percentiles using axis "//Naxs1// ', ' // Naxs2 // "' !!"
        CALL ErrMsg(msg, fname, -1)
      END IF
    ELSE
      CALL quantilesR_K(d1*d2*d3, RESHAPE(values, (/d1*d2*d3/)), Npercen, percentiles(1,1,1,:))
    END IF

  END SUBROUTINE percentiles_R_K3D

  SUBROUTINE percentiles_R_K4D(values, axisS, Npercen, d1, d2, d3, d4, percentiles)
  ! Subroutine to compute the percentiles of a 4D R_K array along given set of axis

    IMPLICIT NONE
  
    INTEGER, INTENT(in)                                    :: d1, d2, d3, d4, Npercen
    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)          :: values
    CHARACTER(LEN=*), INTENT(in)                           :: axisS
    REAL(r_k), DIMENSION(d1,d2,d3,d4,Npercen), INTENT(out) :: percentiles

    ! Local
    INTEGER                                                :: i, j, k
    INTEGER                                                :: Lstring, LaxisS, iichar
    CHARACTER(LEN=1000)                                    :: splitaxis
    INTEGER, DIMENSION(1)                                  :: axis1
    CHARACTER(LEN=200), DIMENSION(2)                       :: axis2S
    INTEGER, DIMENSION(2)                                  :: axis2
    CHARACTER(LEN=200), DIMENSION(3)                       :: axis3S
    INTEGER, DIMENSION(3)                                  :: axis3
    CHARACTER(LEN=200), DIMENSION(4)                       :: axis4S
    CHARACTER(LEN=1)                                       :: Naxs1, Naxs2, Naxs3

!!!!!!! Variables
! d1,d2: length of the 2D dimensions
! values: values to use to compute the percentiles
! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes)
! Npercen: number of percentiles
! percentiles: percentiles of the daata

    fname = 'percentiles_R_K3D'

    LaxisS = LEN_TRIM(axisS)
    iichar = numberTimes(axisS(1:LaxisS), ':')

    splitaxis = ''
    splitaxis(1:LaxisS) = axisS(1:LaxisS)

    percentiles = 0.

    PRINT *,'iichar:', iichar, axisS(1:LaxisS)

    IF (iichar == 0) THEN
      READ(axisS,'(I1)')axis1(1)
    ELSE IF (iichar == 1) THEN
      CALL split(splitaxis, ':', 2, axis2S)
      DO i=1,2
        READ(axis2S(i), '(I1)')axis2(i)
      END DO
    ELSE IF (iichar == 2) THEN
      CALL split(splitaxis, ':', 3, axis3S)
      DO i=1,3
        READ(axis3S(i), '(I1)')axis3(i)
      END DO
    ELSE IF (iichar == 3) THEN
      CALL split(splitaxis, ':', 4, axis4S)
    ELSE
      READ(Naxs1,'(A1)')iichar
      msg = "' rank 4 values can not compute percentiles using " // Naxs1 // "' number of axis !!"
      CALL ErrMsg(msg, fname, -1)
    END IF

    IF (TRIM(axisS) == 'all') iichar = 4
  
    IF (iichar == 0) THEN
      ! Might be a better way, but today I can't think it !!
      IF (axis1(1) == 1) THEN
        DO i=1, d2
          DO j=1, d3
            DO k=1, d4
              CALL quantilesR_K(d1, values(:,i,j,k), Npercen, percentiles(1,i,j,k,:))
            END DO
          END DO
        END DO
      ELSE IF (axis1(1) == 2) THEN
        DO i=1, d1
          DO j=1, d3
            DO k=1, d4
              CALL quantilesR_K(d2, values(i,:,j,k), Npercen, percentiles(i,1,j,k,:))
            END DO
          END DO
        END DO
      ELSE IF (axis1(1) == 3) THEN
        DO i=1, d1
          DO j=1, d2
            DO k=1, d4
              CALL quantilesR_K(d3, values(i,j,:,k), Npercen, percentiles(i,j,1,k,:))
            END DO
          END DO
        END DO
      ELSE IF (axis1(1) == 4) THEN
        DO i=1, d1
          DO j=1, d2
            DO k=1, d3
              CALL quantilesR_K(d4, values(i,j,k,:), Npercen, percentiles(i,j,k,1,:))
            END DO
          END DO
        END DO
      ELSE
        WRITE(Naxs1,'(A1)')axis1(1)
        msg = "' rank 3 values can not compute percentiles using axis " // Naxs1 // "' !!"
        CALL ErrMsg(msg, fname, -1)
      END IF
    ELSE IF (iichar == 1) THEN
      ! Might be a better way, but today I can't think it !!
      IF (axis2(1) == 1 .AND. axis2(2) == 2) THEN
        DO i=1, d3
          DO j=1, d4
            CALL quantilesR_K(d1*d2, RESHAPE(values(:,:,i,j), (/d1*d2/)), Npercen,                    &
              percentiles(1,1,i,j,:))
          END DO
        END DO
      ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3) THEN
        DO i=1, d2
          DO j=1, d4
            CALL quantilesR_K(d1*d3, RESHAPE(values(:,i,:,j), (/d1*d3/)), Npercen,                    &
              percentiles(1,i,1,j,:))
          END DO
        END DO
      ELSE IF (axis2(1) == 1 .AND. axis2(2) == 4) THEN
        DO i=1, d2
          DO j=1, d3
            CALL quantilesR_K(d1*d4, RESHAPE(values(:,i,j,:), (/d1*d4/)), Npercen,                    &
              percentiles(1,i,j,1,:))
          END DO
        END DO
      ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3) THEN
        DO i=1, d1
          DO j=1, d4
            CALL quantilesR_K(d2*d3, RESHAPE(values(i,:,:,j), (/d2*d3/)), Npercen,                    &
              percentiles(i,1,1,j,:))
          END DO
        END DO
      ELSE IF (axis2(1) == 2 .AND. axis2(2) == 4) THEN
        DO i=1, d1
          DO j=1, d3
            CALL quantilesR_K(d2*d4, RESHAPE(values(i,:,j,:), (/d2*d4/)), Npercen,                    &
              percentiles(i,1,j,1,:))
          END DO
        END DO
      ELSE IF (axis2(1) == 3 .AND. axis2(2) == 4) THEN
        DO i=1, d1
          DO j=1, d2
            CALL quantilesR_K(d3*d4, RESHAPE(values(i,j,:,:), (/d3*d4/)), Npercen,                    &
              percentiles(i,j,1,1,:))
          END DO
        END DO
      ELSE
        WRITE(Naxs1,'(A1)')axis2(1)
        WRITE(Naxs2,'(A1)')axis2(2)
        msg="' rank 4 values can not compute percentiles using axis "//Naxs1// ', ' // Naxs2 // "' !!"
        CALL ErrMsg(msg, fname, -1)
      END IF
    ELSE IF (iichar == 2) THEN
      IF (axis2(1) == 1 .AND. axis2(2) == 2 .AND. axis3(3) == 3) THEN
        DO i=1, d4
          CALL quantilesR_K(d1*d2*d3, RESHAPE(values(:,:,:,i), (/d1*d2*d3/)), Npercen,                &
            percentiles(1,1,1,i,:))
        END DO
      ELSE IF (axis2(1) == 1 .AND. axis2(2) == 2 .AND. axis3(3) == 4) THEN
        DO i=1, d3
          CALL quantilesR_K(d1*d2*d4, RESHAPE(values(:,:,i,:), (/d1*d2*d4/)), Npercen,                &
            percentiles(1,1,i,1,:))
        END DO
      ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3 .AND. axis3(3) == 4) THEN
        DO i=1, d2
          CALL quantilesR_K(d1*d3*d4, RESHAPE(values(:,i,:,:), (/d1*d3*d4/)), Npercen,                &
            percentiles(1,i,1,1,:))
        END DO
      ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3 .AND. axis3(3) == 4) THEN
        DO i=1, d1
          CALL quantilesR_K(d2*d3*d4, RESHAPE(values(i,:,:,:), (/d2*d3*d4/)), Npercen,                &
            percentiles(i,1,1,1,:))
        END DO
      ELSE
        WRITE(Naxs1,'(A1)')axis3(1)
        WRITE(Naxs2,'(A1)')axis3(2)
        WRITE(Naxs3,'(A1)')axis3(2)
        msg="' rank 4 values can not compute percentiles using axis "// Naxs1 // ', ' // Naxs2 //     &
          ', ' // Naxs3 //"' !!"
        CALL ErrMsg(msg, fname, -1)
      END IF
    ELSE
      CALL quantilesR_K(d1*d2*d3*d4, RESHAPE(values, (/d1*d2*d3*d4/)), Npercen, percentiles(1,1,1,1,:))
    END IF

  END SUBROUTINE percentiles_R_K4D

  REAL(r_k) FUNCTION distanceRK(pointA, pointB)
  ! Function to provide the distance between two points

    IMPLICIT NONE

    REAL(r_k), DIMENSION(2), INTENT(in)                  :: pointA, pointB

!!!!!!! Variables
! pointA, B: couple of points to compute the distance between them

    fname = 'distanceRK'

    distanceRK = SQRT( (pointB(1)-pointA(1))**2 + (pointB(2)-pointA(2))**2 )

  END FUNCTION distanceRK

  REAL(r_k) FUNCTION shoelace_area_polygon(Nvertex, poly)
  ! Computing the area of a polygon using sholace formula
  ! FROM: https://en.wikipedia.org/wiki/Shoelace_formula

    IMPLICIT NONE

      INTEGER, INTENT(in)                                :: Nvertex
      REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)        :: poly

! Local
      INTEGER                                            :: i
      REAL(r_k)                                          :: areapos, areaneg

!!!!!!! Variables
! Nvertex: number of vertices of the polygon
! poly: coordinates of the vertex of the polygon (sorted)

    fname = 'shoelace_area_polygon'

    areapos = 0.
    areaneg = 0.

    DO i=1, Nvertex-1
      areapos = areapos + poly(i,1)*poly(i+1,2)
      areaneg = areaneg + poly(i+1,1)*poly(i,2)
    END DO

    areapos = areapos + poly(Nvertex,1)*poly(1,2)
    areaneg = areaneg + poly(1,1)*poly(Nvertex,2)

    shoelace_area_polygon = 0.5*(areapos - areaneg)

  END FUNCTION shoelace_area_polygon

  SUBROUTINE intersection_2Dlines(lineA, lineB, intersect, ptintersect)
  ! Subroutine to provide the intersection point between two lines on the plane using Cramer's method

    IMPLICIT NONE

    REAL(r_k), DIMENSION(2,2), INTENT(in)                  :: lineA, lineB
    LOGICAL, INTENT(out)                                   :: intersect
    REAL(r_k), DIMENSION(2), INTENT(out)                   :: ptintersect

! Local
    REAL(r_k), DIMENSION(2)                                :: segmentA, segmentB
    REAL(r_k)                                              :: a11, a12, a21, a22, z1, z2
    REAL(r_k)                                              :: det, detX, detY
    LOGICAL                                                :: axisAx, axisBx, axisAy, axisBy

!!!!!!! Variables
! lineA: couple of coordinates for the line A
! lineB: couple of coordinates for the line B
! intersect: whether two lines intersect
! ptintersect: point of intersection [(0,0) if they do not intersect]

    fname = 'intersection_2Dlines'

    axisAx = .FALSE.
    axisAy = .FALSE.
    axisBx = .FALSE.
    axisBy = .FALSE.
    ! Setting segment parameters y = A + B*x
    IF (lineA(2,1) /= lineA(1,1)) THEN
      segmentA(2) = (lineA(2,2)-lineA(1,2))/(lineA(2,1)-lineA(1,1))
      ! This might be to ask too much... ?
      !IF ( (lineA(1,1)*segmentA(2) - lineA(1,2)) /= (lineA(2,1)*segmentA(2) - lineA(2,2)) ) THEN
      !  PRINT *,'A = y1 - x2*B = ', lineA(1,2) - lineA(1,1)*segmentA(2)
      !  PRINT *,'A = y2 - x2*B = ', lineA(2,2) - lineA(2,1)*segmentA(2)
      !  msg = 'Wrong calculation of parameter A, for lineA'
      !  CALL ErrMSg(msg, fname, -1)
      !END IF
      segmentA(1) = lineA(1,2) - lineA(1,1)*segmentA(2)
      a11 = segmentA(2)
      a12 = -oneRK
      z1 = -segmentA(1)
      IF (lineA(2,2) == lineA(1,2)) axisAx = .TRUE.
    ELSE
      ! lineA || y-axis
      axisAy = .TRUE.
    END IF

    IF (lineB(2,1) /= lineB(1,1)) THEN
      segmentB(2) = (lineB(2,2)-lineB(1,2))/(lineB(2,1)-lineB(1,1))
      ! This might be to ask too much... ?
      !IF ( (lineB(1,1)*segmentB(2) - lineB(1,2)) /= (lineB(2,1)*segmentB(2) - lineB(2,2)) ) THEN
      !  PRINT *,'A = x1*B -y1 = ', lineB(1,1)*segmentB(2) - lineB(1,2)
      !  PRINT *,'A = x2*B -y2 = ', lineB(2,1)*segmentB(2) - lineB(2,2)
      !  msg = 'Wrong calculation of parameter A, for lineB'
      !  CALL ErrMSg(msg, fname, -1)
      !END IF
      segmentB(1) = lineB(1,2) - lineB(1,1)*segmentB(2)
      a21 = segmentB(2)
      a22 = -oneRK
      z2 = -segmentB(1)
      IF (lineB(2,2) == lineB(1,2)) axisBx = .TRUE.
    ELSE
      ! lineB || y-axis
      axisBy = .TRUE.
    END IF
    ! Cramer's method
    ! a11 = B1; a12 = -1
    ! a21 = B2; a22 = -1
    ! z1 = -A1
    ! z2 = -A2
    ! (a11 a12)(x) (z1)
    ! (a21 a22)(y) (z2)
    ! -------- ------ ----- ---- --- -- -
    ! det = (a11*a22-a12*a21)
    ! detX = (z1*a22-z2*a21)
    ! detY = (a11*z1-a12*z2)
    ! ptintercept = (detX/det, detY/det)

    ! Cases when some of the lines are parallel to any given axis
!    PRINT *,'          axisAx', axisAx, 'axisAy', axisAy, 'axisBx', axisBx, 'axisBy', axisBy
    IF (axisAx .OR. axisAy .OR. axisBx .OR. axisBy) THEN
      IF (axisAx) THEN
        IF (axisBy) THEN
          intersect = .TRUE.
          ptintersect(1) = lineB(1,1)
          ptintersect(2) = lineA(1,2)
        ELSE
          intersect = .TRUE.
          ptintersect(1) = (lineA(1,2)-segmentB(1))/segmentB(2)
          ptintersect(2) = lineA(1,2)
        END IF
      END IF
      IF (axisAy) THEN
        IF (axisBy) THEN
          intersect = .TRUE.
          ptintersect(1) = lineA(1,1)
          ptintersect(2) = lineB(1,2)
        ELSE
          intersect = .TRUE.
          ptintersect(1) = lineA(1,1)
          ptintersect(2) = segmentB(1) + lineA(1,1)*segmentB(2)
        END IF
      END IF
      IF (axisBx) THEN
        IF (axisAy) THEN
          intersect = .TRUE.
          ptintersect(1) = lineA(1,1)
          ptintersect(2) = lineB(1,2)
        ELSE
          intersect = .TRUE.
          ptintersect(1) = (lineB(1,2)-segmentA(1))/segmentA(2)
          ptintersect(2) = lineB(1,2)
        END IF
      END IF
      IF (axisBy) THEN
        IF (axisAx) THEN
          intersect = .TRUE.
          ptintersect(1) = lineB(1,1)
          ptintersect(2) = lineA(1,2)
        ELSE
          intersect = .TRUE.
          ptintersect(1) = lineB(1,1)
          ptintersect(2) = segmentA(1) + lineB(1,1)*segmentA(2)
        END IF
      END IF
    ELSE
      det = (a11*a22-a12*a21)
      ! Parallel lines !
      IF (det == zeroRK) THEN
        intersect = .FALSE.
        ptintersect = zeroRK
      ELSE
        intersect = .TRUE.
        detX = (z1*a22-z2*a12)
        detY = (a11*z2-a21*z1)

        ptintersect(1) = detX/det
        ptintersect(2) = detY/det
      END IF
    END IF

  END SUBROUTINE intersection_2Dlines

!refs: 
!https://www.mathopenref.com/heronsformula.html
!https://math.stackexchange.com/questions/1406340/intersect-area-of-two-polygons-in-cartesian-plan
!http://www.cap-lore.com/MathPhys/IP/
!http://www.cap-lore.com/MathPhys/IP/IL.html
!https://www.element84.com/blog/determining-the-winding-of-a-polygon-given-as-a-set-of-ordered-points
!https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order
!https://en.wikipedia.org/wiki/Shoelace_formula
!https://en.wikipedia.org/wiki/Winding_number
!https://en.wikipedia.org/wiki/Simple_polygon
!https://en.wikipedia.org/wiki/Polygon#Properties
!https://en.wikipedia.org/wiki/Convex_polygon
!https://en.wikipedia.org/wiki/Jordan_curve_theorem
!https://www.sangakoo.com/ca/temes/metode-de-cramer
!https://www.geogebra.org/m/pw4QHFYT

  SUBROUTINE intersectfaces(faceA, faceB, intersect, intersectpt)
  ! Subroutine to provide if two faces of two polygons intersect
  ! AFTER: http://www.cap-lore.com/MathPhys/IP/IL.html
  !   A: faceA(1,:)
  !   B: faceA(2,:)
  !   C: faceB(1,:)
  !   D: faceB(2,:)

    IMPLICIT NONE

    REAL(r_k), DIMENSION(2,2), INTENT(in)                :: faceA, faceB
    INTEGER, INTENT(out)                                 :: intersect
    REAL(r_k), DIMENSION(2), INTENT(out)                 :: intersectpt

! Local
    REAL(r_k)                                            :: Axmin, Aymin, Axmax, Aymax
    REAL(r_k)                                            :: Bxmin, Bymin, Bxmax, Bymax
    REAL(r_k)                                            :: areaABD, areaACD, areaBDC, areaDAB
    REAL(r_k), DIMENSION(3,2)                            :: triangle
    LOGICAL                                              :: Lintersect

!!!!!!! Variables
! faceA/B: coordinates of faces A and B to determine if they intersect
! intersect: integer to say if they intersect (0, no-intersect, +/-1 intersect)
! intersectpt: point where faces intersect [(0,0) otherwise]

    fname = 'intersectfaces'

!    PRINT *,'     ' // TRIM(fname) // ' ________'
!    PRINT *,'            faceA:', faceA(1,:), ';',faceA(2,:)
!    PRINT *,'            faceB:', faceB(1,:), ';',faceB(2,:)

    Axmin = MINVAL(faceA(:,1))
    Axmax = MAXVAL(faceA(:,1))
    Aymin = MINVAL(faceA(:,2))
    Aymax = MAXVAL(faceA(:,2))
    Bxmin = MINVAL(faceB(:,1))
    Bxmax = MAXVAL(faceB(:,1))
    Bymin = MINVAL(faceB(:,2))
    Bymax = MAXVAL(faceB(:,2))

    ! No intersection
    IF ( (Axmax <= Bxmin) .OR. (Axmin >= Bxmax) .OR. (Aymax <= Bymin) .OR. (Aymin >= Bymax) ) THEN
      intersect = 0
      intersectpt = zeroRK
    ELSE
      ! Triangle ABD
      triangle(1,:) = faceA(1,:)
      triangle(2,:) = faceA(2,:)
      triangle(3,:) = faceB(2,:)
      areaABD = shoelace_area_polygon(3, triangle)
  
      ! Triangle ACD
      triangle(1,:) = faceA(1,:)
      triangle(2,:) = faceB(1,:)
      triangle(3,:) = faceB(2,:)
      areaACD = shoelace_area_polygon(3, triangle)

      ! Triangle BDC
      triangle(1,:) = faceA(2,:)
      triangle(2,:) = faceB(2,:)
      triangle(3,:) = faceB(1,:)
      areaBDC = shoelace_area_polygon(3, triangle)

      ! Triangle DAB
      triangle(1,:) = faceB(2,:)
      triangle(2,:) = faceA(1,:)
      triangle(3,:) = faceA(2,:)
      areaDAB = shoelace_area_polygon(3, triangle)

      IF (areaABD>zeroRK .AND. areaACD>zeroRK .AND. areaBDC>zeroRK .AND. areaDAB>zeroRK) THEN
        intersect = INT(ABS(areaABD)/areaABD)
        CALL intersection_2Dlines(faceA, faceB, Lintersect, intersectpt)
      ELSE IF (areaABD<zeroRK .AND. areaACD<zeroRK .AND. areaBDC<zeroRK .AND. areaDAB<zeroRK) THEN
        intersect = INT(ABS(areaABD)/areaABD)
        CALL intersection_2Dlines(faceA, faceB, Lintersect, intersectpt)
      ELSE
        intersect = 0
        intersectpt = zeroRK
      END IF
!      PRINT *,'     intersect faces: areaABD',areaABD, 'areaACD', areaACD, 'areaBDC',areaBDC, 'areaDAB',areaDAB, 'prod', &
!        areaABD*areaACD*areaBDC*areaDAB, 'L:', areaABD*areaACD*areaBDC*areaDAB > zeroRK, 'I', intersect

    END IF

  END SUBROUTINE intersectfaces

  LOGICAL FUNCTION poly_has_point(Nvertex, polygon, point)
  ! Function to determine if a polygon has already a given point as one of its vertex

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: Nvertex
    REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)          :: polygon
    REAL(r_k), DIMENSION(2), INTENT(in)                  :: point

! Local
    INTEGER                                              :: iv
    REAL(r_k), DIMENSION(2)                              :: diff

!!!!!!! Vertrex
! Nvertex: number of vertexs of the polygon
! polygon: vertexs of the polygon
! point: point to look for its ownership into the polygon

    fname = 'poly_has_point'

    poly_has_point = .FALSE.
    DO iv=1, Nvertex
      diff = polygon(iv,:)-point
      IF ( (diff(1) == zeroRK) .AND. (diff(2) == zeroRK)) THEN
        poly_has_point = .TRUE.
        EXIT
      END IF
    END DO

  END FUNCTION poly_has_point

  SUBROUTINE join_polygon(NvertexA, NvertexB, NvertexAB, polyA, polyB, Ncoinvertex, coinpoly)
  ! Subroutine to join two polygons
  ! AFTER: http://www.cap-lore.com/MathPhys/IP/ and http://www.cap-lore.com/MathPhys/IP/IL.html

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: NvertexA, NvertexB, NvertexAB
    REAL(r_k), DIMENSION(NvertexA,2), INTENT(in)         :: polyA
    REAL(r_k), DIMENSION(NvertexB,2), INTENT(in)         :: polyB
    INTEGER, INTENT(out)                                 :: Ncoinvertex
    REAL(r_k), DIMENSION(NvertexAB,2), INTENT(out)        :: coinpoly

! Local
    INTEGER                                              :: iA, iB, icoin, ii
    REAL(r_k), DIMENSION(2,2)                            :: face1, face2
    INTEGER                                              :: intersct
    REAL(r_k), DIMENSION(2)                              :: ptintersct


!!!!!!! variables
! NvertexA: number of vertexs polygon A
! NvertexB: number of vertexs polygon B
! polyA: pairs of coordinates for the polygon A (clockwise)
! polyB: pairs of coordinates for the polygon B (clockwise)
! Ncoinvertex: number of vertexes for the coincident polygon
! coinpoly: pairs of coordinates for the coincident polygon (clockwise)

    fname = 'join_polygon'

    icoin = 0
    coinpoly = 0.

    ! First, include that vertex which do not lay within any polygon
    DO iA=1, NvertexA
      !PRINT *, '  iA:', iA, ':', polyA(iA,:), point_inside(polyA(iA,:), NvertexB, polyB)
      IF (.NOT. point_inside(polyA(iA,:), NvertexB, polyB)) THEN
        icoin = icoin + 1
        coinpoly(icoin,:) = polyA(iA,:)
      END IF
    END DO

    DO iB=1, NvertexB
      !PRINT *, '  iB:', iB, ':', polyB(iB,:), point_inside(polyB(iB,:), NvertexA, polyA)
      IF (.NOT. point_inside(polyB(iB,:), NvertexA, polyA)) THEN
        icoin = icoin + 1
        coinpoly(icoin,:) = polyB(iB,:)
      END IF
    END DO

    DO iA=1, NvertexA
      ! Getting couple of vertexs from polyA and polyB
      IF (iA /= NvertexA) THEN
        face1(1,:) = polyA(iA,:)
        face1(2,:) = polyA(iA+1,:)
      ELSE
        face1(1,:) = polyA(iA,:)
        face1(2,:) = polyA(1,:)
      END IF
      DO iB=1, NvertexB
        IF (iB /= NvertexB) THEN
          face2(1,:) = polyB(iB,:)
          face2(2,:) = polyB(iB+1,:)
        ELSE
          face2(1,:) = polyB(iB,:)
          face2(2,:) = polyB(1,:)
        END IF
        
        ! Compute areas of the four possible triangles. Introduce the coincident vertexs not included
        CALL intersectfaces(face1, face2, intersct, ptintersct)
        !PRINT *,iA,':',face1(1,:),';',face1(2,:), '=', iB, face2(1,:),';',face2(2,:), '<>', intersct,':', ptintersct
        IF (intersct == 1) THEN
          IF (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct) ) THEN
            icoin = icoin + 1
            coinpoly(icoin,:) = ptintersct
          END IF
        ELSE IF (intersct == -1) THEN
          IF (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct) ) THEN
            icoin = icoin + 1
            coinpoly(icoin,:) = ptintersct
          END IF
        END IF

      END DO
    END DO
    Ncoinvertex = icoin

  END SUBROUTINE join_polygon

  SUBROUTINE sort_polygon(Nvertex, polygon, sense, Nnewvertex, newpoly)
  ! Subroutine to sort a polygon using its center as average of the coordinates and remove duplicates
  !    Should be used the centroid instead, but by now let do it simple
  !    https://en.wikipedia.org/wiki/Centroid


    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: Nvertex, sense
    REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)          :: polygon
    INTEGER, INTENT(out)                                 :: Nnewvertex
    REAL(r_k), DIMENSION(Nvertex,2), INTENT(out)         :: newpoly

! Local
    INTEGER                                              :: iv, j
    REAL(r_k)                                            :: vang
    REAL(r_k), DIMENSION(2)                              :: center
    REAL(r_k), DIMENSION(Nvertex)                        :: angles
    REAL(r_k), DIMENSION(Nvertex,2)                      :: sortpoly

!!!!!!! Variables
! Nvertex: number of vertices
! polygon: coordinates of the vertices of the polygon
! sense: sens of sorting thepolygon (1: clockwise, -1: anti-clockwise)
! sortpoly: sorted polygon
! Nnewvertex: number of vertices new polygon
! newpoly: sorted and duplicate removed polygon

    fname = 'sort_polygon'

    ! To be substituted by centroid calculation (which requires already sorted vetexs...)
    center(1) = SUM(polygon(:,1))/Nvertex
    center(2) = SUM(polygon(:,2))/Nvertex

    DO iv=1, Nvertex
      angles(iv) = ATAN2(polygon(iv,2)-center(2),polygon(iv,1)-center(1))
    END DO
    CALL sortR_K(angles, Nvertex)

    sortpoly = zeroRK
    DO iv=1, Nvertex
      DO j=1, Nvertex
        vang = ATAN2(polygon(j,2)-center(2), polygon(j,1)-center(1))
        IF (angles(iv) == vang) THEN
          IF (sense == -1) THEN
            sortpoly(iv,:) = polygon(j,:)
          ELSE
            sortpoly(Nvertex-iv+1,:) = polygon(j,:)
          END IF
          EXIT
        END IF
      END DO
    END DO

    newpoly(1,:) = sortpoly(1,:)
    j = 1
    DO iv=2, Nvertex
      IF (.NOT.poly_has_point(j,newpoly(1:j,:),sortpoly(iv,:)) ) THEN
        j = j+1
        newpoly(j,:) = sortpoly(iv,:)
      END IF
    END DO
    Nnewvertex = j

  END SUBROUTINE sort_polygon

  LOGICAL FUNCTION point_inside(point, Nvertex, polygon)
  ! Function to determine if a given point is inside a polygon providing its sorted vertices
  ! FROM: https://en.wikipedia.org/wiki/Point_in_polygon

    IMPLICIT NONE

    REAL(r_k), DIMENSION(2), INTENT(in)                  :: point
    INTEGER, INTENT(in)                                  :: Nvertex
    REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)          :: polygon

    ! Local
    INTEGER                                              :: iv, Nintersect
    INTEGER                                              :: cross
    REAL(r_k)                                            :: xmin
    REAL(r_k), DIMENSION(2)                              :: crosspoint
    REAL(r_k), DIMENSION(2,2)                            :: face1, face2
    REAL(r_k), DIMENSION(Nvertex)                        :: abovebelow

!!!!!!! Variables
! point: point to look for
! Nvertrex: number of vertices of a polygon
! polygon: vertices of a polygon

    fname = 'point_inside'

    xmin = MINVAL(polygon(:,1))

    ! Looking for the intersection with the ray
    Nintersect = 0
    face1(1,:) = (/ xmin-0.5, point(2) /)
    face1(2,:) = (/ point(1), point(2) /)

    DO iv = 1, Nvertex
      IF (iv /= Nvertex) THEN
        face2(1,:) = polygon(iv,:)
        face2(2,:) = polygon(iv+1,:)
      ELSE
        face2(1,:) = polygon(iv,:)
        face2(2,:) = polygon(1,:)
      END IF
      CALL intersectfaces(face1, face2, cross, crosspoint)
      IF (cross /= 0) THEN
        Nintersect = Nintersect + 1
        abovebelow(Nintersect) = iv
      END IF
    END DO

    IF (MOD(Nintersect,2) == 0) THEN
      point_inside = .FALSE.
    ELSE
      point_inside = .TRUE.
    END IF

  END FUNCTION point_inside

  LOGICAL FUNCTION point_in_face(pt, Nvertex, poly)
  ! Function to determine if a given point is on a face of a polygon

    IMPLICIT NONE

    REAL(r_k), DIMENSION(2), INTENT(in)                  :: pt
    INTEGER, INTENT(in)                                  :: Nvertex
    REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)          :: poly
! Local
    INTEGER                                              :: iv
    REAL(r_k)                                            :: ix, ex, iy, ey, tmpv
    REAL(r_k)                                            :: dx, dy, A, B

!!!!!!! Variables
! pt: point to look for
! Nvertex: Number of vertices of the polygon
! poly: polygon
    fname = 'point_in_face'

    point_in_face = .FALSE.
    DO iv=1, Nvertex
      IF (iv < Nvertex) THEN
        ix = poly(iv,1)
        ex = poly(iv+1,1)
        iy = poly(iv,2)
        ey = poly(iv+1,2)
      ELSE
        ix = poly(iv,1)
        ex = poly(1,1)
        iy = poly(iv,2)
        ey = poly(1,2)
      END IF    
      dx = ex - ix
      dy = ey - iy

      IF (dx == zeroRK) THEN
        IF (pt(1) == ix) THEN
          IF ( (iy < ey) .AND. (pt(2) >= iy) .AND. pt(2) <= ey) THEN
            point_in_face = .TRUE.
            EXIT
          ELSE IF ( (iy > ey) .AND. (pt(2) >= ey) .AND. pt(2) <= iy) THEN
            point_in_face = .TRUE.
            EXIT
          END IF
        END IF
      ELSE
        IF (dy == zeroRK) THEN
          IF (pt(2) == iy) THEN
            IF ((ix < ex) .AND. (pt(1) >= ix) .AND. pt(1) <= ex) THEN
              point_in_face = .TRUE.
              EXIT            
            ELSE IF ((ix > ex) .AND. (pt(1) >= ex) .AND. pt(1) <= ix) THEN
              point_in_face = .TRUE.
              EXIT
            END IF
          END IF
        ELSE
          A = iy
          B = (ey-iy)/(ex-ix)
          IF (A+B*(pt(1)-ix) == pt(2)) THEN
            point_in_face = .TRUE.
            EXIT
          END IF
        END IF
      END IF
    END DO

  END FUNCTION point_in_face

  SUBROUTINE coincident_polygon(NvertexA, NvertexB, NvertexAB, polyA, polyB, Ncoinvertex, coinpoly)
  ! Subroutine to provide the intersection polygon between two polygons
  ! AFTER: http://www.cap-lore.com/MathPhys/IP/ and http://www.cap-lore.com/MathPhys/IP/IL.html

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: NvertexA, NvertexB, NvertexAB
    REAL(r_k), DIMENSION(NvertexA,2), INTENT(in)         :: polyA
    REAL(r_k), DIMENSION(NvertexB,2), INTENT(in)         :: polyB
    INTEGER, INTENT(out)                                 :: Ncoinvertex
    REAL(r_k), DIMENSION(NvertexAB,2), INTENT(out)       :: coinpoly

! Local
    INTEGER                                              :: iA, iB, icoin, ii
    REAL(r_k), DIMENSION(2,2)                            :: face1, face2
    INTEGER                                              :: intersct
    REAL(r_k), DIMENSION(2)                              :: ptintersct

!!!!!!! variables
! NvertexA: number of vertexs polygon A
! NvertexB: number of vertexs polygon B
! polyA: pairs of coordinates for the polygon A (clockwise)
! polyB: pairs of coordinates for the polygon B (clockwise)
! Ncoinvertex: number of vertexes for the coincident polygon
! coinpoly: pairs of coordinates for the coincident polygon (clockwise)

    fname = 'coincident_polygon'

    icoin = 0
    coinpoly = 0.
    ! First, include that vertex which lay within any polygon
    DO iA=1, NvertexA
      IF (point_inside(polyA(iA,:), NvertexB, polyB)) THEN
        icoin = icoin + 1
        coinpoly(icoin,:) = polyA(iA,:)
      END IF
      IF (point_in_face(polyA(iA,:), NvertexB, polyB)) THEN
        icoin = icoin + 1
        coinpoly(icoin,:) = polyA(iA,:)
      END IF
    END DO

    DO iB=1, NvertexB
      IF (point_inside(polyB(iB,:), NvertexA, polyA)) THEN
        icoin = icoin + 1
        coinpoly(icoin,:) = polyB(iB,:)
      END IF
      IF (point_in_face(polyB(iB,:), NvertexA, polyA)) THEN
        icoin = icoin + 1
        coinpoly(icoin,:) = polyB(iB,:)
      END IF
    END DO

    ! Look interesections
    DO iA=1, NvertexA
      ! Getting couple of vertexs from polyA and polyB
      IF (iA /= NvertexA) THEN
        face1(1,:) = polyA(iA,:)
        face1(2,:) = polyA(iA+1,:)
      ELSE
        face1(1,:) = polyA(iA,:)
        face1(2,:) = polyA(1,:)
      END IF
      DO iB=1, NvertexB
        IF (iB /= NvertexB) THEN
          face2(1,:) = polyB(iB,:)
          face2(2,:) = polyB(iB+1,:)
        ELSE
          face2(1,:) = polyB(iB,:)
          face2(2,:) = polyB(1,:)
        END IF
        
        ! Compute areas of the four possible triangles. Introduce the coincident vertexs not included
        CALL intersectfaces(face1, face2, intersct, ptintersct)
        !PRINT *,iA,':',face1(1,:),';',face1(2,:), '=', iB, face2(1,:),';',face2(2,:), '<>', intersct,':', ptintersct
        IF ((intersct /= 0) .AND. (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct)) ) THEN
          icoin = icoin + 1
          coinpoly(icoin,:) = ptintersct
        END IF

      END DO
    END DO
    Ncoinvertex = icoin

  END SUBROUTINE coincident_polygon

  SUBROUTINE grid_within_polygon(NvertexA, polygonA, dx, dy, dxy, xCvals, yCvals, Nvertexmax, xBvals, &
    yBvals, Ngridsin, gridsin)
  ! Subroutine to determine which grid cells from a matrix lay inside a polygon

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: NvertexA, dx, dy, dxy, Nvertexmax
    REAL(r_k), DIMENSION(NvertexA,2), INTENT(in)         :: polygonA
    REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: xCvals, yCvals
    REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in)   :: xBvals, yBvals
    INTEGER, INTENT(out)                                 :: Ngridsin
    INTEGER, DIMENSION(dxy,2), INTENT(out)               :: gridsin

! Local
    INTEGER                                              :: ix, iy, iv
    REAL(r_k), DIMENSION(2)                              :: centergrid, vertex
    LOGICAL, DIMENSION(dx,dy)                            :: within

!!!!!!! Variables
! NvertexA: Number of vertices of the polygin to find the grids
! polygonA: ordered vertices of the polygon
! dx, dy: shape of the matrix with the grid points
! xCvals, yCvals: coordinates of the center of the grid cells
! Nvertexmax: Maximum number of vertices of the grid cells
! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex)
! Ngridsin: number of grids with some extension within the polygon
! gridsin: grids within the polygin
! percentages: percentages of area of each of the grids within the polygon

    fname = 'spacepercen_within_reg'

    Ngridsin = 0
    gridsin = 0
    within = .FALSE.
    DO ix = 1, dx
      DO iy = 1, dy
        IF (.NOT.within(ix,iy)) THEN
          centergrid = (/ xCvals(ix,iy), yCvals(ix,iy) /)
          ! By grid center
          IF (point_inside(centergrid, NvertexA, polygonA)) THEN
            Ngridsin = Ngridsin + 1
            ! Getting coordinates
            gridsin(Ngridsin,1) = ix
            gridsin(Ngridsin,2) = iy
            within(ix,iy) = .TRUE.
            CYCLE
          END IF

          ! Getting grid vertices
          DO iv=1, Nvertexmax
            IF (.NOT.within(ix,iy)) THEN
              IF (xBvals(ix,iy,iv) /= fillvalI) THEN
                vertex = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /)
                IF (point_inside(vertex, NvertexA, polygonA)) THEN
                  Ngridsin = Ngridsin + 1
                  ! Getting coordinates
                  gridsin(Ngridsin,1) = ix
                  gridsin(Ngridsin,2) = iy
                  within(ix,iy) = .TRUE.
                  CYCLE
                END IF
              END IF
            END IF
          END DO

        END IF
      END DO
    END DO

  END SUBROUTINE grid_within_polygon

  SUBROUTINE spacepercen_within_reg(NvertexA, polygonA, dx, dy, Nvertexmax, xBvals, yBvals,           &
    Ngridsin, gridsin, strict, percentages)
  ! Subroutine to compute the percentage of a series of grid cells which are encompassed by a polygon
  ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x 
  !   and y axes. 

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: NvertexA, dx, dy, Nvertexmax
    REAL(r_k), DIMENSION(NvertexA,2), INTENT(in)         :: polygonA
    REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in)   :: xBvals, yBvals
    INTEGER, INTENT(in)                                  :: Ngridsin
    INTEGER, DIMENSION(Ngridsin,2), INTENT(in)           :: gridsin
    LOGICAL, INTENT(in)                                  :: strict
    REAL(r_k), DIMENSION(Ngridsin), INTENT(out)          :: percentages

! Local
   INTEGER                                               :: ig, iv, ix, iy
   INTEGER                                               :: Nvertex, NvertexAgrid, Ncoin, Nsort
   CHARACTER(len=20)                                     :: DS
   REAL(r_k)                                             :: areapoly, areagpoly, totarea, totpercent
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                :: vertexgrid, icoinpoly, coinpoly,          &
     sortpoly, poly

!!!!!!! Variables
! NvertexA: Number of vertices of the polygin to find the grids
! polygonA: ordered vertices of the polygon
! dx, dy: shape of the matrix with the grid points
! xCvals, yCvals: coordinates of the center of the grid cells
! Nvertexmax: Maximum number of vertices of the grid cells
! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex)
! Ngridsin: number of grids with some extension within the polygon
! gridsin: grids within the polygon
! strict: give an error if the area of the polygon is not fully covered
! percentages: percentages of area of each of the grids within the polygon

    fname = 'spacepercen_within_reg'

    percentages = zeroRK
    totpercent = zeroRK
    totarea = zeroRK

    areapoly = shoelace_area_polygon(NvertexA, polygonA)

    DO ig = 1, Ngridsin
      ix = gridsin(ig,1)
      iy = gridsin(ig,2)

      ! Getting grid vertices
      Nvertex = 0
      DO iv=1, Nvertexmax
        IF (xBvals(ix,iy,iv) /= fillvalI) THEN
          Nvertex = Nvertex + 1
        END IF
      END DO
      IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
      ALLOCATE(vertexgrid(Nvertex,2))
      vertexgrid(:,1) = xBvals(ix,iy,1:Nvertex)
      vertexgrid(:,2) = yBvals(ix,iy,1:Nvertex)

      ! Getting common vertices
      NvertexAgrid = NvertexA*Nvertex*2
      IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
      ALLOCATE(icoinpoly(NvertexAgrid,2))
      CALL coincident_polygon(NvertexA, Nvertex, NvertexAgrid, polygonA, vertexgrid, Ncoin, icoinpoly)

      IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly)
      ALLOCATE(coinpoly(Ncoin,2))
      DO iv=1, Ncoin
        coinpoly(iv,:) = icoinpoly(iv,:)
      END DO

      IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly)
      ALLOCATE(sortpoly(Ncoin,2))
      CALL sort_polygon(Ncoin, coinpoly, 1, Nsort, sortpoly)

      IF (ALLOCATED(poly)) DEALLOCATE(poly)
      ALLOCATE(poly(Nsort,2))
      DO iv=1, Nsort
        poly(iv,:) = sortpoly(iv,:)
      END DO

      areagpoly = shoelace_area_polygon(Nsort, poly)
      IF (INT(LOG10(EPSILON(totpercent))) < 12) THEN
        totarea = totarea + ABS(areagpoly)
        percentages(ig) = ABS(areagpoly / areapoly)
! f2py does not like it!
!      ELSE
!        totarea = totarea + DABS(areagpoly)
!        percentages(ig) = DABS(areagpoly / areapoly)
      END IF
      totpercent = totpercent + percentages(ig)
    END DO

    IF (INT(LOG10(EPSILON(totpercent))) < 12) THEN
      IF (strict .AND. ABS(totpercent - oneRK) > epsilonRK) THEN
        PRINT *, 'totarea:', totarea, ' area polygon:', areapoly
        PRINT *, 'totpercent:', totpercent, ' oneRK:', oneRK, ' diff:', totpercent - oneRK
        WRITE(DS,'(F20.8)')ABS(totpercent - oneRK)
        msg = 'sum of all grid space percentages does not cover (' // TRIM(DS) // ') all polygon'
        CALL ErrMsg(msg, fname, -1)
      END IF
    ELSE
! f2py does not like it!
!      IF (strict .AND. ABS(totpercent - oneRK) > epsilonRK) THEN
!        PRINT *, 'totarea:', totarea, ' area polygon:', areapoly
!        PRINT *, 'totpercent:', totpercent, ' oneRK:', oneRK, ' diff:', totpercent - oneRK
!        WRITE(DS,'(F20.16)')ABS(totpercent - oneRK)
!        msg = 'sum of all grid space percentages does not cover (' // TRIM(DS) // ') all polygon'
!        CALL ErrMsg(msg, fname, -1)
!      END IF
    END IF

    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
    IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
    IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly)
    IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly)
    IF (ALLOCATED(poly)) DEALLOCATE(poly)

  END SUBROUTINE spacepercen_within_reg

  SUBROUTINE grid_spacepercen_within_reg(NvertexA, polygonA, dx, dy, Nvertexmax, xBvals, yBvals,      &
    Ngridsin, gridsin, strict, gridspace, percentages)
  ! Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed 
  !   by a polygon
  ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x 
  !   and y axes. 

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: NvertexA, dx, dy, Nvertexmax
    REAL(r_k), DIMENSION(NvertexA,2), INTENT(in)         :: polygonA
    REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in)   :: xBvals, yBvals
    INTEGER, INTENT(in)                                  :: Ngridsin
    INTEGER, DIMENSION(Ngridsin,2), INTENT(in)           :: gridsin
    LOGICAL, INTENT(in)                                  :: strict
    REAL(r_k), DIMENSION(Ngridsin), INTENT(out)          :: gridspace, percentages

! Local
   INTEGER                                               :: ig, iv, ix, iy
   INTEGER                                               :: Nvertex, NvertexAgrid, Ncoin, Nsort
   CHARACTER(len=20)                                     :: DS
   REAL(r_k)                                             :: areapoly, areagpoly
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                :: vertexgrid, icoinpoly, coinpoly,          &
     sortpoly, poly

!!!!!!! Variables
! NvertexA: Number of vertices of the polygon to find the grids
! polygonA: ordered vertices of the polygon
! dx, dy: shape of the matrix with the grid points
! xCvals, yCvals: coordinates of the center of the grid cells
! Nvertexmax: Maximum number of vertices of the grid cells
! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex)
! Ngridsin: number of grids with some extension within the polygon
! gridsin: grids within the polygon
! strict: give an error if the area of the polygon is not fully covered
! gridspace: area of each of the grids
! percentages: percentages of grid area of each of the grids within the polygon

    fname = 'grid_spacepercen_within_reg'

    gridspace = zeroRK
    percentages = zeroRK

    DO ig = 1, Ngridsin
      ix = gridsin(ig,1)
      iy = gridsin(ig,2) 

     ! Getting grid vertices
      Nvertex = 0
      DO iv=1, Nvertexmax
        IF (xBvals(ix,iy,iv) /= fillvalI) THEN
          Nvertex = Nvertex + 1
        END IF
      END DO
      IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
      ALLOCATE(vertexgrid(Nvertex,2))
      vertexgrid(:,1) = xBvals(ix,iy,1:Nvertex)
      vertexgrid(:,2) = yBvals(ix,iy,1:Nvertex)
      areapoly = shoelace_area_polygon(Nvertex, vertexgrid)

      ! Getting common vertices
      NvertexAgrid = NvertexA*Nvertex*2
      IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
      ALLOCATE(icoinpoly(NvertexAgrid,2))
      CALL coincident_polygon(NvertexA, Nvertex, NvertexAgrid, polygonA, vertexgrid, Ncoin, icoinpoly)

      IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly)
      ALLOCATE(coinpoly(Ncoin,2))
      DO iv=1, Ncoin
        coinpoly(iv,:) = icoinpoly(iv,:)
      END DO

      IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly)
      ALLOCATE(sortpoly(Ncoin,2))
      CALL sort_polygon(Ncoin, coinpoly, 1, Nsort, sortpoly)

      IF (ALLOCATED(poly)) DEALLOCATE(poly)
      ALLOCATE(poly(Nsort,2))
      DO iv=1, Nsort
        poly(iv,:) = sortpoly(iv,:)
      END DO

      areagpoly = shoelace_area_polygon(Nsort, poly)
      gridspace(ig) = ABS(areapoly)
      percentages(ig) = ABS(areagpoly / areapoly)
    END DO

    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
    IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
    IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly)
    IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly)
    IF (ALLOCATED(poly)) DEALLOCATE(poly)

  END SUBROUTINE grid_spacepercen_within_reg

  SUBROUTINE grid_spacepercen_within_reg_providing_polys(NvertexA, polygonA, dx, dy, Nvertexmax,      &
    xBvals, yBvals, Ngridsin, gridsin, strict, Nmaxver2, Ncoinpoly, ccoinpoly, gridspace, percentages)
  ! Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed 
  !   by a polygon providing coordinates of the resultant polygons
  ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x 
  !   and y axes. 

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: NvertexA, dx, dy, Nvertexmax, Nmaxver2
    REAL(r_k), DIMENSION(NvertexA,2), INTENT(in)         :: polygonA
    REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in)   :: xBvals, yBvals
    INTEGER, INTENT(in)                                  :: Ngridsin
    INTEGER, DIMENSION(Ngridsin,2), INTENT(in)           :: gridsin
    LOGICAL, INTENT(in)                                  :: strict
    INTEGER, DIMENSION(Ngridsin), INTENT(out)            :: Ncoinpoly
    REAL(r_k), DIMENSION(Ngridsin,Nmaxver2,2),                                                        &
      INTENT(out)                                        :: ccoinpoly
    REAL(r_k), DIMENSION(Ngridsin), INTENT(out)          :: gridspace, percentages

! Local
   INTEGER                                               :: ig, iv, ix, iy
   INTEGER                                               :: Nvertex, NvertexAgrid, Ncoin, Nsort
   CHARACTER(len=20)                                     :: DS
   REAL(r_k)                                             :: areapoly, areagpoly
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                :: vertexgrid, icoinpoly, coinpoly,          &
     sortpoly, poly

!!!!!!! Variables
! NvertexA: Number of vertices of the polygon to find the grids
! polygonA: ordered vertices of the polygon
! dx, dy: shape of the matrix with the grid points
! xCvals, yCvals: coordinates of the center of the grid cells
! Nvertexmax: Maximum number of vertices of the grid cells
! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex)
! Ngridsin: number of grids with some extension within the polygon
! gridsin: grids within the polygon
! strict: give an error if the area of the polygon is not fully covered
! Nmaxver2: maximum possible number of vertices of the coincident polygon
! Ncoinpoly: number of vertices of the coincident polygon
! coinpoly: coordinates of the vertices of the coincident polygon
! gridspace: area of each of the grids
! percentages: percentages of grid area of each of the grids within the polygon

    fname = 'grid_spacepercen_within_reg_providing_polys'

    gridspace = zeroRK
    percentages = zeroRK

    DO ig = 1, Ngridsin
      ix = gridsin(ig,1)
      iy = gridsin(ig,2) 

     ! Getting grid vertices
      Nvertex = 0
      DO iv=1, Nvertexmax
        IF (xBvals(ix,iy,iv) /= fillvalI) THEN
          Nvertex = Nvertex + 1
        END IF
      END DO
      IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
      ALLOCATE(vertexgrid(Nvertex,2))
      vertexgrid(:,1) = xBvals(ix,iy,1:Nvertex)
      vertexgrid(:,2) = yBvals(ix,iy,1:Nvertex)
      areapoly = shoelace_area_polygon(Nvertex, vertexgrid)

      ! Getting common vertices
      NvertexAgrid = NvertexA*Nvertex*2
      IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
      ALLOCATE(icoinpoly(NvertexAgrid,2))
      CALL coincident_polygon(NvertexA, Nvertex, NvertexAgrid, polygonA, vertexgrid, Ncoin, icoinpoly)

      IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly)
      ALLOCATE(coinpoly(Ncoin,2))
      DO iv=1, Ncoin
        coinpoly(iv,:) = icoinpoly(iv,:)
      END DO

      IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly)
      ALLOCATE(sortpoly(Ncoin,2))
      CALL sort_polygon(Ncoin, coinpoly, 1, Nsort, sortpoly)

      IF (ALLOCATED(poly)) DEALLOCATE(poly)
      ALLOCATE(poly(Nsort,2))
      DO iv=1, Nsort
        poly(iv,:) = sortpoly(iv,:)
      END DO

      areagpoly = shoelace_area_polygon(Nsort, poly)
      Ncoinpoly(ig)= Nsort
      ccoinpoly(ig,:,:) = poly(:,:)
      gridspace(ig) = ABS(areapoly)
      percentages(ig) = ABS(areagpoly / areapoly)
    END DO

    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
    IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
    IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly)
    IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly)
    IF (ALLOCATED(poly)) DEALLOCATE(poly)

  END SUBROUTINE grid_spacepercen_within_reg_providing_polys

  SUBROUTINE spacepercen(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals, xBBvals, yBBvals,      &
    dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Ngridsin, gridsin, areas, percentages)
  ! Subroutine to compute the space-percentages of a series of grid cells (B) into another series of 
  !   grid-cells (A)
  ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x 
  !   and y axes. 

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dxA, dyA, NAvertexmax
    INTEGER, INTENT(in)                                  :: dxB, dyB, NBvertexmax, dxyB
    REAL(r_k), DIMENSION(dxA,dyA), INTENT(in)            :: xCAvals, yCAvals
    REAL(r_k), DIMENSION(dxB,dyB), INTENT(in)            :: xCBvals, yCBvals
    REAL(r_k), DIMENSION(dxA,dyA,NAvertexmax), INTENT(in):: xBAvals, yBAvals
    REAL(r_k), DIMENSION(dxB,dyB,NBvertexmax), INTENT(in):: xBBvals, yBBvals
    LOGICAL, INTENT(in)                                  :: strict
    INTEGER, DIMENSION(dxA,dyA), INTENT(out)             :: Ngridsin
    INTEGER, DIMENSION(dxA,dyA,dxyB,2), INTENT(out)      :: gridsin
    REAL(r_k), DIMENSION(dxA,dyA), INTENT(out)           :: areas
    REAL(r_k), DIMENSION(dxA,dyA,dxyB), INTENT(out)      :: percentages

! Local
   INTEGER                                               :: iv, ix, iy
   INTEGER                                               :: Nvertex
   INTEGER, ALLOCATABLE, DIMENSION(:,:)                  :: poinsin
   CHARACTER(len=20)                                     :: IS
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                :: vertexgrid

!!!!!!! Variables
! dxA, dyA: shape of the matrix with the grid points A
! xCAvals, yCAvals: coordinates of the center of the grid cells A
! NAvertexmax: Maximum number of vertices of the grid cells A
! xBAvals, yBAvals: coordinates of th vertices of the grid cells A (-99999 for no vertex)
! dxB, dyB: shape of the matrix with the grid points B
! xCBvals, yCBvals: coordinates of the center of the grid cells B
! NBvertexmax: Maximum number of vertices of the grid cells B
! xBBvals, yBBvals: coordinates of th vertices of the grid cells B (-99999 for no vertex)
! strict: give an error if the area of the polygon is not fully covered
! Ngridsin: number of grids from grid B with some extension within the grid cell A
! gridsin: indices of B grids within the grids of A
! areas: areas of the polygons
! percentages: percentages of area of cells B of each of the grids within the grid cell A

    fname = 'spacepercen'

    DO ix = 1, dxA
      DO iy = 1, dyA

        ! Getting grid vertices
        Nvertex = 0
        DO iv=1, NAvertexmax
          IF (xBAvals(ix,iy,iv) /= fillval64) THEN
           Nvertex = Nvertex + 1
          END IF
        END DO
        IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
        ALLOCATE(vertexgrid(Nvertex,2))
        vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex)
        vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex)
  
        CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xCBvals, yCBvals,            &
          NBvertexmax, xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,:,:))
   
        IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)
        ALLOCATE(poinsin(Ngridsin(ix,iy),2))

        DO iv=1, Ngridsin(ix,iy)
          poinsin(iv,1) = gridsin(ix,iy,iv,1)
          poinsin(iv,2) = gridsin(ix,iy,iv,2)
        END DO

        areas(ix,iy) = shoelace_area_polygon(Nvertex, vertexgrid)
        CALL spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals, yBBvals,     &
          Ngridsin(ix,iy), poinsin, strict, percentages(ix,iy,:))

      END DO
    END DO

    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
    IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)

  END SUBROUTINE spacepercen

  SUBROUTINE grid_spacepercen(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals, xBBvals, yBBvals, &
    dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Ngridsin, gridsin,  areas2D, areas,   &
    percentages)
  ! Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another 
  !   series of grid-cells (A) porviding coincident polygons
  ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x 
  !   and y axes. 

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dxA, dyA, NAvertexmax
    INTEGER, INTENT(in)                                  :: dxB, dyB, NBvertexmax, dxyB
    REAL(r_k), DIMENSION(dxA,dyA), INTENT(in)            :: xCAvals, yCAvals
    REAL(r_k), DIMENSION(dxB,dyB), INTENT(in)            :: xCBvals, yCBvals
    REAL(r_k), DIMENSION(dxA,dyA,NAvertexmax), INTENT(in):: xBAvals, yBAvals
    REAL(r_k), DIMENSION(dxB,dyB,NBvertexmax), INTENT(in):: xBBvals, yBBvals
    LOGICAL, INTENT(in)                                  :: strict
    INTEGER, DIMENSION(dxA,dyA), INTENT(out)             :: Ngridsin
    INTEGER, DIMENSION(dxA,dyA,dxyB,2), INTENT(out)      :: gridsin
    REAL(r_k), DIMENSION(dxB,dyB), INTENT(out)           :: areas2D
    REAL(r_k), DIMENSION(dxA,dyA,dxyB), INTENT(out)      :: areas,percentages

! Local
   INTEGER                                               :: iv, ix, iy
   INTEGER                                               :: Nvertex, Nptin
   INTEGER, ALLOCATABLE, DIMENSION(:,:)                  :: poinsin
   CHARACTER(len=20)                                     :: IS
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                :: vertexgrid

!!!!!!! Variables
! dxA, dyA: shape of the matrix with the grid points A
! xCAvals, yCAvals: coordinates of the center of the grid cells A
! NAvertexmax: Maximum number of vertices of the grid cells A
! xBAvals, yBAvals: coordinates of th vertices of the grid cells A (-99999 for no vertex)
! dxB, dyB: shape of the matrix with the grid points B
! xCBvals, yCBvals: coordinates of the center of the grid cells B
! NBvertexmax: Maximum number of vertices of the grid cells B
! xBBvals, yBBvals: coordinates of th vertices of the grid cells B (-99999 for no vertex)
! strict: give an error if the area of the polygon is not fully covered
! Ngridsin: number of grids from grid B with some extension within the grid cell A
! gridsin: indices of B grids within the grids of A
! areas2D: areas of the grids as 2D matrix in the original shape
! areas: areas of cells B of each of the grids inside the grid cell A
! percentages: percentages of area of cells B of each of the grids inside the grid cell A

    fname = 'grid_spacepercen'

    areas2D = zeroRK
    areas = zeroRK
    percentages = zeroRK

    DO ix = 1, dxA
      DO iy = 1, dyA

        ! Getting grid vertices
        Nvertex = 0
        DO iv=1, NAvertexmax
          IF (xBAvals(ix,iy,iv) /= fillval64) THEN
           Nvertex = Nvertex + 1
          END IF
        END DO
        IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
        ALLOCATE(vertexgrid(Nvertex,2))
        vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex)
        vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex)
  
        CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xCBvals, yCBvals,            &
          NBvertexmax, xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,1:dxyB,:))
   
        IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)
        ALLOCATE(poinsin(Ngridsin(ix,iy),2))

        DO iv=1, Ngridsin(ix,iy)
          poinsin(iv,1) = gridsin(ix,iy,iv,1)
          poinsin(iv,2) = gridsin(ix,iy,iv,2)
        END DO

        Nptin = Ngridsin(ix,iy)
        CALL grid_spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals,        &
          yBBvals, Ngridsin(ix,iy), poinsin, strict, areas(ix,iy,1:Nptin), percentages(ix,iy,1:Nptin))

        ! Filling areas
        DO iv = 1, Ngridsin(ix,iy)
          IF (areas2D(poinsin(iv,1), poinsin(iv,2)) == zeroRK) THEN
            areas2D(poinsin(iv,1), poinsin(iv,2)) = areas(ix,iy,iv)
          END IF
        END DO

      END DO
    END DO

    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
    IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)

  END SUBROUTINE grid_spacepercen

  SUBROUTINE grid_spacepercen_providing_polys(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals,   &
    xBBvals, yBBvals, dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Nmaxvercoin,        &
    Nvercoinpolys, vercoinpolys, Ngridsin, gridsin,  areas, percentages)
  ! Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another 
  !   series of grid-cells (A) providing coincident polygons
  ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x 
  !   and y axes. 

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dxA, dyA, NAvertexmax
    INTEGER, INTENT(in)                                  :: dxB, dyB, NBvertexmax, dxyB
    INTEGER, INTENT(in)                                  :: Nmaxvercoin
    REAL(r_k), DIMENSION(dxA,dyA), INTENT(in)            :: xCAvals, yCAvals
    REAL(r_k), DIMENSION(dxB,dyB), INTENT(in)            :: xCBvals, yCBvals
    REAL(r_k), DIMENSION(dxA,dyA,NAvertexmax), INTENT(in):: xBAvals, yBAvals
    REAL(r_k), DIMENSION(dxB,dyB,NBvertexmax), INTENT(in):: xBBvals, yBBvals
    LOGICAL, INTENT(in)                                  :: strict
    INTEGER, DIMENSION(dxA,dyA,dxyB,Nmaxvercoin),                                                     &
      INTENT(out)                                        :: Nvercoinpolys
    REAL(r_k), DIMENSION(dxA,dyA,dxyB,Nmaxvercoin,2),                                                 &
      INTENT(out)                                        :: vercoinpolys
    INTEGER, DIMENSION(dxA,dyA), INTENT(out)             :: Ngridsin
    INTEGER, DIMENSION(dxA,dyA,dxyB,2), INTENT(out)      :: gridsin
    REAL(r_k), DIMENSION(dxB,dyB), INTENT(out)           :: areas
    REAL(r_k), DIMENSION(dxA,dyA,dxyB), INTENT(out)      :: percentages

! Local
   INTEGER                                               :: iv, ix, iy
   INTEGER                                               :: Nvertex
   INTEGER, ALLOCATABLE, DIMENSION(:,:)                  :: poinsin
   CHARACTER(len=20)                                     :: IS
   REAL(r_k), ALLOCATABLE, DIMENSION(:)                  :: pareas
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                :: vertexgrid

!!!!!!! Variables
! dxA, dyA: shape of the matrix with the grid points A
! xCAvals, yCAvals: coordinates of the center of the grid cells A
! NAvertexmax: Maximum number of vertices of the grid cells A
! xBAvals, yBAvals: coordinates of th vertices of the grid cells A (-99999 for no vertex)
! dxB, dyB: shape of the matrix with the grid points B
! xCBvals, yCBvals: coordinates of the center of the grid cells B
! NBvertexmax: Maximum number of vertices of the grid cells B
! xBBvals, yBBvals: coordinates of th vertices of the grid cells B (-99999 for no vertex)
! strict: give an error if the area of the polygon is not fully covered
! Nvercoinpolys: number of vertices of the coincident polygon of each grid
! coinpolys: of vertices of the coincident polygon of each grid
! Ngridsin: number of grids from grid B with some extension within the grid cell A
! gridsin: indices of B grids within the grids of A
! areas: areas of the grids
! percentages: percentages of area of cells B of each of the grids inside the grid cell A

    fname = 'grid_spacepercen_providing_polys'

    areas = zeroRK

    DO ix = 1, dxA
      DO iy = 1, dyA

        ! Getting grid vertices
        Nvertex = 0
        DO iv=1, NAvertexmax
          IF (xBAvals(ix,iy,iv) /= fillval64) THEN
           Nvertex = Nvertex + 1
          END IF
        END DO
        IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
        ALLOCATE(vertexgrid(Nvertex,2))
        vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex)
        vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex)
  
        CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xCBvals, yCBvals,            &
          NBvertexmax, xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,:,:))
   
        IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)
        ALLOCATE(poinsin(Ngridsin(ix,iy),2))
        IF (ALLOCATED(pareas)) DEALLOCATE(pareas)
        ALLOCATE(pareas(Ngridsin(ix,iy)))

        DO iv=1, Ngridsin(ix,iy)
          poinsin(iv,1) = gridsin(ix,iy,iv,1)
          poinsin(iv,2) = gridsin(ix,iy,iv,2)
        END DO

        CALL grid_spacepercen_within_reg_providing_polys(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, &
          xBBvals, yBBvals, Ngridsin(ix,iy), poinsin, strict, Nmaxvercoin, Nvercoinpolys(ix,iy,:,:), &
          vercoinpolys(ix,iy,:,:,:), pareas, percentages(ix,iy,:))

        ! Filling areas
        DO iv = 1, Ngridsin(ix,iy)
          IF (areas(poinsin(iv,1), poinsin(iv,2)) == zeroRK) THEN
            areas(poinsin(iv,1), poinsin(iv,2)) = pareas(iv)
          END IF
        END DO

      END DO
    END DO

    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
    IF (ALLOCATED(pareas)) DEALLOCATE(pareas)
    IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)

  END SUBROUTINE grid_spacepercen_providing_polys

  SUBROUTINE unique_matrixRK2D(dx, dy, dxy, matrix2D, Nunique, unique)
  ! Subroutine to provide the unique values within a 2D RK matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dx, dy, dxy
    REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: matrix2D
    INTEGER, INTENT(out)                                 :: Nunique
    REAL(r_k), DIMENSION(dxy), INTENT(out)               :: unique

! Local
    INTEGER                                              :: ix, iy, iu, minvalv
    LOGICAL                                              :: single
    REAL(r_k), ALLOCATABLE, DIMENSION(:)                 :: uniques


!!!!!!! Variables
! dx, dy: dimensions of the matrix
! dxy: dx*dy, maximum possible amount of different values
! matrix2D: matgrix of values
! Nunique: amount of unique values 
! unique: sorted from minimum to maximum vector with the unique values

    fname = 'unique_matrixRK2D'

    minvalv = MINVAL(matrix2D)

    Nunique = 1
    unique(1) = minvalv
    DO ix= 1, dx
      DO iy= 1, dy
        single = .TRUE.
        DO iu = 1, Nunique
          IF (matrix2D(ix,iy) == unique(iu)) THEN
            single = .FALSE.
            EXIT
          END IF
        END DO
        IF (single) THEN
          Nunique = Nunique + 1
          unique(Nunique) = matrix2D(ix,iy)
        END IF
      END DO
    END DO
    IF (ALLOCATED(uniques)) DEALLOCATE(uniques)
    ALLOCATE(uniques(Nunique)) 
    uniques(1:Nunique) = unique(1:Nunique)
    
    CALL sortR_K(uniques(1:Nunique), Nunique)
    unique(1:Nunique) = uniques(1:Nunique)

  END SUBROUTINE unique_matrixRK2D

  SUBROUTINE spaceweightstats(varin, Ngridsin, gridsin, percentages, stats, varout, dxA, dyA, dxB,    &
    dyB, maxNgridsin, Lstats)
  ! Subroutine to compute an spatial statistics value from a matrix B into a matrix A using weights

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dxA, dyA, dxB, dyB, maxNgridsin, Lstats
    CHARACTER(len=*), INTENT(in)                         :: stats
    INTEGER, DIMENSION(dxA,dyA), INTENT(in)              :: Ngridsin
    INTEGER, DIMENSION(dxA,dyA,maxNgridsin,2), INTENT(in):: gridsin
    REAL(r_k), DIMENSION(dxB,dyB), INTENT(in)            :: varin
    REAL(r_k), DIMENSION(dxA,dyA,maxNgridsin), INTENT(in):: percentages
    REAL(r_k), DIMENSION(dxA,dyA,Lstats), INTENT(out)    :: varout

! Local
    INTEGER                                              :: ix, iy, iv, ic, iu, ii, jj
    INTEGER                                              :: Ncounts
    CHARACTER(len=3)                                     :: val1S, val2S
    CHARACTER(len=30)                                    :: val3S
    REAL(r_k)                                            :: val1, val2
    REAL(r_k), DIMENSION(Lstats)                         :: icounts

!!!!!!! Variables
! dxA, dyA: length of dimensions of matrix A
! dxB, dyB: length of dimensions of matrix B
! maxNgridsin: maximum number of grid points from B to be used to compute into a grid of matrix A
! Lstats: length of the dimension of the statistics
! varin: variable from matrix B to be used
! Ngridsin: number of grids from matrix B for each grid of matrix A
! gridsin: coordinates of grids of matrix B for each grid of matrix A
! percentages: weights as percentages of space of grid in matrix A covered by grid of matrix B
! stats: name of the spatial statistics to compute inside each grid of matrix A using values from 
!     matrix B. Avaialbe ones:
!   'min': minimum value
!   'max': maximum value
!   'mean': space weighted mean value
!   'mean2': space weighted quadratic mean value
!   'stddev': space weighted standard deviation value
!   'count': percentage of the space of matrix A covered by each different value of matrix B
! varout: output statistical variable

    fname = 'spaceweightstats'

    ! Let's be efficvient?
    statn: SELECT CASE(TRIM(stats))
      CASE('min')
        varout = fillVal64
        DO ix=1, dxA
          DO iy=1, dyA
            DO iv=1, Ngridsin(ix,iy)
              ii = gridsin(ix,iy,iv,1)
              jj = gridsin(ix,iy,iv,2)
              IF (varin(ii,jj) < varout(ix,iy,Lstats)) varout(ix,iy,1) = varin(ii,jj)
            END DO
          END DO
        END DO
      CASE('max')
        varout = -fillVal64
        DO ix=1, dxA
          DO iy=1, dyA
            DO iv=1, Ngridsin(ix,iy)
              ii = gridsin(ix,iy,iv,1)
              jj = gridsin(ix,iy,iv,2)
              IF (varin(ii,jj) > varout(ix,iy,Lstats)) varout(ix,iy,1) = varin(ii,jj)
            END DO
          END DO
        END DO
      CASE('mean')
        varout = zeroRK
        DO ix=1, dxA
          DO iy=1, dyA
            DO iv=1, Ngridsin(ix,iy)
              ii = gridsin(ix,iy,iv,1)
              jj = gridsin(ix,iy,iv,2)
              varout(ix,iy,1) = varout(ix,iy,1) + varin(ii,jj)*percentages(ix,iy,iv)
            END DO
          END DO
        END DO
      CASE('mean2')
        varout = zeroRK
        DO ix=1, dxA
          DO iy=1, dyA
            DO iv=1, Ngridsin(ix,iy)
              ii = gridsin(ix,iy,iv,1)
              jj = gridsin(ix,iy,iv,2)
              varout(ix,iy,1) = varout(ix,iy,1) + percentages(ix,iy,iv)*(varin(ii,jj))**2
            END DO
            varout(ix,iy,1) = varout(ix,iy,1) / Ngridsin(ix,iy)
          END DO
        END DO
      CASE('stddev')
        varout = zeroRK
        DO ix=1, dxA
          DO iy=1, dyA
            val1 = zeroRK
            val2 = zeroRK
            DO iv=1, Ngridsin(ix,iy)
              ii = gridsin(ix,iy,iv,1)
              jj = gridsin(ix,iy,iv,2)
              val1 = val1 + varin(ii,jj)*percentages(ix,iy,iv)
              val2 = val2 + percentages(ix,iy,iv)*(varin(ii,jj))**2
            END DO
            varout(ix,iy,1) = SQRT(val2 - val1**2)
          END DO
        END DO
      CASE('count')
        CALL unique_matrixRK2D(dxB, dyB, dxB*dyB, varin, Ncounts, icounts)
        IF (Lstats /= Ncounts) THEN
          PRINT *,'  ' // TRIM(fname) // 'provided:', Lstats
          PRINT *,'  ' // TRIM(fname) // 'found:', Ncounts, ' :', icounts
          WRITE(val1S,'(I3)')Lstats
          WRITE(val2S,'(I3)')Ncounts
          msg = "for 'count' different amount of passed categories: " // TRIM(val1S) //               &
            ' and found ' // TRIM(val2S) 
          CALL ErrMsg(msg, fname, -1)
        END IF
        varout = zeroRK
        DO ix=1, dxA
          DO iy=1, dyA
            DO iv=1, Ngridsin(ix,iy)
              ii = gridsin(ix,iy,iv,1)
              jj = gridsin(ix,iy,iv,2)
              ic = Index1DArrayR_K(icounts, Ncounts, varin(ii,jj))
              IF (ic == -1) THEN
                WRITE(val3S,'(f30.20)')varin(ii,jj)
                msg = "value '" // val3S // "' for 'count' not found"
                CALL ErrMSg(msg, fname, -1)
              ELSE
                varout(ix,iy,ic) = varout(ix,iy,ic) + percentages(ix,iy,iv)
              END IF
            END DO
          END DO
        END DO
      CASE DEFAULT
        msg = "statisitcs '" // TRIM(stats) // "' not ready !!" // CHAR(44) // " available ones: " // &
          "'min', 'max', 'mean', 'mean2', 'stddev', 'count'" 
        CALL ErrMsg(msg, fname, -1)
    END SELECT statn

  END SUBROUTINE spaceweightstats

  SUBROUTINE multi_spaceweightstats_in1DRKno_slc3v3(varin, idv, Ngridsin, gridsin, percentages,       &
    varout, di1, ds1, ds2, ds3, maxNgridsin)
  ! Subroutine to compute an spatial statistics value from a 1D RK matrix without running one into a 
  !   matrix of 3-variables slices of rank 3 using spatial weights

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: di1, idv, ds1, ds2, ds3
    INTEGER, INTENT(in)                                  :: maxNgridsin
    INTEGER, DIMENSION(ds1,ds2,ds3), INTENT(in)          :: Ngridsin
    INTEGER, DIMENSION(ds1,ds2,ds3,maxNgridsin,2),                                                    &
      INTENT(in)                                         :: gridsin
    REAL(r_k), DIMENSION(di1), INTENT(in)                :: varin
    REAL(r_k), INTENT(in),                                                                            &
      DIMENSION(ds1,ds2,ds3,maxNgridsin)                 :: percentages
    REAL(r_k), DIMENSION(ds1,ds2,ds3,7), INTENT(out)     :: varout

! Local
    INTEGER                                              :: i1, i2, i3, s1, s2, s3, iv
    INTEGER                                              :: ii3, ss1, ss2, ss3
    INTEGER                                              :: Ncounts, Nin
    INTEGER, DIMENSION(1)                                :: dmaxvarin
    CHARACTER(len=3)                                     :: val1S, val2S
    CHARACTER(len=30)                                    :: val3S
    REAL(r_k)                                            :: minv, maxv, meanv, mean2v, stdv, medv
    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: pin
    INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: gin
    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: svin
    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: vin

!!!!!!! Variables
! di1: length of dimensions of the 1D matrix of values
! ds[1-3]: length of dimensions of matrix with the slices
! maxNgridsin: maximum number of grid points from the 3D matrix in any slice
! varin: 1D RK variable to be used
! idv: which dimension of the sliced grids coincide with the dimension of 1D varin
! Ngridsin: number of grids from 3D RK matrix for each slice
! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices
! percentages: weights as percentages of space of 3D RK matrix for each slice
!!!!!
! Available spatial statistics to compute inside each slice using values from 3D RK matrix
!   'min': minimum value
!   'max': maximum value
!   'mean': space weighted mean value
!   'mean2': space weighted quadratic mean value
!   'stddev': space weighted standard deviation value
!   'median': median value
!   'count': percentage of the space of matrix A covered by each different value of matrix B
! varout: output statistical variable

    fname = 'multi_spaceweightstats_in1DRKno_slc3v3'

    varout = fillval64

    ss1 = 8 + 1
    ss2 = 5 + 1
    ss3 = 3 + 1
    ii3 = 1 + 1

    dmaxvarin = UBOUND(varin)

    ! Let's be efficient?
    varout = fillVal64
    DO s1 =1, ds1
      DO s2 =1, ds2
        DO s3 =1, ds3
          Nin = Ngridsin(s1,s2,s3)
          ! Computing along d3
          IF (Nin > 1) THEN
            IF (ALLOCATED(gin)) DEALLOCATE(gin)
            ALLOCATE(gin(Nin,2))
            IF (ALLOCATED(pin)) DEALLOCATE(pin)
            ALLOCATE(pin(Nin))
            IF (ALLOCATED(vin)) DEALLOCATE(vin)
            ALLOCATE(vin(Nin))
            IF (ALLOCATED(svin)) DEALLOCATE(svin)
            ALLOCATE(svin(Nin))
            gin = gridsin(s1,s2,s3,1:Nin,:)
            pin = percentages(s1,s2,s3,1:Nin)

            ! Getting the values
            DO iv=1, Nin
              i1 = gin(iv,idv)
              vin(iv) = varin(i1)
            END DO
            minv = fillVal64
            maxv = -fillVal64
            meanv = zeroRK
            mean2v = zeroRK
            stdv = zeroRK
            minv = MINVAL(vin)
            maxv = MAXVAL(vin)
            meanv = SUM(vin*pin)
            mean2v = SUM(vin**2*pin)
            DO iv=1,Nin
              stdv = stdv + ( (meanv - vin(iv))*pin(iv) )**2
            END DO
            stdv = SQRT(stdv)
            svin = vin(:)
            CALL SortR_K(svin, Nin)
            medv = svin(INT(Nin/2))
            varout(s1,s2,s3,1) = minv
            varout(s1,s2,s3,2) = maxv
            varout(s1,s2,s3,3) = meanv
            varout(s1,s2,s3,4) = mean2v
            varout(s1,s2,s3,5) = stdv
            varout(s1,s2,s3,6) = medv
            varout(s1,s2,s3,7) = Nin*1.
          ELSE
            i1 = gridsin(s1,s2,s3,1,idv)
            IF (i1 > 0 .AND. i1 <= dmaxvarin(1)) THEN
              varout(s1,s2,s3,1) = varin(i1)
              varout(s1,s2,s3,2) = varin(i1)
              varout(s1,s2,s3,3) = varin(i1)
              varout(s1,s2,s3,4) = varin(i1)*varin(i1)
              varout(s1,s2,s3,5) = zeroRK
              varout(s1,s2,s3,6) = varin(i1)
              varout(s1,s2,s3,7) = Nin*1.
            END IF
          END IF
        END DO
      END DO
    END DO

    IF (ALLOCATED(gin)) DEALLOCATE(gin)
    IF (ALLOCATED(pin)) DEALLOCATE(pin)
    IF (ALLOCATED(vin)) DEALLOCATE(vin)
    IF (ALLOCATED(svin)) DEALLOCATE(svin)
    
    RETURN

  END SUBROUTINE multi_spaceweightstats_in1DRKno_slc3v3

  SUBROUTINE multi_spaceweightstats_in2DRKno_slc3v3(varin, Ngridsin, gridsin, percentages, varout,     &
    di1, di2, ds1, ds2, ds3, maxNgridsin)
  ! Subroutine to compute an spatial statistics value from a 2D RK matrix without running one into a 
  !   matrix of 3-variables slices of rank 3 using spatial weights

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: di1, di2, ds1, ds2, ds3
    INTEGER, INTENT(in)                                  :: maxNgridsin
    INTEGER, DIMENSION(ds1,ds2,ds3), INTENT(in)          :: Ngridsin
    INTEGER, DIMENSION(ds1,ds2,ds3,maxNgridsin,2),                                                    &
      INTENT(in)                                         :: gridsin
    REAL(r_k), DIMENSION(di1,di2), INTENT(in)            :: varin
    REAL(r_k), INTENT(in),                                                                            &
      DIMENSION(ds1,ds2,ds3,maxNgridsin)                 :: percentages
    REAL(r_k), DIMENSION(ds1,ds2,ds3,7), INTENT(out)     :: varout

! Local
    INTEGER                                              :: i1, i2, i3, s1, s2, s3, iv
    INTEGER                                              :: ii3, ss1, ss2, ss3
    INTEGER                                              :: Ncounts, Nin
    CHARACTER(len=3)                                     :: val1S, val2S
    CHARACTER(len=30)                                    :: val3S
    REAL(r_k)                                            :: minv, maxv, meanv, mean2v, stdv, medv
    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: pin
    INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: gin
    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: svin
    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: vin

!!!!!!! Variables
! di1, di2: length of dimensions of the 2D matrix of values
! ds[1-3]: length of dimensions of matrix with the slices
! maxNgridsin: maximum number of grid points from the 3D matrix in any slice
! varin: 2D RK variable to be used
! Ngridsin: number of grids from 3D RK matrix for each slice
! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices
! percentages: weights as percentages of space of 3D RK matrix for each slice
!!!!!
! Available spatial statistics to compute inside each slice using values from 3D RK matrix
!   'min': minimum value
!   'max': maximum value
!   'mean': space weighted mean value
!   'mean2': space weighted quadratic mean value
!   'stddev': space weighted standard deviation value
!   'median': median value
!   'count': percentage of the space of matrix A covered by each different value of matrix B
! varout: output statistical variable

    fname = 'multi_spaceweightstats_in2DRKno_slc3v3'

    varout = fillval64

    ss1 = 8 + 1
    ss2 = 5 + 1
    ss3 = 3 + 1
    ii3 = 1 + 1

    ! Let's be efficient?
    varout = fillVal64
    DO s1 =1, ds1
      DO s2 =1, ds2
        DO s3 =1, ds3
          Nin = Ngridsin(s1,s2,s3)
          ! Computing along d3
          IF (Nin > 1) THEN
            IF (ALLOCATED(gin)) DEALLOCATE(gin)
            ALLOCATE(gin(Nin,2))
            IF (ALLOCATED(pin)) DEALLOCATE(pin)
            ALLOCATE(pin(Nin))
            IF (ALLOCATED(vin)) DEALLOCATE(vin)
            ALLOCATE(vin(Nin))
            IF (ALLOCATED(svin)) DEALLOCATE(svin)
            ALLOCATE(svin(Nin))
            gin = gridsin(s1,s2,s3,1:Nin,:)
            pin = percentages(s1,s2,s3,1:Nin)

            ! Getting the values
            DO iv=1, Nin
              i1 = gin(iv,1)
              i2 = gin(iv,2)
              vin(iv) = varin(i1,i2)
            END DO
            minv = fillVal64
            maxv = -fillVal64
            meanv = zeroRK
            mean2v = zeroRK
            stdv = zeroRK
            minv = MINVAL(vin)
            maxv = MAXVAL(vin)
            meanv = SUM(vin*pin)
            mean2v = SUM(vin**2*pin)
            DO iv=1,Nin
              stdv = stdv + ( (meanv - vin(iv))*pin(iv) )**2
            END DO
            stdv = SQRT(stdv)
            svin = vin(:)
            CALL SortR_K(svin, Nin)
            medv = svin(INT(Nin/2))
            varout(s1,s2,s3,1) = minv
            varout(s1,s2,s3,2) = maxv
            varout(s1,s2,s3,3) = meanv
            varout(s1,s2,s3,4) = mean2v
            varout(s1,s2,s3,5) = stdv
            varout(s1,s2,s3,6) = medv
            varout(s1,s2,s3,7) = Nin*1.
          ELSE
            i1 = gridsin(s1,s2,s3,1,1)
            i2 = gridsin(s1,s2,s3,1,2)
            varout(s1,s2,s3,1) = varin(i1,i2)
            varout(s1,s2,s3,2) = varin(i1,i2)
            varout(s1,s2,s3,3) = varin(i1,i2)
            varout(s1,s2,s3,4) = varin(i1,i2)*varin(i1,i2)
            varout(s1,s2,s3,5) = zeroRK
            varout(s1,s2,s3,6) = varin(i1,i2)
            varout(s1,s2,s3,7) = Nin*1.
          END IF
        END DO
      END DO
    END DO

    IF (ALLOCATED(gin)) DEALLOCATE(gin)
    IF (ALLOCATED(pin)) DEALLOCATE(pin)
    IF (ALLOCATED(vin)) DEALLOCATE(vin)
    IF (ALLOCATED(svin)) DEALLOCATE(svin)
    
    RETURN

  END SUBROUTINE multi_spaceweightstats_in2DRKno_slc3v3

  SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v3(varin, Ngridsin, gridsin, percentages, varout,     &
    di1, di2, di3, ds1, ds2, ds3, maxNgridsin)
  ! Subroutine to compute an spatial statistics value from a 3D RK matrix using 3rd dimension as 
  !   running one into a matrix of 3-variables slices of rank 3 using spatial weights

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: di1, di2, di3, ds1, ds2, ds3
    INTEGER, INTENT(in)                                  :: maxNgridsin
    INTEGER, DIMENSION(ds1,ds2,ds3), INTENT(in)          :: Ngridsin
    INTEGER, DIMENSION(ds1,ds2,ds3,maxNgridsin,2),                                                    &
      INTENT(in)                                         :: gridsin
    REAL(r_k), DIMENSION(di1,di2,di3), INTENT(in)        :: varin
    REAL(r_k), INTENT(in),                                                                            &
      DIMENSION(ds1,ds2,ds3,maxNgridsin)                 :: percentages
    REAL(r_k), DIMENSION(ds1,ds2,ds3,di3,7), INTENT(out) :: varout

! Local
    INTEGER                                              :: i1, i2, i3, s1, s2, s3, iv
    INTEGER                                              :: ii3, ss1, ss2, ss3
    INTEGER                                              :: Ncounts, Nin
    CHARACTER(len=3)                                     :: val1S, val2S
    CHARACTER(len=30)                                    :: val3S
    REAL(r_k)                                            :: minv, maxv, meanv, mean2v, stdv, medv
    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: pin
    INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: gin
    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: svin
    REAL(r_k), DIMENSION(:,:), ALLOCATABLE               :: vin

!!!!!!! Variables
! di1, di2, di3: length of dimensions of the 3D matrix of values
! ds[1-3]: length of dimensions of matrix with the slices
! maxNgridsin: maximum number of grid points from the 3D matrix in any slice
! varin: 3D RK variable to be used
! Ngridsin: number of grids from 3D RK matrix for each slice
! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices
! percentages: weights as percentages of space of 3D RK matrix for each slice
!!!!!
! Available spatial statistics to compute inside each slice using values from 3D RK matrix
!   'min': minimum value
!   'max': maximum value
!   'mean': space weighted mean value
!   'mean2': space weighted quadratic mean value
!   'stddev': space weighted standard deviation value
!   'median': median value
!   'count': percentage of the space of matrix A covered by each different value of matrix B
! varout: output statistical variable

    fname = 'multi_spaceweightstats_in3DRK3_slc3v3'

    varout = fillval64

    ss1 = 8 + 1
    ss2 = 5 + 1
    ss3 = 3 + 1
    ii3 = 1 + 1

    ! Let's be efficient?
    varout = fillVal64
    DO s1 =1, ds1
      DO s2 =1, ds2
        DO s3 =1, ds3
          Nin = Ngridsin(s1,s2,s3)
          ! Computing along d3
          IF (Nin > 1) THEN
            IF (ALLOCATED(gin)) DEALLOCATE(gin)
            ALLOCATE(gin(Nin,2))
            IF (ALLOCATED(pin)) DEALLOCATE(pin)
            ALLOCATE(pin(Nin))
            IF (ALLOCATED(vin)) DEALLOCATE(vin)
            ALLOCATE(vin(Nin,di3))
            IF (ALLOCATED(svin)) DEALLOCATE(svin)
            ALLOCATE(svin(Nin))
            gin = gridsin(s1,s2,s3,1:Nin,:)
            pin = percentages(s1,s2,s3,1:Nin)

            ! Getting the values
            DO iv=1, Nin
              i1 = gin(iv,1)
              i2 = gin(iv,2)
              vin(iv,:) = varin(i1,i2,:)
            END DO
            DO i3=1, di3
              minv = fillVal64
              maxv = -fillVal64
              meanv = zeroRK
              mean2v = zeroRK
              stdv = zeroRK
              minv = MINVAL(vin(:,i3))
              maxv = MAXVAL(vin(:,i3))
              meanv = SUM(vin(:,i3)*pin)
              mean2v = SUM(vin(:,i3)**2*pin)
              DO iv=1,Nin
                stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2
              END DO
              stdv = SQRT(stdv)
              svin = vin(:,i3)
              CALL SortR_K(svin, Nin)
              medv = svin(INT(Nin/2))
              varout(s1,s2,s3,i3,1) = minv
              varout(s1,s2,s3,i3,2) = maxv
              varout(s1,s2,s3,i3,3) = meanv
              varout(s1,s2,s3,i3,4) = mean2v
              varout(s1,s2,s3,i3,5) = stdv
              varout(s1,s2,s3,i3,6) = medv
              varout(s1,s2,s3,i3,7) = Nin*1.
            END DO
          ELSE
            i1 = gridsin(s1,s2,s3,1,1)
            i2 = gridsin(s1,s2,s3,1,2)
            varout(s1,s2,s3,:,1) = varin(i1,i2,:)
            varout(s1,s2,s3,:,2) = varin(i1,i2,:)
            varout(s1,s2,s3,:,3) = varin(i1,i2,:)
            varout(s1,s2,s3,:,4) = varin(i1,i2,:)*varin(i1,i2,:)
            varout(s1,s2,s3,:,5) = zeroRK
            varout(s1,s2,s3,:,6) = varin(i1,i2,:)
            varout(s1,s2,s3,:,7) = Nin*1.
          END IF
        END DO
      END DO
    END DO

    IF (ALLOCATED(gin)) DEALLOCATE(gin)
    IF (ALLOCATED(pin)) DEALLOCATE(pin)
    IF (ALLOCATED(vin)) DEALLOCATE(vin)
    IF (ALLOCATED(svin)) DEALLOCATE(svin)
    
    RETURN

  END SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v3

  SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v4(varin, Ngridsin, gridsin, percentages, varout,     &
    di1, di2, di3, ds1, ds2, ds3, ds4, maxNgridsin)
  ! Subroutine to compute an spatial statistics value from a 3D RK matrix using 3rd dimension as 
  !   running one into a matrix of 3-variables slices of rank 4 using spatial weights

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: di1, di2, di3, ds1, ds2, ds3, ds4
    INTEGER, INTENT(in)                                  :: maxNgridsin
    INTEGER, DIMENSION(ds1,ds2,ds3,ds4), INTENT(in)      :: Ngridsin
    INTEGER, INTENT(in),                                                                              &
      DIMENSION(ds1,ds2,ds3,ds4,maxNgridsin,2)           :: gridsin
    REAL(r_k), DIMENSION(di1,di2,di3), INTENT(in)        :: varin
    REAL(r_k), INTENT(in),                                                                            &
      DIMENSION(ds1,ds2,ds3,ds4,maxNgridsin)             :: percentages
    REAL(r_k), DIMENSION(ds1,ds2,ds3,ds4,di3,7),                                                      &
      INTENT(out)                                        :: varout

! Local
    INTEGER                                              :: i1, i2, i3, s1, s2, s3, s4, iv
    INTEGER                                              :: Ncounts, Nin
    CHARACTER(len=3)                                     :: val1S, val2S
    CHARACTER(len=30)                                    :: val3S
    REAL(r_k)                                            :: minv, maxv, meanv, mean2v, stdv, medv
    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: pin
    INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: gin
    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: svin
    REAL(r_k), DIMENSION(:,:), ALLOCATABLE               :: vin

!!!!!!! Variables
! di1, di2, di3: length of dimensions of the 3D matrix of values
! ds[1-4]: length of dimensions of matrix with the slices
! maxNgridsin: maximum number of grid points from the 3D matrix in any slice
! varin: 3D RK variable to be used
! Ngridsin: number of grids from 3D RK matrix for each slice
! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices
! percentages: weights as percentages of space of 3D RK matrix for each slice
!!!!!
! Available spatial statistics to compute inside each slice using values from 3D RK matrix
!   'min': minimum value
!   'max': maximum value
!   'mean': space weighted mean value
!   'mean2': space weighted quadratic mean value
!   'stddev': space weighted standard deviation value
!   'median': median value
!   'count': percentage of the space of matrix A covered by each different value of matrix B
! varout: output statistical variable

    fname = 'multi_spaceweightstats_in3DRK3_slc3v4'

    varout = fillval64

    ! Let's be efficient?
    varout = fillVal64
    DO s1 =1, ds1
      DO s2 =1, ds2
        DO s3 =1, ds3
          DO s4 =1, ds4
            Nin = Ngridsin(s1,s2,s3,s4)
            IF (Nin > 1) THEN
              IF (ALLOCATED(gin)) DEALLOCATE(gin)
              ALLOCATE(gin(Nin,2))
              IF (ALLOCATED(pin)) DEALLOCATE(pin)
              ALLOCATE(pin(Nin))
              IF (ALLOCATED(vin)) DEALLOCATE(vin)
              ALLOCATE(vin(Nin,di3))
              IF (ALLOCATED(svin)) DEALLOCATE(svin)
              ALLOCATE(svin(Nin))
              gin = gridsin(s1,s2,s3,s4,1:Nin,:)
              pin = percentages(s1,s2,s3,s4,1:Nin)

              ! Getting the values
              DO iv=1, Nin
                i1 = gin(iv,1)
                i2 = gin(iv,2)
                vin(iv,:) = varin(i1,i2,:)
              END DO
              ! Computing along d3
              DO i3=1, di3
                minv = fillVal64
                maxv = -fillVal64
                meanv = zeroRK
                mean2v = zeroRK
                stdv = zeroRK

                minv = MINVAL(vin(:,i3))
                maxv = MAXVAL(vin(:,i3))
                meanv = SUM(vin(:,i3)*pin)
                mean2v = SUM(vin(:,i3)**2*pin)  
                DO iv=1,Nin
                  stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2
                END DO
                stdv = SQRT(stdv)
                svin = vin(:,i3)
                CALL SortR_K(svin, Nin)
                medv = svin(INT(Nin/2))
                varout(s1,s2,s3,s4,i3,1) = minv
                varout(s1,s2,s3,s4,i3,2) = maxv
                varout(s1,s2,s3,s4,i3,3) = meanv
                varout(s1,s2,s3,s4,i3,4) = mean2v
                varout(s1,s2,s3,s4,i3,5) = stdv
                varout(s1,s2,s3,s4,i3,6) = medv
                varout(s1,s2,s3,s4,i3,7) = Nin*1.
              END DO
            ELSE
                i1 = gridsin(s1,s2,s3,s4,1,1)
                i2 = gridsin(s1,s2,s3,s4,1,2)
                varout(s1,s2,s3,s4,:,1) = varin(i1,i2,:)
                varout(s1,s2,s3,s4,:,2) = varin(i1,i2,:)
                varout(s1,s2,s3,s4,:,3) = varin(i1,i2,:)
                varout(s1,s2,s3,s4,:,4) = varin(i1,i2,:)*varin(i1,i2,:)
                varout(s1,s2,s3,s4,:,5) = zeroRK
                varout(s1,s2,s3,s4,:,6) = varin(i1,i2,:)
                varout(s1,s2,s3,s4,:,7) = Nin*1.
            END IF
          END DO
        END DO
      END DO
    END DO

    IF (ALLOCATED(gin)) DEALLOCATE(gin)
    IF (ALLOCATED(pin)) DEALLOCATE(pin)
    IF (ALLOCATED(vin)) DEALLOCATE(vin)
    IF (ALLOCATED(svin)) DEALLOCATE(svin)
    
    RETURN

  END SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v4

  SUBROUTINE multi_index_mat2DI(d1, d2, d12, mat, value, Nindices, indices)
  ! Subroutine to provide the indices of the different locations of a value inside a 2D integer matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d12
    INTEGER, DIMENSION(d1,d2), INTENT(in)                :: mat
    INTEGER,INTENT(in)                                   :: value
    INTEGER, INTENT(out)                                 :: Nindices
    INTEGER, DIMENSION(2,d12), INTENT(out)               :: indices

! Local
    INTEGER                                              :: i,j
    INTEGER                                              :: Ncounts1D, icount1D
    INTEGER, DIMENSION(d2)                               :: diffmat1D

    !!!!!!! Variables
    ! d1, d2: shape of the 2D matrix
    ! mat: 2D matrix
    ! value: value to be looking for
    ! Nindices: number of times value found within matrix
    ! indices: indices of the found values

    fname = 'multi_index_mat2DI'

    Nindices = 0
    indices = 0
    DO i=1, d1
      diffmat1D = mat(i,:) - value
      IF (ANY(diffmat1D == 0)) THEN
        Ncounts1D = COUNT(diffmat1D == 0)
        icount1D = 0
        DO j=1, d2
          IF (diffmat1D(j) == 0) THEN
            Nindices = Nindices + 1
            indices(1,Nindices) = i
            indices(2,Nindices) = j
            icount1D = icount1D + 1
            IF (icount1D == Ncounts1D) EXIT
          END IF
        END DO
      END IF
    END DO

  END SUBROUTINE multi_index_mat2DI

  SUBROUTINE multi_index_mat3DI(d1, d2, d3, d123, mat, value, Nindices, indices)
  ! Subroutine to provide the indices of the different locations of a value inside a 3D integer matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3, d123
    INTEGER, DIMENSION(d1,d2,d3), INTENT(in)             :: mat
    INTEGER, INTENT(in)                                  :: value
    INTEGER, INTENT(out)                                 :: Nindices
    INTEGER, DIMENSION(3,d123), INTENT(out)              :: indices

! Local
    INTEGER                                              :: i,j,k
    INTEGER                                              :: Ncounts1D, icount1D
    INTEGER                                              :: Ncounts2D, icount2D
    INTEGER, DIMENSION(d2,d3)                            :: diffmat2D
    INTEGER, DIMENSION(d3)                               :: diffmat1D

    !!!!!!! Variables
    ! d1, d2, d3: shape of the 3D matrix
    ! mat: 3D matrix
    ! value: value to be looking for
    ! Nindices: number of times value found within matrix
    ! indices: indices of the found values

    fname = 'multi_index_mat3DI'

    Nindices = 0
    indices = 0
    DO i=1, d1
      diffmat2D = mat(i,:,:) - value
      IF (ANY(diffmat2D == 0)) THEN
        Ncounts2D = COUNT(diffmat2D == 0)
        icount2D = 0
        DO j=1, d2
          diffmat1D = mat(i,j,:) - value
          IF (ANY(diffmat1D == 0)) THEN
            Ncounts1D = COUNT(diffmat1D == 0)
            icount1D = 0
            DO k=1, d3
              IF (diffmat1D(k) == 0) THEN
                Nindices = Nindices + 1
                indices(1,Nindices) = i
                indices(2,Nindices) = j
                indices(3,Nindices) = k
                icount1D = icount1D + 1
                IF (icount1D == Ncounts1D) EXIT
              END IF
            END DO
            icount2D = icount2D + icount1D
            IF (icount2D == Ncounts2D) EXIT
          END IF
        END DO
      END IF
    END DO

  END SUBROUTINE multi_index_mat3DI

  SUBROUTINE multi_index_mat4DI(d1, d2, d3, d4, d1234, mat, value, Nindices, indices)
  ! Subroutine to provide the indices of the different locations of a value inside a 4D integer matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4, d1234
    INTEGER, DIMENSION(d1,d2,d3,d4), INTENT(in)          :: mat
    INTEGER, INTENT(in)                                  :: value
    INTEGER, INTENT(out)                                 :: Nindices
    INTEGER, DIMENSION(4,d1234), INTENT(out)             :: indices

! Local
    INTEGER                                              :: i,j,k,l
    INTEGER                                              :: Ncounts1D, icount1D
    INTEGER                                              :: Ncounts2D, icount2D
    INTEGER                                              :: Ncounts3D, icount3D
    INTEGER, DIMENSION(d2,d3,d4)                         :: diffmat3D
    INTEGER, DIMENSION(d3,d4)                            :: diffmat2D
    INTEGER, DIMENSION(d4)                               :: diffmat1D

    !!!!!!! Variables
    ! d1, d2, d3, d4: shape of the 4D matrix
    ! mat: 4D matrix
    ! value: value to be looking for
    ! Nindices: number of times value found within matrix
    ! indices: indices of the found values

    fname = 'multi_index_mat4DI'

    Nindices = 0
    indices = 0
    DO i=1, d1
      diffmat3D = mat(i,:,:,:) - value
      IF (ANY(diffmat3D == 0)) THEN
        Ncounts3D = COUNT(diffmat3D == 0)
        icount3D = 0
        DO j=1, d2
          diffmat2D = mat(i,j,:,:) - value
          IF (ANY(diffmat2D == 0)) THEN
            Ncounts2D = COUNT(diffmat2D == 0)
            icount2D = 0
            DO k=1, d3
              diffmat1D = mat(i,j,k,:) - value
              IF (ANY(diffmat1D == 0)) THEN
                Ncounts1D = COUNT(diffmat1D == 0)
                icount1D = 0
                DO l=1, d4
                  IF (diffmat1D(l) == 0) THEN
                    Nindices = Nindices + 1
                    indices(1,Nindices) = i
                    indices(2,Nindices) = j
                    indices(3,Nindices) = k
                    indices(4,Nindices) = l
                    icount1D = icount1D + 1
                    IF (icount1D == Ncounts1D) EXIT
                  END IF
                END DO
              icount2D = icount2D + icount1D
              IF (icount2D == Ncounts2D) EXIT
              END IF
            END DO
            icount3D = icount3D + icount1D
            IF (icount3D == Ncounts3D) EXIT
          END IF
        END DO
      END IF
    END DO

  END SUBROUTINE multi_index_mat4DI

  SUBROUTINE multi_index_mat2DRK(d1, d2, d12, mat, value, Nindices, indices)
  ! Subroutine to provide the indices of the different locations of a value inside a 2D RK matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d12
    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: mat
    REAL(r_k),INTENT(in)                                 :: value
    INTEGER, INTENT(out)                                 :: Nindices
    INTEGER, DIMENSION(2,d12), INTENT(out)               :: indices

! Local
    INTEGER                                              :: i,j
    INTEGER                                              :: Ncounts1D, icount1D
    REAL(r_k), DIMENSION(d2)                             :: diffmat1D

    !!!!!!! Variables
    ! d1, d2: shape of the 2D matrix
    ! mat: 2D matrix
    ! value: value to be looking for
    ! Nindices: number of times value found within matrix
    ! indices: indices of the found values

    fname = 'multi_index_mat2DRK'

    Nindices = 0
    indices = 0
    DO i=1, d1
      diffmat1D = mat(i,:) - value
      IF (ANY(diffmat1D == zeroRK)) THEN
        Ncounts1D = COUNT(diffmat1D == zeroRK)
        icount1D = 0
        DO j=1, d2
          IF (diffmat1D(j) == zeroRK) THEN
            Nindices = Nindices + 1
            indices(1,Nindices) = i
            indices(2,Nindices) = j
            icount1D = icount1D + 1
            IF (icount1D == Ncounts1D) EXIT
          END IF
        END DO
      END IF
    END DO

  END SUBROUTINE multi_index_mat2DRK

  SUBROUTINE multi_index_mat3DRK(d1, d2, d3, d123, mat, value, Nindices, indices)
  ! Subroutine to provide the indices of the different locations of a value inside a 3D RK matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3, d123
    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: mat
    REAL(r_k),INTENT(in)                                 :: value
    INTEGER, INTENT(out)                                 :: Nindices
    INTEGER, DIMENSION(3,d123), INTENT(out)              :: indices

! Local
    INTEGER                                              :: i,j,k
    INTEGER                                              :: Ncounts1D, icount1D
    INTEGER                                              :: Ncounts2D, icount2D
    REAL(r_k), DIMENSION(d2,d3)                          :: diffmat2D
    REAL(r_k), DIMENSION(d3)                             :: diffmat1D

    !!!!!!! Variables
    ! d1, d2, d3: shape of the 3D matrix
    ! mat: 3D matrix
    ! value: value to be looking for
    ! Nindices: number of times value found within matrix
    ! indices: indices of the found values

    fname = 'multi_index_mat3DRK'

    Nindices = 0
    indices = 0
    DO i=1, d1
      diffmat2D = mat(i,:,:) - value
      IF (ANY(diffmat2D == zeroRK)) THEN
        Ncounts2D = COUNT(diffmat2D == zeroRK)
        icount2D = 0
        DO j=1, d2
          diffmat1D = mat(i,j,:) - value
          IF (ANY(diffmat1D == zeroRK)) THEN
            Ncounts1D = COUNT(diffmat1D == zeroRK)
            icount1D = 0
            DO k=1, d3
              IF (diffmat1D(k) == zeroRK) THEN
                Nindices = Nindices + 1
                indices(1,Nindices) = i
                indices(2,Nindices) = j
                indices(3,Nindices) = k
                icount1D = icount1D + 1
                IF (icount1D == Ncounts1D) EXIT
              END IF
            END DO
            icount2D = icount2D + icount1D
            IF (icount2D == Ncounts2D) EXIT
          END IF
        END DO
      END IF
    END DO

  END SUBROUTINE multi_index_mat3DRK

  SUBROUTINE multi_index_mat4DRK(d1, d2, d3, d4, d1234, mat, value, Nindices, indices)
  ! Subroutine to provide the indices of the different locations of a value inside a 4D RK matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4, d1234
    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: mat
    REAL(r_k),INTENT(in)                                 :: value
    INTEGER, INTENT(out)                                 :: Nindices
    INTEGER, DIMENSION(4,d1234), INTENT(out)             :: indices

! Local
    INTEGER                                              :: i,j,k,l
    INTEGER                                              :: Ncounts1D, icount1D
    INTEGER                                              :: Ncounts2D, icount2D
    INTEGER                                              :: Ncounts3D, icount3D
    REAL(r_k), DIMENSION(d2,d3,d4)                       :: diffmat3D
    REAL(r_k), DIMENSION(d3,d4)                          :: diffmat2D
    REAL(r_k), DIMENSION(d4)                             :: diffmat1D

    !!!!!!! Variables
    ! d1, d2, d3, d4: shape of the 4D matrix
    ! mat: 4D matrix
    ! value: value to be looking for
    ! Nindices: number of times value found within matrix
    ! indices: indices of the found values

    fname = 'multi_index_mat4DRK'

    Nindices = 0
    indices = 0
    DO i=1, d1
      diffmat3D = mat(i,:,:,:) - value
      IF (ANY(diffmat3D == zeroRK)) THEN
        Ncounts3D = COUNT(diffmat3D == zeroRK)
        icount3D = 0
        DO j=1, d2
          diffmat2D = mat(i,j,:,:) - value
          IF (ANY(diffmat2D == zeroRK)) THEN
            Ncounts2D = COUNT(diffmat2D == zeroRK)
            icount2D = 0
            DO k=1, d3
              diffmat1D = mat(i,j,k,:) - value
              IF (ANY(diffmat1D == zeroRK)) THEN
                Ncounts1D = COUNT(diffmat1D == zeroRK)
                icount1D = 0
                DO l=1, d4
                  IF (diffmat1D(l) == zeroRK) THEN
                    Nindices = Nindices + 1
                    indices(1,Nindices) = i
                    indices(2,Nindices) = j
                    indices(3,Nindices) = k
                    indices(4,Nindices) = l
                    icount1D = icount1D + 1
                    IF (icount1D == Ncounts1D) EXIT
                  END IF
                END DO
              icount2D = icount2D + icount1D
              IF (icount2D == Ncounts2D) EXIT
              END IF
            END DO
            icount3D = icount3D + icount1D
            IF (icount3D == Ncounts3D) EXIT
          END IF
        END DO
      END IF
    END DO

  END SUBROUTINE multi_index_mat4DRK

  SUBROUTINE coincident_list_2Dcoords(NpointsA, pointsA, NpointsB, pointsB, Npoints, points, inpA,    &
    inpB)
  ! Subroutine to determine which 2D points of an A list are also found in a B list

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: NpointsA, NpointsB
    INTEGER, DIMENSION(NpointsA,2), INTENT(in)           :: pointsA
    INTEGER, DIMENSION(NpointsB,2), INTENT(in)           :: pointsB
    INTEGER, INTENT(out)                                 :: Npoints
    INTEGER, DIMENSION(NpointsA,2), INTENT(out)          :: points
    INTEGER, DIMENSION(NpointsA), INTENT(out)            :: inpA, inpB

    ! Local
    INTEGER                                              :: iA, iB

!!!!!!! Variables
! NpointsA: Number of points of the list A
! pointsA: points of the list A
! NpointsB: Number of points of the list B
! pointsB: points of the list B
! Npoints: Number of coincident points
! points: coincident points
! inpA: coincident points list A
! inpB: coincident points list B


    fname = 'coincident_list_2Dcoords'

    Npoints = 0
    points = 0
    inpA = 0
    inpB = 0

    DO iA = 1, NpointsA
      DO iB = 1, NpointsB
        IF ( (pointsA(iA,1) == pointsB(iB,1)) .AND. (pointsA(iA,2) == pointsB(iB,2)) ) THEN
          Npoints = Npoints + 1
          points(Npoints,1) = pointsA(iA,1)
          points(Npoints,2) = pointsA(iA,2)
          inpA(Npoints) = iA
          inpB(Npoints) = iB
          EXIT
        END IF

      END DO
    END DO

  END SUBROUTINE coincident_list_2Dcoords

  SUBROUTINE coincident_gridsin2D_old(dxA, dyA, dxyA, NpointsA, pointsA, dxB, dyB, dxyB, NpointsB,    &
    pointsB, Npoints, points, inpointsA, inpointsB)
  ! Subroutine to determine which lists of 2D gridsin points of an A list are also found in a B list

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dxA, dyA, dxyA
    INTEGER, INTENT(in)                                  :: dxB, dyB, dxyB
    INTEGER, DIMENSION(dxA, dyA), INTENT(in)             :: NpointsA
    INTEGER, DIMENSION(dxB, dyB), INTENT(in)             :: NpointsB
    INTEGER, DIMENSION(dxA, dyA, dxyA, 2), INTENT(in)    :: pointsA
    INTEGER, DIMENSION(dxB, dyB, dxyB, 2), INTENT(in)    :: pointsB
    INTEGER, DIMENSION(dxA, dyA, dxB, dyB), INTENT(out)  :: Npoints
    INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA, 2),                                                  &
      INTENT(out)                                        :: points
    INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA),                                                     &
      INTENT(out)                                        :: inpointsA
    INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA),                                                     &
      INTENT(out)                                        :: inpointsB

    ! Local
    INTEGER                                              :: ixA, iyA, ixB, iyB, iv, ii
    INTEGER                                              :: NA, NB
    INTEGER, DIMENSION(dxyA)                             :: ptsA, ptsB
    INTEGER, DIMENSION(dxyA, 2)                          :: pts


!!!!!!! Variables
! dxA, dyA: 2D shape of the list A
! NpointsA: 2D Number of points of the list A
! pointsA: points of the list A
! dxB, dyB: 2D shape of the list B
! NpointsB: 2D Number of points of the list B
! pointsB: points of the list B
! Npoints: Number of coincident points
! points: coincident points
! inpointsA: coincident points list A
! inpointsB: coincident points list B

    fname = 'coincident_gridsin2D_old'

    Npoints = 0
    points = 0
    inpointsA = 0
    inpointsB = 0

    DO ixA=1, dxA
      DO iyA=1, dyA
        NA = NpointsA(ixA,iyA)
        DO ixB=1, dxB
          DO iyB=1, dyB
            NB = NpointsB(ixB,iyB)
            pts = -1
            CALL coincident_list_2Dcoords(NA, pointsA(ixA,iyA,1:NA,:), NB, pointsB(ixB,iyB,1:NB,:),   &
              Npoints(ixA,iyA,ixB,iyB), pts(1:NA,:), ptsA, ptsB)
            DO iv = 1, Npoints(ixA,iyA,ixB,iyB)
              points(ixA,iyA,ixB,iyB,iv,1) = pts(iv,1)
              points(ixA,iyA,ixB,iyB,iv,2) = pts(iv,2)
              inpointsA(ixA,iyA,ixB,iyB,iv) = ptsA(iv)
              inpointsB(ixA,iyA,ixB,iyB,iv) = ptsB(iv)
            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE coincident_gridsin2D_old

  SUBROUTINE coincident_gridsin2D(dxA, dyA, dxyA, NpointsA, pointsA, dxB, dyB, dxyB, NpointsB,        &
    pointsB, Npoints, points, inpointsA, inpointsB)
  ! Subroutine to determine which lists of 2D gridsin points of an A list are also found in a B list

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dxA, dyA, dxyA
    INTEGER, INTENT(in)                                  :: dxB, dyB, dxyB
    INTEGER, DIMENSION(dxA, dyA), INTENT(in)             :: NpointsA
    INTEGER, DIMENSION(dxB, dyB), INTENT(in)             :: NpointsB
    INTEGER, DIMENSION(dxA, dyA, dxyA, 2), INTENT(in)    :: pointsA
    INTEGER, DIMENSION(dxB, dyB, dxyB, 2), INTENT(in)    :: pointsB
    INTEGER, DIMENSION(dxA, dyA, dxB, dyB), INTENT(out)  :: Npoints
    INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA, 2),                                                  &
      INTENT(out)                                        :: points
    INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA),                                                     &
      INTENT(out)                                        :: inpointsA
    INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA),                                                     &
      INTENT(out)                                        :: inpointsB

    ! Local
    INTEGER                                              :: ixA, iyA, ixB, iyB, iv, iv1, iv2
    INTEGER                                              :: NA, NB
    INTEGER, DIMENSION(dxyA)                             :: ptsA, ptsB
    INTEGER, DIMENSION(dxyA, 2)                          :: pts


!!!!!!! Variables
! dxA, dyA: 2D shape of the list A
! NpointsA: 2D Number of points of the list A
! pointsA: points of the list A
! dxB, dyB: 2D shape of the list B
! NpointsB: 2D Number of points of the list B
! pointsB: points of the list B
! Npoints: Number of coincident points
! points: coincident points
! inpointsA: coincident points list A
! inpointsB: coincident points list B

    fname = 'coincident_gridsin2D'

    Npoints = 0
    points = 0
    inpointsA = 0
    inpointsB = 0

    DO ixA=1, dxA
      DO iyA=1, dyA
        NA = NpointsA(ixA,iyA)
        DO ixB=1, dxB
          DO iyB=1, dyB
            NB = NpointsB(ixB,iyB)
            iv = 0
            DO iv1=1, NA
              DO iv2=1, NB
                IF ( (pointsA(ixA,iyA,iv1,1) == pointsB(ixB,iyB,iv2,1)) .AND.                         &
                  (pointsA(ixA,iyA,iv1,2) == pointsB(ixB,iyB,iv2,2)) ) THEN
                  iv = iv + 1
                  points(ixA,iyA,ixB,iyB,iv,1) = pointsA(ixA,iyA,iv1,1)
                  points(ixA,iyA,ixB,iyB,iv,2) = pointsA(ixA,iyA,iv1,2)
                  inpointsA(ixA,iyA,ixB,iyB,iv) = iv1
                  inpointsB(ixA,iyA,ixB,iyB,iv) = iv2
                END IF
              END DO
            END DO
            Npoints(ixA,iyA,ixB,iyB) = iv
          END DO
        END DO
      END DO
    END DO    

  END SUBROUTINE coincident_gridsin2D

END MODULE module_scientific