Changeset 1614 in lmdz_wrf

Timestamp:

Sep 1, 2017, 12:21:34 AM (7 years ago)

Author:

lfita

Message:

Fixing last issues on polygons detection
Adding direct acces to the 'module_scientific.f90' as `module_ForSci'

Location:

trunk/tools

Files:

• Makefile.llamp (modified) (1 diff)
• module_scientific.f90 (modified) (21 diffs)

trunk/tools/Makefile.llamp

@@ -25 +25 @@
 diagsrcs = module_definitions.f90 module_generic.f90 module_scientific.f90 module_ForDiagnosticsVars.f90 module_ForDiagnostics.f90
 intsrcs = module_definitions.f90 module_generic.f90 module_scientific.f90 module_ForInterpolate.f90
-scisrcs = module_definitions.f90 module_generic.f90 module_scientific.f90 module_ForInterpolate.f90
+scisrcs = module_definitions.f90 module_generic.f90 module_scientific.f90
 
 ####### ###### ##### #### ### ## #

trunk/tools/module_scientific.f90

@@ -52 +52 @@
 
   polys = -1
+  Npoly = 0
 
   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))
+    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
 
@@ -80 +87 @@
   INTEGER                                                :: ierr
   INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: borders
-  LOGICAL, DIMENSION(dx,dy)                              :: isborder
+  LOGICAL, DIMENSION(dx,dy)                              :: isborder, isbordery
   INTEGER, DIMENSION(:,:,:), ALLOCATABLE                 :: paths
   INTEGER                                                :: Npath
   INTEGER, DIMENSION(:), ALLOCATABLE                     :: Nptpaths
   INTEGER, DIMENSION(2)                                  :: xtrx, xtry, meanpth
-  INTEGER                                                :: Nvertx
+  INTEGER                                                :: Nvertx, Npts
   INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: vertxs, points
   LOGICAL, DIMENSION(:), ALLOCATABLE                     :: isin
@@ -99 +106 @@
   polys = -1
 
-  Nppt = dx*dy
+  Nppt = dx*dy/10
+  Npts = dx*dy
+
   IF (ALLOCATED(borders)) DEALLOCATE(borders)
   ALLOCATE(borders(Nppt,2), STAT=ierr)
@@ -117 +126 @@
   ! Filling with the points of all the space
   IF (ALLOCATED(points)) DEALLOCATE(points)
-  ALLOCATE(points(Nppt,2), STAT=ierr)
+  ALLOCATE(points(Npts,2), STAT=ierr)
   msg = "Problems allocating matrix 'points'"
   CALL ErrMsg(msg, fname, ierr)
@@ -130 +139 @@
   END DO
 
-  CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder)
+  CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery)
   CALL paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
 
@@ -142 +151 @@
 
     IF (ALLOCATED(isin)) DEALLOCATE(isin)
-    ALLOCATE(isin(Nppt), STAT=ierr)
+    ALLOCATE(isin(Npts), STAT=ierr)
     msg = "Problems allocating matrix 'isin'"
     CALL ErrMsg(msg, fname, ierr)
@@ -164 +173 @@
     END IF
  
-    CALL gridpoints_InsidePolygon(dx, dy, boolmat, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,  &
-      xtrx, xtry, vertxs, Nppt, points, isin)
+    CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,&
+      xtrx, xtry, vertxs, Npts, points, isin)
 
     ! Filling polygons
@@ -175 +184 @@
       END DO
     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
@@ -306 +324 @@
 ! Local
   INTEGER                                                :: i,j,ip,ix,iy
-  INTEGER                                                :: Nintersecs, isvertex
+  INTEGER                                                :: Nintersecs, isvertex, ispath
   INTEGER                                                :: ierr
   LOGICAL, DIMENSION(:,:), ALLOCATABLE                   :: halo_brdr
+  INTEGER                                                :: Nbrbrdr
 
 !!!!!!! Variables
@@ -316 +335 @@
 ! 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
@@ -345 +365 @@
 
       ! It is a border point?
-      IF (isbrdr(ix,iy)) THEN
+      ispath = index_list_coordsI(Npath, path, (/ix,iy/))
+      IF (isbrdr(ix,iy) .AND. (ispath /= -1)) THEN
         inside(ip) = .TRUE.
         CYCLE
@@ -351 +372 @@
 
       ! Looking along y-axis
-      DO j=1,iy
+      ! 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. ( isvertex == -1)) Nintersecs = Nintersecs + 1 
+
+        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.
@@ -398 +432 @@
 
   ! Label of the search
-  lclock = (/ 'N ', 'NE', 'E ', 'SE', 'S ', 'SW', 'W ', 'NW' /)
+  lclock = (/ 'W ', 'NW', 'N ', 'NE', 'E ', 'SE', 'S ', 'SW' /)
   ! 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 /)
+  !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
@@ -434 +468 @@
 END SUBROUTINE look_clockwise_borders
 
-SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr)
+SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
 ! Subroutine to provide the borders of a logical array (interested in .TRUE.)
 
@@ -442 +476 @@
   LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: Lmat
   INTEGER, DIMENSION(dxy,2), INTENT(out)                 :: brdrs
-  LOGICAL, DIMENSION(dx,dy), INTENT(out)                 :: isbrdr
+  LOGICAL, DIMENSION(dx,dy), INTENT(out)                 :: isbrdr, isbrdry
 
 ! Local
@@ -453 +487 @@
 ! 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'
@@ -491 +526 @@
     END IF
   END DO
+
+  isbrdry = isbrdr
 
   ! Border as that when looking on x-axis points with Lmat(i) /= Lmat(i+1)
@@ -514 +551 @@
           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
@@ -519 +557 @@
           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
@@ -578 +641 @@
     
   IF (dbg) THEN
-    PRINT *,'  isborder ______'
-    DO i=1,dx
-      PRINT *,isborder(i,:)
-    END DO
-
     PRINT *,'    borders _______'
     DO i=1,Nbrdr
@@ -651 +709 @@
           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
+            ! 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
-              ! Continung unfinished path
-              iip = paths(ipath,Nptspath,1)
-              ijp = paths(ipath,Nptspath,2)
-              iipth = index_list_coordsI(Nppt, borders, (/iip,ijp/))
-            ELSE
+            END DO
+            IF (.NOT.found) THEN
               ! Looking for the next available border point for the new path
               DO i=1,Nbrdr