Index: trunk/tools/module_scientific.f90
===================================================================
--- trunk/tools/module_scientific.f90	(revision 2553)
+++ trunk/tools/module_scientific.f90	(revision 2554)
@@ -4941,4 +4941,5 @@
     REAL(r_k), DIMENSION(2,2)                            :: face, faceA
     LOGICAL, DIMENSION(dx,dy)                            :: within
+    REAL(r_k), DIMENSION(:,:), ALLOCATABLE               :: gridpoly
 
 !!!!!!! Variables
@@ -4990,82 +4991,93 @@
           END DO
 
-          ! Getting grid faces
-          DO iv=1, Nvertexmax-1
-            IF (.NOT.within(ix,iy)) THEN
-              IF (xBvals(ix,iy,iv) /= fillvalI) THEN
-                face(1,:) = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /)
-                face(2,:) = (/ xBvals(ix,iy,iv+1), yBvals(ix,iy,iv+1) /)
-                DO ivA=1, NvertexA-1
-                  faceA(1,:) = (/ polygonA(ivA,1), polygonA(ivA,2) /)
-                  faceA(2,:) = (/polygonA(ivA+1,1), polygonA(ivA+1,2) /)
-                  CALL intersectfaces(face, faceA, cross, crosspoint)
-                  IF (cross /= 0) THEN
-                    Ngridsin = Ngridsin + 1
-                    ! Getting coordinates
-                    gridsin(Ngridsin,1) = ix
-                    gridsin(Ngridsin,2) = iy
-                    within(ix,iy) = .TRUE.
-                    CYCLE
-                  END IF
-                END DO
-                IF (.NOT.within(ix,iy)) THEN
-                  faceA(1,:) = (/ polygonA(NvertexA,1), polygonA(NvertexA,2) /)
-                  faceA(2,:) = (/polygonA(1,1), polygonA(1,2) /)
-                  CALL intersectfaces(face, faceA, cross, crosspoint)
-                  IF (cross /= 0) THEN
-                    Ngridsin = Ngridsin + 1
-                    ! Getting coordinates
-                    gridsin(Ngridsin,1) = ix
-                    gridsin(Ngridsin,2) = iy
-                    within(ix,iy) = .TRUE.
-                  END IF
-                END IF
-              END IF
-            END IF
-          END DO
-          iv = Nvertexmax
-          IF (.NOT.within(ix,iy)) THEN
-            IF (xBvals(ix,iy,iv) /= fillvalI) THEN
-              face(1,:) = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /)
-              face(2,:) = (/ xBvals(ix,iy,1), yBvals(ix,iy,1) /)
-              DO ivA=1, NvertexA-1
-                faceA(1,:) = (/ polygonA(ivA,1), polygonA(ivA,2) /)
-                faceA(2,:) = (/polygonA(ivA+1,1), polygonA(ivA+1,2) /)
-                CALL intersectfaces(face, faceA, cross, crosspoint)
-                IF (cross /= 0) THEN
-                  Ngridsin = Ngridsin + 1
-                  ! Getting coordinates
-                  gridsin(Ngridsin,1) = ix
-                  gridsin(Ngridsin,2) = iy
-                  within(ix,iy) = .TRUE.
-                  CYCLE
-                END IF
-              END DO
-              IF (.NOT.within(ix,iy)) THEN
-                faceA(1,:) = (/ polygonA(NvertexA,1), polygonA(NvertexA,2) /)
-                faceA(2,:) = (/polygonA(1,1), polygonA(1,2) /)
-                CALL intersectfaces(face, faceA, cross, crosspoint)
-                IF (cross /= 0) THEN
-                  Ngridsin = Ngridsin + 1
-                  ! Getting coordinates
-                  gridsin(Ngridsin,1) = ix
-                  gridsin(Ngridsin,2) = iy
-                  within(ix,iy) = .TRUE.
-                END IF
-              END IF
-            END IF
-          END IF
+!          ! Getting grid faces
+!          DO iv=1, Nvertexmax-1
+!            IF (.NOT.within(ix,iy)) THEN
+!              IF (xBvals(ix,iy,iv) /= fillvalI) THEN
+!                face(1,:) = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /)
+!                face(2,:) = (/ xBvals(ix,iy,iv+1), yBvals(ix,iy,iv+1) /)
+!                DO ivA=1, NvertexA-1
+!                  faceA(1,:) = (/ polygonA(ivA,1), polygonA(ivA,2) /)
+!                  faceA(2,:) = (/polygonA(ivA+1,1), polygonA(ivA+1,2) /)
+!                  CALL intersectfaces(face, faceA, cross, crosspoint)
+!                  IF (cross /= 0) THEN
+!                    Ngridsin = Ngridsin + 1
+!                    ! Getting coordinates
+!                    gridsin(Ngridsin,1) = ix
+!                    gridsin(Ngridsin,2) = iy
+!                    within(ix,iy) = .TRUE.
+!                    CYCLE
+!                  END IF
+!                END DO
+!                IF (.NOT.within(ix,iy)) THEN
+!                  faceA(1,:) = (/ polygonA(NvertexA,1), polygonA(NvertexA,2) /)
+!                  faceA(2,:) = (/polygonA(1,1), polygonA(1,2) /)
+!                  CALL intersectfaces(face, faceA, cross, crosspoint)
+!                  IF (cross /= 0) THEN
+!                    Ngridsin = Ngridsin + 1
+!                    ! Getting coordinates
+!                    gridsin(Ngridsin,1) = ix
+!                    gridsin(Ngridsin,2) = iy
+!                    within(ix,iy) = .TRUE.
+!                  END IF
+!                END IF
+!              END IF
+!            END IF
+!          END DO
+!          iv = Nvertexmax
+!          IF (.NOT.within(ix,iy)) THEN
+!            IF (xBvals(ix,iy,iv) /= fillvalI) THEN
+!              face(1,:) = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /)
+!              face(2,:) = (/ xBvals(ix,iy,1), yBvals(ix,iy,1) /)
+!              DO ivA=1, NvertexA-1
+!                faceA(1,:) = (/ polygonA(ivA,1), polygonA(ivA,2) /)
+!                faceA(2,:) = (/polygonA(ivA+1,1), polygonA(ivA+1,2) /)
+!                CALL intersectfaces(face, faceA, cross, crosspoint)
+!                IF (cross /= 0) THEN
+!                  Ngridsin = Ngridsin + 1
+!                  ! Getting coordinates
+!                  gridsin(Ngridsin,1) = ix
+!                  gridsin(Ngridsin,2) = iy
+!                  within(ix,iy) = .TRUE.
+!                  CYCLE
+!                END IF
+!              END DO
+!              IF (.NOT.within(ix,iy)) THEN
+!                faceA(1,:) = (/ polygonA(NvertexA,1), polygonA(NvertexA,2) /)
+!                faceA(2,:) = (/polygonA(1,1), polygonA(1,2) /)
+!                CALL intersectfaces(face, faceA, cross, crosspoint)
+!                IF (cross /= 0) THEN
+!                  Ngridsin = Ngridsin + 1
+!                  ! Getting coordinates
+!                  gridsin(Ngridsin,1) = ix
+!                  gridsin(Ngridsin,2) = iy
+!                  within(ix,iy) = .TRUE.
+!               END IF
+!              END IF
+!            END IF
+!          END IF
 
           ! Inverting process, look polygon within grid
           IF (.NOT.within(ix,iy)) THEN
-            DO iv = 1, NvertexA
-              vertex = polygonA(iv,:)
-              IF (point_inside(vertex, Nvertexmax, xBvals(ix,iy,:))) THEN
-                Ngridsin = Ngridsin + 1
-                ! Getting coordinates
-                gridsin(Ngridsin,1) = ix
-                gridsin(Ngridsin,2) = iy
-                within(ix,iy) = .TRUE.
-                CYCLE
+            IF (ALLOCATED(gridpoly)) DEALLOCATE(gridpoly)
+            ALLOCATE(gridpoly(Nvertexmax,2))
+            DO iv = 1, Nvertexmax
+              gridpoly(iv,1) = xBvals(ix,iy,iv)
+              gridpoly(iv,2) = yBvals(ix,iy,iv)
+            END DO
+            print *,ix,iy,'gridpoly :', gridpoly
+            DO ivA = 1, NvertexA
+              IF (.NOT.within(ix,iy)) THEN
+                 vertex = polygonA(ivA,:)
+                 print *,'  ', ivA, 'vertex :', vertex
+                 PRINT *,'  ',point_inside(vertex, Nvertexmax, gridpoly)
+                 IF (point_inside(vertex, Nvertexmax, gridpoly)) THEN
+                   Ngridsin = Ngridsin + 1
+                   ! Getting coordinates
+                   gridsin(Ngridsin,1) = ix
+                   gridsin(Ngridsin,2) = iy
+                   within(ix,iy) = .TRUE.
+                   CYCLE
+                 END IF
               END IF
             END DO
@@ -5074,4 +5086,6 @@
       END DO
     END DO
+
+    IF (ALLOCATED(gridpoly)) DEALLOCATE(gridpoly)
 
   END SUBROUTINE grid_within_polygon
