Changeset 2289 in lmdz_wrf for trunk

Timestamp:

Jan 25, 2019, 8:31:33 PM (7 years ago)

Author:

lfita

Message:

Adding:

• `coincident_gridsin2D': Subroutine to determine which lists of 2D gridsin points of an A list are also found in a B list
• `coincident_list_2Dcoords': Subroutine to determine which 2D points of an A list are also found in a B list
• `grid_spacepercen_providing_polys': Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another series of grid-cells (A) providing coincident polygons
• `grid_spacepercen_within_reg_providing_polys': Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed by a polygon providing coordinates of the resultant polygons

File:

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

trunk/tools/module_scientific.f90

@@ -13 +13 @@
 ! coincidence_poly_area: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
 !   to a map of integer polygons filtrered by a minimal area
+! coincident_gridsin2D: Subroutine to determine which lists of 2D gridsin points of an A list are also found in a B list
+! coincident_list_2Dcoords: Subroutine to determine which 2D points of an A list are also found in a B list
 ! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
 ! coincident_polygon: Subroutine to provide the intersection polygon between two polygons
@@ -22 +24 @@
 ! 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_providing_polys: Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another 
+!   series of grid-cells (A) providing coincident polygons
 ! grid_spacepercen_within_reg: Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed 
 !   by a polygon
+! grid_spacepercen_within_reg_providing_polys: Subroutine to compute the percentage of grid space of a series of grid 
+!   cells which are encompassed by a polygon providing coordinates of the resultant polygons
 ! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon 
 !   following ray casting algorithm
@@ -4893 +4899 @@
 
       ! Getting common vertices
-      NvertexAgrid = NvertexA*Nvertex
+      NvertexAgrid = NvertexA*Nvertex*2
       IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
       ALLOCATE(icoinpoly(NvertexAgrid,2))
@@ -5014 +5020 @@
 
       ! Getting common vertices
-      NvertexAgrid = NvertexA*Nvertex
+      NvertexAgrid = NvertexA*Nvertex*2
       IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly)
       ALLOCATE(icoinpoly(NvertexAgrid,2))
@@ -5047 +5053 @@
 
   END SUBROUTINE grid_spacepercen_within_reg
+
+  SUBROUTINE grid_spacepercen_within_reg_providing_polys(NvertexA, polygonA, dx, dy, Nvertexmax,      &
+    xBvals, yBvals, Ngridsin, gridsin, strict, Nmaxver2, Ncoinpoly, ccoinpoly, gridspace, percentages)
+  ! Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed 
+  !   by a polygon providing coordinates of the resultant polygons
+  ! 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, Nmaxver2
+    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
+    INTEGER, DIMENSION(Ngridsin), INTENT(out)            :: Ncoinpoly
+    REAL(r_k), DIMENSION(Ngridsin,Nmaxver2,2),                                                        &
+      INTENT(out)                                        :: ccoinpoly
+    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
+! Nmaxver2: maximum possible number of vertices of the coincident polygon
+! Ncoinpoly: number of vertices of the coincident polygon
+! coinpoly: coordinates of the vertices of the coincident polygon
+! 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_providing_polys'
+
+    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*2
+      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)
+      Ncoinpoly(ig)= Nsort
+      ccoinpoly(ig,:,:) = 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_providing_polys
 
   SUBROUTINE spacepercen(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals, xBBvals, yBBvals,      &
@@ -5132 +5241 @@
 
   SUBROUTINE grid_spacepercen(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals, xBBvals, yBBvals, &
-    dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Ngridsin, gridsin, areas, 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) which lay inside another 
-  !   series of grid-cells (A)
+  !   series of grid-cells (A) porviding coincident polygons
   ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x 
   !   and y axes. 
@@ -5207 +5316 @@
         END DO
 
-        CALL grid_spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals,       &
+        CALL grid_spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals,        &
           yBBvals, Ngridsin(ix,iy), poinsin, strict, pareas, percentages(ix,iy,:))
 
@@ -5225 +5334 @@
 
   END SUBROUTINE grid_spacepercen
+
+  SUBROUTINE grid_spacepercen_providing_polys(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals,   &
+    xBBvals, yBBvals, dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Nmaxvercoin,        &
+    Nvercoinpolys, vercoinpolys, 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) providing coincident polygons
+  ! 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
+    INTEGER, INTENT(in)                                  :: Nmaxvercoin
+    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,dxyB,Nmaxvercoin),                                                     &
+      INTENT(out)                                        :: Nvercoinpolys
+    REAL(r_k), DIMENSION(dxA,dyA,dxyB,Nmaxvercoin,2),                                                 &
+      INTENT(out)                                        :: vercoinpolys
+    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
+! Nvercoinpolys: number of vertices of the coincident polygon of each grid
+! coinpolys: of vertices of the coincident polygon of each grid
+! 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_providing_polys'
+
+    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,:,:))
+   
+        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_providing_polys(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, &
+          xBBvals, yBBvals, Ngridsin(ix,iy), poinsin, strict, Nmaxvercoin, Nvercoinpolys(ix,iy,:,:), &
+          vercoinpolys(ix,iy,:,:,:), 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
+    END DO
+
+    IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
+    IF (ALLOCATED(pareas)) DEALLOCATE(pareas)
+    IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)
+
+  END SUBROUTINE grid_spacepercen_providing_polys
 
   SUBROUTINE unique_matrixRK2D(dx, dy, dxy, matrix2D, Nunique, unique)
@@ -5768 +5981 @@
   END SUBROUTINE multi_index_mat4DRK
 
+  SUBROUTINE coincident_list_2Dcoords(NpointsA, pointsA, NpointsB, pointsB, Npoints, points, inpA,    &
+    inpB)
+  ! Subroutine to determine which 2D points of an A list are also found in a B list
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: NpointsA, NpointsB
+    INTEGER, DIMENSION(NpointsA,2), INTENT(in)           :: pointsA
+    INTEGER, DIMENSION(NpointsB,2), INTENT(in)           :: pointsB
+    INTEGER, INTENT(out)                                 :: Npoints
+    INTEGER, DIMENSION(NpointsA,2), INTENT(out)          :: points
+    INTEGER, DIMENSION(NpointsA), INTENT(out)            :: inpA, inpB
+
+    ! Local
+    INTEGER                                              :: iA, iB
+
+!!!!!!! Variables
+! NpointsA: Number of points of the list A
+! pointsA: points of the list A
+! NpointsB: Number of points of the list B
+! pointsB: points of the list B
+! Npoints: Number of coincident points
+! points: coincident points
+! inpA: coincident points list A
+! inpB: coincident points list B
+
+
+    fname = 'coincident_list_2Dcoords'
+
+    Npoints = 0
+    points = 0
+    inpA = 0
+    inpB = 0
+
+    DO iA = 1, NpointsA
+      DO iB = 1, NpointsB
+        IF ( (pointsA(iA,1) == pointsB(iB,1)) .AND. (pointsA(iA,2) == pointsB(iB,2)) ) THEN
+          Npoints = Npoints + 1
+          points(Npoints,:) = pointsA(iA,:)
+          inpA(Npoints) = iA
+          inpB(Npoints) = iB
+          EXIT
+        END IF
+
+      END DO
+    END DO
+
+  END SUBROUTINE coincident_list_2Dcoords
+
+  SUBROUTINE coincident_gridsin2D(dxA, dyA, dxyA, NpointsA, pointsA, dxB, dyB, dxyB, NpointsB,        &
+    pointsB, Npoints, points, inpointsA, inpointsB)
+  ! Subroutine to determine which lists of 2D gridsin points of an A list are also found in a B list
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dxA, dyA, dxyA
+    INTEGER, INTENT(in)                                  :: dxB, dyB, dxyB
+    INTEGER, DIMENSION(dxA, dyA), INTENT(in)             :: NpointsA
+    INTEGER, DIMENSION(dxB, dyB), INTENT(in)             :: NpointsB
+    INTEGER, DIMENSION(dxA, dyA, dxyA, 2), INTENT(in)    :: pointsA
+    INTEGER, DIMENSION(dxA, dyA, dxyB, 2), INTENT(in)    :: pointsB
+    INTEGER, DIMENSION(dxA, dyA, dxB, dyB), INTENT(out)  :: Npoints
+    INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA, 2),                                                  &
+      INTENT(out)                                        :: points
+    INTEGER, DIMENSION(dxA, dyA, dxyA), INTENT(out)      :: inpointsA
+    INTEGER, DIMENSION(dxB, dyB, dxyA), INTENT(out)      :: inpointsB
+
+    ! Local
+    INTEGER                                              :: ixA, iyA, ixB, iyB
+    INTEGER                                              :: NA, NB
+
+!!!!!!! Variables
+! dxA, dyA: 2D shape of the list A
+! NpointsA: 2D Number of points of the list A
+! pointsA: points of the list A
+! dxB, dyB: 2D shape of the list B
+! NpointsB: 2D Number of points of the list B
+! pointsB: points of the list B
+! Npoints: Number of coincident points
+! points: coincident points
+! inpointsA: coincident points list A
+! inpointsB: coincident points list B
+
+    fname = 'coincident_gridsin2D'
+
+    Npoints = 0
+    points = 0
+    inpointsA = 0
+    inpointsB = 0
+
+    DO ixA=1, dxA
+      DO iyA=1, dyA
+        NA = NpointsA(ixA,iyA)
+        DO ixB=1, dxB
+          DO iyB=1, dyB
+            NB = NpointsA(ixB,iyB)
+            CALL coincident_list_2Dcoords(NA, pointsA(ixA,iyA,1:NA,:), NB, pointsB(ixB,iyB,1:NB,:),   &
+              Npoints(ixA,iyA,ixB,iyB), points(ixA,iyA,ixB,iyB,:,:), inpointsA(ixA,iyA,:),            &
+              inpointsB(ixB,iyB,:))
+          END DO
+        END DO
+      END DO
+    END DO
+
+  END SUBROUTINE coincident_gridsin2D
+
 END MODULE module_scientific