Changeset 2554 in lmdz_wrf for trunk/tools/module_scientific.f90


Ignore:
Timestamp:
May 24, 2019, 8:57:12 PM (6 years ago)
Author:
lfita
Message:

Working version with reversing search from polygon to grid in `gridpoint_within'

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/module_scientific.f90

    r2553 r2554  
    49414941    REAL(r_k), DIMENSION(2,2)                            :: face, faceA
    49424942    LOGICAL, DIMENSION(dx,dy)                            :: within
     4943    REAL(r_k), DIMENSION(:,:), ALLOCATABLE               :: gridpoly
    49434944
    49444945!!!!!!! Variables
    … …  
    49904991          END DO
    49914992
    4992           ! Getting grid faces
    4993           DO iv=1, Nvertexmax-1
    4994             IF (.NOT.within(ix,iy)) THEN
    4995               IF (xBvals(ix,iy,iv) /= fillvalI) THEN
    4996                 face(1,:) = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /)
    4997                 face(2,:) = (/ xBvals(ix,iy,iv+1), yBvals(ix,iy,iv+1) /)
    4998                 DO ivA=1, NvertexA-1
    4999                   faceA(1,:) = (/ polygonA(ivA,1), polygonA(ivA,2) /)
    5000                   faceA(2,:) = (/polygonA(ivA+1,1), polygonA(ivA+1,2) /)
    5001                   CALL intersectfaces(face, faceA, cross, crosspoint)
    5002                   IF (cross /= 0) THEN
    5003                     Ngridsin = Ngridsin + 1
    5004                     ! Getting coordinates
    5005                     gridsin(Ngridsin,1) = ix
    5006                     gridsin(Ngridsin,2) = iy
    5007                     within(ix,iy) = .TRUE.
    5008                     CYCLE
    5009                   END IF
    5010                 END DO
    5011                 IF (.NOT.within(ix,iy)) THEN
    5012                   faceA(1,:) = (/ polygonA(NvertexA,1), polygonA(NvertexA,2) /)
    5013                   faceA(2,:) = (/polygonA(1,1), polygonA(1,2) /)
    5014                   CALL intersectfaces(face, faceA, cross, crosspoint)
    5015                   IF (cross /= 0) THEN
    5016                     Ngridsin = Ngridsin + 1
    5017                     ! Getting coordinates
    5018                     gridsin(Ngridsin,1) = ix
    5019                     gridsin(Ngridsin,2) = iy
    5020                     within(ix,iy) = .TRUE.
    5021                   END IF
    5022                 END IF
    5023               END IF
    5024             END IF
    5025           END DO
    5026           iv = Nvertexmax
    5027           IF (.NOT.within(ix,iy)) THEN
    5028             IF (xBvals(ix,iy,iv) /= fillvalI) THEN
    5029               face(1,:) = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /)
    5030               face(2,:) = (/ xBvals(ix,iy,1), yBvals(ix,iy,1) /)
    5031               DO ivA=1, NvertexA-1
    5032                 faceA(1,:) = (/ polygonA(ivA,1), polygonA(ivA,2) /)
    5033                 faceA(2,:) = (/polygonA(ivA+1,1), polygonA(ivA+1,2) /)
    5034                 CALL intersectfaces(face, faceA, cross, crosspoint)
    5035                 IF (cross /= 0) THEN
    5036                   Ngridsin = Ngridsin + 1
    5037                   ! Getting coordinates
    5038                   gridsin(Ngridsin,1) = ix
    5039                   gridsin(Ngridsin,2) = iy
    5040                   within(ix,iy) = .TRUE.
    5041                   CYCLE
    5042                 END IF
    5043               END DO
    5044               IF (.NOT.within(ix,iy)) THEN
    5045                 faceA(1,:) = (/ polygonA(NvertexA,1), polygonA(NvertexA,2) /)
    5046                 faceA(2,:) = (/polygonA(1,1), polygonA(1,2) /)
    5047                 CALL intersectfaces(face, faceA, cross, crosspoint)
    5048                 IF (cross /= 0) THEN
    5049                   Ngridsin = Ngridsin + 1
    5050                   ! Getting coordinates
    5051                   gridsin(Ngridsin,1) = ix
    5052                   gridsin(Ngridsin,2) = iy
    5053                   within(ix,iy) = .TRUE.
    5054                 END IF
    5055               END IF
    5056             END IF
    5057           END IF
     4993!          ! Getting grid faces
     4994!          DO iv=1, Nvertexmax-1
     4995!            IF (.NOT.within(ix,iy)) THEN
     4996!              IF (xBvals(ix,iy,iv) /= fillvalI) THEN
     4997!                face(1,:) = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /)
     4998!                face(2,:) = (/ xBvals(ix,iy,iv+1), yBvals(ix,iy,iv+1) /)
     4999!                DO ivA=1, NvertexA-1
     5000!                  faceA(1,:) = (/ polygonA(ivA,1), polygonA(ivA,2) /)
     5001!                  faceA(2,:) = (/polygonA(ivA+1,1), polygonA(ivA+1,2) /)
     5002!                  CALL intersectfaces(face, faceA, cross, crosspoint)
     5003!                  IF (cross /= 0) THEN
     5004!                    Ngridsin = Ngridsin + 1
     5005!                    ! Getting coordinates
     5006!                    gridsin(Ngridsin,1) = ix
     5007!                    gridsin(Ngridsin,2) = iy
     5008!                    within(ix,iy) = .TRUE.
     5009!                    CYCLE
     5010!                  END IF
     5011!                END DO
     5012!                IF (.NOT.within(ix,iy)) THEN
     5013!                  faceA(1,:) = (/ polygonA(NvertexA,1), polygonA(NvertexA,2) /)
     5014!                  faceA(2,:) = (/polygonA(1,1), polygonA(1,2) /)
     5015!                  CALL intersectfaces(face, faceA, cross, crosspoint)
     5016!                  IF (cross /= 0) THEN
     5017!                    Ngridsin = Ngridsin + 1
     5018!                    ! Getting coordinates
     5019!                    gridsin(Ngridsin,1) = ix
     5020!                    gridsin(Ngridsin,2) = iy
     5021!                    within(ix,iy) = .TRUE.
     5022!                  END IF
     5023!                END IF
     5024!              END IF
     5025!            END IF
     5026!          END DO
     5027!          iv = Nvertexmax
     5028!          IF (.NOT.within(ix,iy)) THEN
     5029!            IF (xBvals(ix,iy,iv) /= fillvalI) THEN
     5030!              face(1,:) = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /)
     5031!              face(2,:) = (/ xBvals(ix,iy,1), yBvals(ix,iy,1) /)
     5032!              DO ivA=1, NvertexA-1
     5033!                faceA(1,:) = (/ polygonA(ivA,1), polygonA(ivA,2) /)
     5034!                faceA(2,:) = (/polygonA(ivA+1,1), polygonA(ivA+1,2) /)
     5035!                CALL intersectfaces(face, faceA, cross, crosspoint)
     5036!                IF (cross /= 0) THEN
     5037!                  Ngridsin = Ngridsin + 1
     5038!                  ! Getting coordinates
     5039!                  gridsin(Ngridsin,1) = ix
     5040!                  gridsin(Ngridsin,2) = iy
     5041!                  within(ix,iy) = .TRUE.
     5042!                  CYCLE
     5043!                END IF
     5044!              END DO
     5045!              IF (.NOT.within(ix,iy)) THEN
     5046!                faceA(1,:) = (/ polygonA(NvertexA,1), polygonA(NvertexA,2) /)
     5047!                faceA(2,:) = (/polygonA(1,1), polygonA(1,2) /)
     5048!                CALL intersectfaces(face, faceA, cross, crosspoint)
     5049!                IF (cross /= 0) THEN
     5050!                  Ngridsin = Ngridsin + 1
     5051!                  ! Getting coordinates
     5052!                  gridsin(Ngridsin,1) = ix
     5053!                  gridsin(Ngridsin,2) = iy
     5054!                  within(ix,iy) = .TRUE.
     5055!               END IF
     5056!              END IF
     5057!            END IF
     5058!          END IF
    50585059
    50595060          ! Inverting process, look polygon within grid
    50605061          IF (.NOT.within(ix,iy)) THEN
    5061             DO iv = 1, NvertexA
    5062               vertex = polygonA(iv,:)
    5063               IF (point_inside(vertex, Nvertexmax, xBvals(ix,iy,:))) THEN
    5064                 Ngridsin = Ngridsin + 1
    5065                 ! Getting coordinates
    5066                 gridsin(Ngridsin,1) = ix
    5067                 gridsin(Ngridsin,2) = iy
    5068                 within(ix,iy) = .TRUE.
    5069                 CYCLE
     5062            IF (ALLOCATED(gridpoly)) DEALLOCATE(gridpoly)
     5063            ALLOCATE(gridpoly(Nvertexmax,2))
     5064            DO iv = 1, Nvertexmax
     5065              gridpoly(iv,1) = xBvals(ix,iy,iv)
     5066              gridpoly(iv,2) = yBvals(ix,iy,iv)
     5067            END DO
     5068            print *,ix,iy,'gridpoly :', gridpoly
     5069            DO ivA = 1, NvertexA
     5070              IF (.NOT.within(ix,iy)) THEN
     5071                 vertex = polygonA(ivA,:)
     5072                 print *,'  ', ivA, 'vertex :', vertex
     5073                 PRINT *,'  ',point_inside(vertex, Nvertexmax, gridpoly)
     5074                 IF (point_inside(vertex, Nvertexmax, gridpoly)) THEN
     5075                   Ngridsin = Ngridsin + 1
     5076                   ! Getting coordinates
     5077                   gridsin(Ngridsin,1) = ix
     5078                   gridsin(Ngridsin,2) = iy
     5079                   within(ix,iy) = .TRUE.
     5080                   CYCLE
     5081                 END IF
    50705082              END IF
    50715083            END DO
    … …  
    50745086      END DO
    50755087    END DO
     5088
     5089    IF (ALLOCATED(gridpoly)) DEALLOCATE(gridpoly)
    50765090
    50775091  END SUBROUTINE grid_within_polygon
Note: See TracChangeset for help on using the changeset viewer.