Changeset 1616 in lmdz_wrf

Timestamp:

Sep 3, 2017, 12:29:48 AM (8 years ago)

Author:

lfita

Message:

Adding:

• `clean_polygons': Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole

File:

• trunk/tools/module_scientific.f90 (modified) (7 diffs)

trunk/tools/module_scientific.f90

@@ -4 +4 @@
 !!!!!!! Functions & subroutines
 ! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.)
+! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
 ! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
 ! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon 
@@ -107 +108 @@
 
   Nppt = dx*dy/10
-  Npts = dx*dy
 
   IF (ALLOCATED(borders)) DEALLOCATE(borders)
@@ -124 +124 @@
   CALL ErrMsg(msg, fname, ierr)
 
-  ! Filling with the points of all the space
+  ! Filling with the points of all the space with .TRUE.
+  Npts = COUNT(boolmat)
+
   IF (ALLOCATED(points)) DEALLOCATE(points)
   ALLOCATE(points(Npts,2), STAT=ierr)
@@ -130 +132 @@
   CALL ErrMsg(msg, fname, ierr)
 
+  ! We only want to localize that points 'inside'
   ip = 1
   DO i=1, dx
     DO j=1, dy
-      points(ip,1) = i
-      points(ip,2) = j
-      ip = ip + 1
+      IF (boolmat(i,j)) THEN
+        points(ip,1) = i
+        points(ip,2) = j
+        ip = ip + 1
+      END IF
     END DO
   END DO
@@ -177 +182 @@
 
     ! Filling polygons
-    ipp = 1
-    DO i=1, dx
-      DO j=1, dy
-        IF (isin(ipp)) polys(i,j) = ip
-        ipp = ipp + 1
-      END DO
+    DO ipp=1, Npts
+      IF (isin(ipp)) polys(points(ipp,1),points(ipp,2)) = ip
     END DO
 
@@ -196 +197 @@
   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)
@@ -206 +210 @@
 
 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)