![(please configure the [header_logo] section in trac.ini)](/LMDZ_WRF/chrome/site/your_project_logo.png){kind=logo}
Changeset 2555 in lmdz_wrf
- Timestamp:
- May 24, 2019, 11:24:47 PM (6 years ago)
- File:
-
- 1 edited
-
trunk/tools/module_scientific.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/tools/module_scientific.f90
r2554 r2555 4991 4991 END DO 4992 4992 4993 !! Getting grid faces4994 !DO iv=1, Nvertexmax-14995 !IF (.NOT.within(ix,iy)) THEN4996 !IF (xBvals(ix,iy,iv) /= fillvalI) THEN4997 !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-15000 !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) THEN5004 !Ngridsin = Ngridsin + 15005 !! Getting coordinates5006 !gridsin(Ngridsin,1) = ix5007 !gridsin(Ngridsin,2) = iy5008 !within(ix,iy) = .TRUE.5009 !CYCLE5010 !END IF5011 !END DO5012 !IF (.NOT.within(ix,iy)) THEN5013 !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) THEN5017 !Ngridsin = Ngridsin + 15018 !! Getting coordinates5019 !gridsin(Ngridsin,1) = ix5020 !gridsin(Ngridsin,2) = iy5021 !within(ix,iy) = .TRUE.5022 !END IF5023 !END IF5024 !END IF5025 !END IF5026 !END DO5027 !iv = Nvertexmax5028 !IF (.NOT.within(ix,iy)) THEN5029 !IF (xBvals(ix,iy,iv) /= fillvalI) THEN5030 !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-15033 !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) THEN5037 !Ngridsin = Ngridsin + 15038 !! Getting coordinates5039 !gridsin(Ngridsin,1) = ix5040 !gridsin(Ngridsin,2) = iy5041 !within(ix,iy) = .TRUE.5042 !CYCLE5043 !END IF5044 !END DO5045 !IF (.NOT.within(ix,iy)) THEN5046 !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) THEN5050 !Ngridsin = Ngridsin + 15051 !! Getting coordinates5052 !gridsin(Ngridsin,1) = ix5053 !gridsin(Ngridsin,2) = iy5054 !within(ix,iy) = .TRUE.5055 !END IF5056 !END IF5057 !END IF5058 !END IF4993 ! 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 5059 5059 5060 5060 ! Inverting process, look polygon within grid … … 5066 5066 gridpoly(iv,2) = yBvals(ix,iy,iv) 5067 5067 END DO 5068 print *,ix,iy,'gridpoly :', gridpoly5069 5068 DO ivA = 1, NvertexA 5070 5069 IF (.NOT.within(ix,iy)) THEN 5071 5070 vertex = polygonA(ivA,:) 5072 print *,' ', ivA, 'vertex :', vertex5073 PRINT *,' ',point_inside(vertex, Nvertexmax, gridpoly)5074 5071 IF (point_inside(vertex, Nvertexmax, gridpoly)) THEN 5075 5072 Ngridsin = Ngridsin + 1
Note: See TracChangeset
for help on using the changeset viewer.
