Changeset 2345 in lmdz_wrf for trunk/tools

Timestamp:

Feb 18, 2019, 4:00:24 PM (6 years ago)

Author:

lfita

Message:

Adding right filling value in variable `range' after homogenization

Location:

trunk/tools

Files:

• module_ForDiagnostics.f90 (modified) (1 diff)
• module_scientific.f90 (modified) (5 diffs)

trunk/tools/module_ForDiagnostics.f90

@@ -1056 +1056 @@
     ! Homogenizing indices of the ranges
     CALL continguos_homogene_zones(d1, d2, iranges, Nranges, ranges)
+    WHERE (ranges == -1)
+      ranges = fillValI
+    END WHERE
 
     RETURN

trunk/tools/module_scientific.f90

@@ -2573 +2573 @@
 
 ! Local
-  INTEGER                                                :: i, j, ip, ipp, Nppt
+  INTEGER                                                :: i, j, k, ip, ipp, Nppt
   INTEGER                                                :: ierr
   INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: borders
@@ -2618 +2618 @@
   ! Filling with the points of all the space with .TRUE.
   Npts = COUNT(boolmat)
-  PRINT * ,'Lluis Npts:', Npts, 'Nppt:', Nppt
+  PRINT *, dx, dy, 'Lluis Npts:', Npts, 'Nppt:', Nppt
 
   IF (ALLOCATED(points)) DEALLOCATE(points)
-  ALLOCATE(points(Npts,2), STAT=ierr)
+  ALLOCATE(points(dy,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(dbg, dx, dy, Nppt, boolmat, borders, isborder, isbordery)
@@ -2654 +2642 @@
 
     isin = .FALSE.
-
-    IF (dbg) THEN
-      PRINT *, '  path:', ip, ' N pts:', Nptpaths(ip)
-      DO j=1, Nptpaths(ip)
-        PRINT *, '      ',j,':',paths(ip,j,:)
-      END DO
-    END IF
 
     borderp = .FALSE.
@@ -2681 +2662 @@
     END IF
  
-    CALL gridpoints_InsidePolygon(dbg, 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
+!    CALL gridpoints_InsidePolygon(dbg, dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:),   &
+!      Nvertx, xtrx, xtry, vertxs, Npts, points, isin)
+
+    ! We only want to localize that points 'inside'
+    Npts = 1
+    DO i=xtrx(1), xtrx(2)
+      DO j=1, dy
+        IF (boolmat(i,j)) THEN
+          points(Npts,1) = i
+          points(Npts,2) = j
+          Npts = Npts + 1
+        END IF
+      END DO
+      CALL gridpoints_InsidePolygon_ray(dbg,dy,borderp(i,:), 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 *,'  path boolmat isborder isborderp polygon (', xtrx(1), ',', xtry(1), ')x(', xtrx(2), &
+          ',', xtry(2), ') _______' , Npts
+        PRINT *,ip,'<>', i,':',boolmat(i,xtry(1):xtry(2)), ' border ', isborder(i,xtry(1):xtry(2)),   &
+          ' isbordery ', borderp(i,xtry(1):xtry(2)), ' polygon ', polys(i,xtry(1):xtry(2))
+      END IF
+
     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
+  PRINT *,'  PRE-clean '
 
   ! Cleaning polygons matrix of not-used and paths around holes
-  CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg)
+  !CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg)
 
   IF (ALLOCATED(borders)) DEALLOCATE (borders)
@@ -3022 +3016 @@
 
 END SUBROUTINE gridpoints_InsidePolygon
+
+
+  SUBROUTINE gridpoints_InsidePolygon_ray(dbg, 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)                                    :: dy,Npath,Nvrtx,Npts
+  LOGICAL, INTENT(in)                                    :: dbg
+  LOGICAL, DIMENSION(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(dy,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
+! dy: y-axis 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_ray'
+
+  ! 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(dy+2), STAT=ierr)
+  msg = "Problems allocating matrix 'halo_brdr'"
+  CALL ErrMsg(msg, fname, ierr)
+  halo_brdr = .FALSE.
+
+  IF (dbg) PRINT *,'Border _______'
+  halo_brdr(2:dy+1) = isbrdr(:)
+  IF (dbg) PRINT *,isbrdr(:)
+
+  inside = .FALSE.
+
+  DO ip=1,dy
+    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(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(j+1) .AND. (ispath /= -1) .AND. (isvertex == -1) ) Nintersecs = Nintersecs + 1 
+        IF (halo_brdr(j+1) .AND. (ispath /= -1) .AND. (halo_brdr(j+1) .EQV. halo_brdr(j+2))) THEN
+          Nbrbrdr = Nbrbrdr + 1
+          IF (dbg) PRINT *,'    ',Nbrbrdr,' Consec brdrs:', halo_brdr(j+1), halo_brdr(j+2), '(',      &
+             ix,j,';', ix,j+1,')', '(', ix,j,';', ix,j+1,')', isbrdr(j), isbrdr(j+1) 
+        ELSE
+          ! Will remove that consecutive borders above 2
+          IF (Nbrbrdr /= 0) THEN
+            IF (dbg) PRINT *, ix,',',iy,';', Nintersecs, '  amount of consecutive borders:', Nbrbrdr, &
+              ' removing:', MAX(Nbrbrdr-1, 0)
+            Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0)
+            Nbrbrdr = 0
+          END IF
+        END IF
+      END DO
+      IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE.
+      IF (dbg) PRINT *,ip,'    point:', ix, iy, 'isbrdr:', isbrdr(1:iy-1), 'y-ray:', halo_brdr(1:iy), &
+        'inside:', inside(ip)
+    END IF
+
+  END DO
+
+  RETURN
+
+END SUBROUTINE gridpoints_InsidePolygon_ray
 
 SUBROUTINE look_clockwise_borders(dx,dy,Nbrdrs,brdrs,gbrdr,isbrdr,ix,iy,dbg,xf,yf,iff)