Changeset 2288 in lmdz_wrf for trunk/tools

Timestamp:

Jan 25, 2019, 3:30:00 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `grid_spacepercen': Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another series of grid-cells (A)
• `grid_spacepercen_within_reg': Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed by a polygon

File:

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

trunk/tools/module_scientific.f90

@@ -20 +20 @@
 ! fill3DR_2Dvec: Subroutine to fill a 3D float matrix from a series of indices from a 2D matrix
 ! grid_within_polygon: Subroutine to determine which grid cells from a matrix lay inside a polygon
+! grid_spacepercen: Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another 
+!   series of grid-cells (A)
+! grid_spacepercen_within_reg: Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed 
+!   by a polygon
 ! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon 
 !   following ray casting algorithm
@@ -4941 +4945 @@
     END IF
 
+    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
+    IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
+    IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly)
+    IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly)
+    IF (ALLOCATED(poly)) DEALLOCATE(poly)
+
   END SUBROUTINE spacepercen_within_reg
 
+  SUBROUTINE grid_spacepercen_within_reg(NvertexA, polygonA, dx, dy, Nvertexmax, xBvals, yBvals,      &
+    Ngridsin, gridsin, strict, gridspace, percentages)
+  ! Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed 
+  !   by a polygon
+  ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x 
+  !   and y axes. 
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: NvertexA, dx, dy, Nvertexmax
+    REAL(r_k), DIMENSION(NvertexA,2), INTENT(in)         :: polygonA
+    REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in)   :: xBvals, yBvals
+    INTEGER, INTENT(in)                                  :: Ngridsin
+    INTEGER, DIMENSION(Ngridsin,2), INTENT(in)           :: gridsin
+    LOGICAL, INTENT(in)                                  :: strict
+    REAL(r_k), DIMENSION(Ngridsin), INTENT(out)          :: gridspace, percentages
+
+! Local
+   INTEGER                                               :: ig, iv, ix, iy
+   INTEGER                                               :: Nvertex, NvertexAgrid, Ncoin, Nsort
+   CHARACTER(len=20)                                     :: DS
+   REAL(r_k)                                             :: areapoly, areagpoly
+   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                :: vertexgrid, icoinpoly, coinpoly,          &
+     sortpoly, poly
+
+!!!!!!! Variables
+! 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
+! xCvals, yCvals: coordinates of the center of the grid cells
+! Nvertexmax: Maximum number of vertices of the grid cells
+! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex)
+! Ngridsin: number of grids with some extension within the polygon
+! gridsin: grids within the polygon
+! strict: give an error if the area of the polygon is not fully covered
+! gridspace: area of each of the grids
+! percentages: percentages of grid area of each of the grids within the polygon
+
+    fname = 'grid_spacepercen_within_reg'
+
+    gridspace = zeroRK
+    percentages = zeroRK
+
+    DO ig = 1, Ngridsin
+      ix = gridsin(ig,1)
+      iy = gridsin(ig,2) 
+
+     ! Getting grid vertices
+      Nvertex = 0
+      DO iv=1, Nvertexmax
+        IF (xBvals(ix,iy,iv) /= fillvalI) THEN
+          Nvertex = Nvertex + 1
+        END IF
+      END DO
+      IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
+      ALLOCATE(vertexgrid(Nvertex,2))
+      vertexgrid(:,1) = xBvals(ix,iy,1:Nvertex)
+      vertexgrid(:,2) = yBvals(ix,iy,1:Nvertex)
+      areapoly = shoelace_area_polygon(Nvertex, vertexgrid)
+
+      ! Getting common vertices
+      NvertexAgrid = NvertexA*Nvertex
+      IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
+      ALLOCATE(icoinpoly(NvertexAgrid,2))
+      CALL coincident_polygon(NvertexA, Nvertex, NvertexAgrid, polygonA, vertexgrid, Ncoin, icoinpoly)
+
+      IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly)
+      ALLOCATE(coinpoly(Ncoin,2))
+      DO iv=1, Ncoin
+        coinpoly(iv,:) = icoinpoly(iv,:)
+      END DO
+
+      IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly)
+      ALLOCATE(sortpoly(Ncoin,2))
+      CALL sort_polygon(Ncoin, coinpoly, 1, Nsort, sortpoly)
+
+      IF (ALLOCATED(poly)) DEALLOCATE(poly)
+      ALLOCATE(poly(Nsort,2))
+      DO iv=1, Nsort
+        poly(iv,:) = sortpoly(iv,:)
+      END DO
+
+      areagpoly = shoelace_area_polygon(Nsort, poly)
+      gridspace(ig) = ABS(areapoly)
+      percentages(ig) = ABS(areagpoly / areapoly)
+    END DO
+
+    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
+    IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
+    IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly)
+    IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly)
+    IF (ALLOCATED(poly)) DEALLOCATE(poly)
+
+  END SUBROUTINE grid_spacepercen_within_reg
+
   SUBROUTINE spacepercen(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals, xBBvals, yBBvals,      &
-    dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Ngridsin, gridsin, percentages)
+    dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Ngridsin, gridsin, areas, percentages)
   ! Subroutine to compute the space-percentages of a series of grid cells (B) into another series of 
   !   grid-cells (A)
@@ -4961 +5066 @@
     INTEGER, DIMENSION(dxA,dyA), INTENT(out)             :: Ngridsin
     INTEGER, DIMENSION(dxA,dyA,dxyB,2), INTENT(out)      :: gridsin
+    REAL(r_k), DIMENSION(dxA,dyA), INTENT(out)           :: areas
     REAL(r_k), DIMENSION(dxA,dyA,dxyB), INTENT(out)      :: percentages
 
@@ -4967 +5073 @@
    INTEGER                                               :: Nvertex
    INTEGER, ALLOCATABLE, DIMENSION(:,:)                  :: poinsin
-   CHARACTER(len=20)                                     :: DS
+   CHARACTER(len=20)                                     :: IS
    REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                :: vertexgrid
 
@@ -4982 +5088 @@
 ! Ngridsin: number of grids from grid B with some extension within the grid cell A
 ! gridsin: indices of B grids within the grids of A
+! areas: areas of the polygons
 ! percentages: percentages of area of cells B of each of the grids within the grid cell A
 
@@ -4989 +5096 @@
       DO iy = 1, dyA
 
-      ! Getting grid vertices
-      Nvertex = 0
-      DO iv=1, NAvertexmax
-        IF (xBAvals(ix,iy,iv) /= fillval64) THEN
-          Nvertex = Nvertex + 1
-        END IF
+        ! Getting grid vertices
+        Nvertex = 0
+        DO iv=1, NAvertexmax
+          IF (xBAvals(ix,iy,iv) /= fillval64) THEN
+           Nvertex = Nvertex + 1
+          END IF
+        END DO
+        IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
+        ALLOCATE(vertexgrid(Nvertex,2))
+        vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex)
+        vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex)
+  
+        CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xCBvals, yCBvals,            &
+          NBvertexmax, xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,:,:))
+   
+        IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)
+        ALLOCATE(poinsin(Ngridsin(ix,iy),2))
+
+        DO iv=1, Ngridsin(ix,iy)
+          poinsin(iv,1) = gridsin(ix,iy,iv,1)
+          poinsin(iv,2) = gridsin(ix,iy,iv,2)
+        END DO
+
+        areas(ix,iy) = shoelace_area_polygon(Nvertex, vertexgrid)
+        CALL spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals, yBBvals,     &
+          Ngridsin(ix,iy), poinsin, strict, percentages(ix,iy,:))
+
       END DO
-      IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
-      ALLOCATE(vertexgrid(Nvertex,2))
-      vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex)
-      vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex)
+    END DO
+
+    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
+    IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)
+
+  END SUBROUTINE spacepercen
+
+  SUBROUTINE grid_spacepercen(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals, xBBvals, yBBvals, &
+    dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Ngridsin, gridsin, areas, percentages)
+  ! Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another 
+  !   series of grid-cells (A)
+  ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x 
+  !   and y axes. 
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dxA, dyA, NAvertexmax
+    INTEGER, INTENT(in)                                  :: dxB, dyB, NBvertexmax, dxyB
+    REAL(r_k), DIMENSION(dxA,dyA), INTENT(in)            :: xCAvals, yCAvals
+    REAL(r_k), DIMENSION(dxB,dyB), INTENT(in)            :: xCBvals, yCBvals
+    REAL(r_k), DIMENSION(dxA,dyA,NAvertexmax), INTENT(in):: xBAvals, yBAvals
+    REAL(r_k), DIMENSION(dxB,dyB,NBvertexmax), INTENT(in):: xBBvals, yBBvals
+    LOGICAL, INTENT(in)                                  :: strict
+    INTEGER, DIMENSION(dxA,dyA), INTENT(out)             :: Ngridsin
+    INTEGER, DIMENSION(dxA,dyA,dxyB,2), INTENT(out)      :: gridsin
+    REAL(r_k), DIMENSION(dxB,dyB), INTENT(out)           :: areas
+    REAL(r_k), DIMENSION(dxA,dyA,dxyB), INTENT(out)      :: percentages
+
+! Local
+   INTEGER                                               :: iv, ix, iy
+   INTEGER                                               :: Nvertex
+   INTEGER, ALLOCATABLE, DIMENSION(:,:)                  :: poinsin
+   CHARACTER(len=20)                                     :: IS
+   REAL(r_k), ALLOCATABLE, DIMENSION(:)                  :: pareas
+   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                :: vertexgrid
+
+!!!!!!! Variables
+! dxA, dyA: shape of the matrix with the grid points A
+! xCAvals, yCAvals: coordinates of the center of the grid cells A
+! NAvertexmax: Maximum number of vertices of the grid cells A
+! xBAvals, yBAvals: coordinates of th vertices of the grid cells A (-99999 for no vertex)
+! dxB, dyB: shape of the matrix with the grid points B
+! xCBvals, yCBvals: coordinates of the center of the grid cells B
+! NBvertexmax: Maximum number of vertices of the grid cells B
+! xBBvals, yBBvals: coordinates of th vertices of the grid cells B (-99999 for no vertex)
+! strict: give an error if the area of the polygon is not fully covered
+! Ngridsin: number of grids from grid B with some extension within the grid cell A
+! gridsin: indices of B grids within the grids of A
+! areas: areas of the grids
+! percentages: percentages of area of cells B of each of the grids inside the grid cell A
+
+    fname = 'grid_spacepercen'
+
+    areas = zeroRK
+
+    DO ix = 1, dxA
+      DO iy = 1, dyA
+
+        ! Getting grid vertices
+        Nvertex = 0
+        DO iv=1, NAvertexmax
+          IF (xBAvals(ix,iy,iv) /= fillval64) THEN
+           Nvertex = Nvertex + 1
+          END IF
+        END DO
+        IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
+        ALLOCATE(vertexgrid(Nvertex,2))
+        vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex)
+        vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex)
   
-      CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xCBvals, yCBvals, NBvertexmax, &
-        xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,:,:))
+        CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xCBvals, yCBvals,            &
+          NBvertexmax, xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,:,:))
    
-      IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)
-      ALLOCATE(poinsin(Ngridsin(ix,iy),2))
-
-      DO iv=1, Ngridsin(ix,iy)
-        poinsin(iv,1) = gridsin(ix,iy,iv,1)
-        poinsin(iv,2) = gridsin(ix,iy,iv,2)
+        IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)
+        ALLOCATE(poinsin(Ngridsin(ix,iy),2))
+        IF (ALLOCATED(pareas)) DEALLOCATE(pareas)
+        ALLOCATE(pareas(Ngridsin(ix,iy)))
+
+        DO iv=1, Ngridsin(ix,iy)
+          poinsin(iv,1) = gridsin(ix,iy,iv,1)
+          poinsin(iv,2) = gridsin(ix,iy,iv,2)
+        END DO
+
+        CALL grid_spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals,       &
+          yBBvals, Ngridsin(ix,iy), poinsin, strict, pareas, percentages(ix,iy,:))
+
+        ! Filling areas
+        DO iv = 1, Ngridsin(ix,iy)
+          IF (areas(poinsin(iv,1), poinsin(iv,2)) == zeroRK) THEN
+            areas(poinsin(iv,1), poinsin(iv,2)) = pareas(iv)
+          END IF
+        END DO
+
       END DO
-
-      CALL spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals, yBBvals,       &
-        Ngridsin(ix,iy), poinsin, strict, percentages(ix,iy,:))
-
     END DO
-    END DO
-
-  END SUBROUTINE spacepercen
+
+    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
+    IF (ALLOCATED(pareas)) DEALLOCATE(pareas)
+    IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)
+
+  END SUBROUTINE grid_spacepercen
 
   SUBROUTINE unique_matrixRK2D(dx, dy, dxy, matrix2D, Nunique, unique)
@@ -5040 +5246 @@
 ! dxy: dx*dy, maximum possible amount of different values
 ! matrix2D: matgrix of values
-! Nunique: amount of unique values
+! Nunique: amount of unique values 
 ! unique: sorted from minimum to maximum vector with the unique values