Changeset 2528 in lmdz_wrf

Timestamp:

May 11, 2019, 8:29:23 PM (6 years ago)

Author:

lfita

Message:

Adding look for coincidence of the faces of the grid and the polygon in `grid_within_polygon'

File:

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

trunk/tools/module_scientific.f90

@@ -4936 +4936 @@
 
 ! Local
-    INTEGER                                              :: ix, iy, iv
-    REAL(r_k), DIMENSION(2)                              :: centergrid, vertex
+    INTEGER                                              :: ix, iy, iv, ivA
+    INTEGER                                              :: cross
+    REAL(r_k), DIMENSION(2)                              :: centergrid, vertex, crosspoint
+    REAL(r_k), DIMENSION(2,2)                            :: face, faceA
     LOGICAL, DIMENSION(dx,dy)                            :: within
 
 !!!!!!! Variables
-! NvertexA: Number of vertices of the polygin to find the grids
+! NvertexA: Number of vertices of the polygon to find the grids
 ! polygonA: ordered vertices of the polygon
 ! dx, dy: shape of the matrix with the grid points
@@ -4956 +4958 @@
     gridsin = 0
     within = .FALSE.
+
     DO ix = 1, dx
       DO iy = 1, dy
@@ -4987 +4990 @@
           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
         END IF
       END DO