Changeset 1612 in lmdz_wrf for trunk

Timestamp:

Aug 30, 2017, 2:54:25 PM (8 years ago)

Author:

lfita

Message:

Adding all the required Fortran code for the path/polygon detection an definition

Location:

trunk/tools

Files:

• module_definitions.f90 (modified) (1 diff)
• module_generic.f90 (modified) (2 diffs)
• module_scientific.f90 (modified) (2 diffs)

trunk/tools/module_definitions.f90

@@ -13 +13 @@
   REAL(r_k)                                              :: fillval64 = 1.e20
   CHARACTER(len=50)                                      :: fname
+  CHARACTER(len=256)                                     :: msg
 
 END MODULE module_definitions

trunk/tools/module_generic.f90

@@ -5 +5 @@
 ! ErrMsg: Subroutine to stop execution and provide an error message
 ! ErrWarnMsg: Function to print error/warning message
+! index_list_coordsI: Function to provide the index of a given coordinate within a list of integer coordinates
 ! Index1DArrayR: Function to provide the first index of a given value inside a 1D real array
 ! Index1DArrayR_K: Function to provide the first index of a given value inside a 1D real(r_k) array
@@ -17 +18 @@
 
   CONTAINS
+
+  INTEGER FUNCTION index_list_coordsI(Ncoords, coords, icoord)
+  ! Function to provide the index of a given coordinate within a list of integer coordinates
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: Ncoords
+    INTEGER, DIMENSION(Ncoords,2), INTENT(in)            :: coords
+    INTEGER, DIMENSION(2), INTENT(in)                    :: icoord
+
+! Local
+    INTEGER, DIMENSION(Ncoords)                          :: dist
+    INTEGER                                              :: i,mindist
+    INTEGER, DIMENSION(1)                                :: iloc
+
+!!!!!!! Variables
+! Ncoords: number of coordinates in the list
+! coords: list of coordinates
+! icoord: coordinate to find
+
+  fname = 'index_list_coordsI'
+
+  dist = (coords(:,1)-icoord(1))**2+(coords(:,2)-icoord(2))**2
+
+  IF (ANY(dist == 0)) THEN
+    iloc = MINLOC(dist)
+    index_list_coordsI = iloc(1)
+  ELSE
+    index_list_coordsI = -1
+  END IF
+
+  END FUNCTION index_list_coordsI
 
   CHARACTER(len=1000) FUNCTION vectorR_KS(d1, vector)

trunk/tools/module_scientific.f90

@@ -3 +3 @@
   
 !!!!!!! Functions & subroutines
+! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.)
 ! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
+! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon 
+!   following ray casting algorithm
+! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders 
+!   (limits of a region)
+! paths_border: Subroutine to search the paths of a border field.
+! path_properties: Subroutine to determine the properties of a path
+! 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!
 ! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
 ! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k'
@@ -17 +26 @@
 
   CONTAINS
+
+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
+
+  DO it=1,dt
+    IF (dbg) THEN
+      PRINT *,'  it:', it, ' bool _______'
+      DO i=1,dx
+        PRINT *,boolmatt(:,:,it)
+      END DO
+    END IF
+    CALL polygons(dbg, dx, dy, boolmatt(:,:,it), polys(:,:,it), Npoly(it))
+  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
+  INTEGER, DIMENSION(:,:,:), ALLOCATABLE                 :: paths
+  INTEGER                                                :: Npath
+  INTEGER, DIMENSION(:), ALLOCATABLE                     :: Nptpaths
+  INTEGER, DIMENSION(2)                                  :: xtrx, xtry, meanpth
+  INTEGER                                                :: Nvertx
+  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
+
+  Nppt = dx*dy
+  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
+  IF (ALLOCATED(points)) DEALLOCATE(points)
+  ALLOCATE(points(Nppt,2), STAT=ierr)
+  msg = "Problems allocating matrix 'points'"
+  CALL ErrMsg(msg, fname, ierr)
+
+  ip = 1
+  DO i=1, dx
+    DO j=1, dy
+      points(ip,1) = i
+      points(ip,2) = j
+      ip = ip + 1
+    END DO
+  END DO
+
+  CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder)
+  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(Nppt), 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, boolmat, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,  &
+      xtrx, xtry, vertxs, Nppt, points, isin)
+
+    ! Filling polygons
+    ipp = 1
+    DO i=1, dx
+      DO j=1, dy
+        IF (isin(ipp)) polys(i,j) = ip
+        ipp = ipp + 1
+      END DO
+    END DO
+
+  END DO
+
+  DEALLOCATE (borders)  
+  DEALLOCATE (Nptpaths)
+  DEALLOCATE (paths)
+  DEALLOCATE (vertxs)
+  DEALLOCATE (points)
+  DEALLOCATE (isin)
+
+  RETURN
+
+END SUBROUTINE 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
+  INTEGER                                                :: ierr
+  LOGICAL, DIMENSION(:,:), ALLOCATABLE                   :: halo_brdr
+
+!!!!!!! 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
+! 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?
+      IF (isbrdr(ix,iy)) THEN
+        inside(ip) = .TRUE.
+        CYCLE
+      END IF
+
+      ! Looking along y-axis
+      DO j=1,iy
+        ! Only counting that borders that are not vertexs
+        isvertex = index_list_coordsI(Npath, vrtxs, (/ix,j/))
+        IF (halo_brdr(ix+1,j+1) .AND. ( isvertex == -1)) Nintersecs = Nintersecs + 1 
+      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 = (/ 'N ', 'NE', 'E ', 'SE', 'S ', 'SW', 'W ', 'NW' /)
+  ! Transformation to apply
+  !spt = (/ (/0,1/), (/1,1/), (/1,0/), (/1,-1/), (/0,-1/), (/-1,-1/), (/-1,0/), (/-1,1/) /)
+  spt(:,1) = (/ 0, 1, 1, 1, 0, -1, -1, -1 /)
+  spt(:,2) = (/ 1, 1, 0, -1, -1, -1, 0, 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)
+! 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
+
+! 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
+
+  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
+
+  ! 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.
+          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.
+          ib=ib+1
+        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 *,'  isborder ______'
+    DO i=1,dx
+      PRINT *,isborder(i,:)
+    END DO
+
+    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 next point would be someone already in the path
+            CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE.,iif,jjf,  &
+              iff)
+            iff = index_list_coordsI(Nptspath, paths(ipath,:,:),(/iif,jjf/))
+            IF (iff /= -1) THEN
+              IF (dbg) THEN
+                PRINT *,'    point already in Path!'
+                PRINT *,'  path:', ipath, '_____'
+                DO i=1, Nptspath
+                  PRINT *,'    ',i,':',paths(ipath,i,:)
+                END DO
+              END IF
+              ! Continung unfinished path
+              iip = paths(ipath,Nptspath,1)
+              ijp = paths(ipath,Nptspath,2)
+              iipth = index_list_coordsI(Nppt, borders, (/iip,ijp/))
+            ELSE
+              ! 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)