#include "swancpp.h" subroutine SwanBpntlist #if !defined SWAN_MODEL ! ! --|-----------------------------------------------------------|-- ! | Delft University of Technology | ! | Faculty of Civil Engineering and Geosciences | ! | Environmental Fluid Mechanics Section | ! | P.O. Box 5048, 2600 GA Delft, The Netherlands | ! | | ! | Programmer: Marcel Zijlema | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 1993-2017 Delft University of Technology ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! A copy of the GNU General Public License is available at ! http://www.gnu.org/copyleft/gpl.html#SEC3 ! or by writing to the Free Software Foundation, Inc., ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! ! Authors ! ! 40.80: Marcel Zijlema ! 40.92: Marcel Zijlema ! 41.14: Nico Booij ! 41.39: Clayton Hiles (Triton Consultants Ltd) ! ! Updates ! ! 40.80, April 2008: New subroutine ! 40.92, June 2008: changes with respect to boundary polygons ! 41.14, July 2010: boundary segments added as output curve ! 41.39, February 2011: modified to better handle occurrences of vertices with only 2 or 3 neighbours ! ! Purpose ! ! Makes list of boundary vertices in ascending order ! - counterclockwise in case of sea/mainland boundaries ! - clockwise in case of island boundaries ! ! Method ! ! The grid contains a number of boundary polygons ! They are by definition closed ! The first boundary polygon refers to sea/mainland boundary and the other polygons refers to island boundaries ! ! The vertices which define the sea/mainland boundary are inserted in the counterclockwise direction ! The vertices which define the island boundary are inserted in the clockwise direction ! ! Modules used ! use ocpcomm4 use SwanGriddata use SwanGridobjects use SwanCompdata use OUTP_DATA ! 41.14 ! implicit none ! ! Local variables ! integer :: icell ! cell index integer, parameter :: idebug=0 ! level of debug output: ! 0 = no output ! 1 = print extra output for debug purposes integer, save :: ient = 0 ! number of entries in this subroutine integer :: iface ! face index integer :: istat ! indicate status of allocation integer :: j ! loop counter integer :: k ! counter integer, dimension(1) :: kx ! location of minimum value in array of x-coordinates of boundary vertices integer, dimension(1) :: ky ! location of minimum value in array of y-coordinates of boundary vertices integer :: m ! loop counter integer :: maxnbp ! maximum number of boundary vertices in set of polygons integer :: nbptot ! total number of boundary vertices integer :: nptemp ! auxiliary integer to store number of points temporarily integer, dimension(3) :: v ! vertices in present cell integer :: v1 ! first vertex of present face integer :: v2 ! second vertex of present face integer :: vc ! considered vertex integer :: vcf ! first considered vertex of a boundary polygon integer :: vn ! next vertex with respect to considered vertex (counterclockwise) integer, dimension(3) :: fc ! cell face ID 41.39 integer :: icntfc ! counts number of faces without finding vc 41.39 integer :: MIP ! number of points on a output curve 41.14 integer :: vmk ! marker value of a boundary point 41.14 integer :: JJ ! counter of points on a curve 41.14 integer :: VM ! index of a boundary part 41.14 integer :: VMMAX ! highest value of VM 41.14 integer :: JBG ! index of a full boundary 41.14 integer :: IP, IPP ! point counter 41.14 integer :: IX ! vertex number 41.14 integer :: ISH ! shift number 41.14 ! integer, dimension(:), allocatable :: blistot ! list of all boundary vertices in ascending order integer, dimension(:), allocatable :: IARR1, IARR2 ! temporary array 41.14 ! real :: d1 ! distance of a point to origin real :: d2 ! distance of another point to origin real :: xp, yp ! coordinates of a boundary point 41.14 ! character(80) :: msgstr ! string to pass message character (len=8) :: PSNAME ! name of output curve 41.14 ! logical :: firstvert ! indicate whether considered vertex is first vertex of boundary polygon logical :: found ! indicates whether a new boundary part was found 41.14 ! type(celltype), dimension(:), pointer :: cell ! datastructure for cells with their attributes type(facetype), dimension(:), pointer :: face ! datastructure for faces with their attributes type(verttype), dimension(:), pointer :: vert ! datastructure for vertices with their attributes ! type(OPSDAT), pointer :: OPSTMP ! 41.14 type XYPT ! 41.14 real :: X, Y ! 41.14 type(XYPT), pointer :: NEXTXY end type XYPT type(XYPT), target :: FRST ! 41.14 type(XYPT), pointer :: CURR, TMP ! 41.14 ! ! Structure ! ! Description of the pseudo code ! ! Source text ! if (ltrace) call strace (ient,'SwanBpntlist') ! ! if list of boundary vertices is already filled, return ! if (allocated(blist)) return ! ! point to vertex, cell and face objects ! vert => gridobject%vert_grid cell => gridobject%cell_grid face => gridobject%face_grid ! vert(:)%atti(BINDX) = 0 vert(:)%atti(BPOL) = 0 nbpt = 0 ! ! determine total number of boundary vertices ! nbptot = count(mask=vert(:)%atti(VMARKER)==1) ! allocate(blistot(0:nbptot)) blistot = 0 ! ! determine first boundary vertex nearest to the origin ! kx = minloc(vert(:)%attr(VERTX), vert(:)%atti(VMARKER)==1) ky = minloc(vert(:)%attr(VERTY), vert(:)%atti(VMARKER)==1) ! if ( kx(1) == ky(1) ) then ! vc = kx(1) ! else ! d1 = sqrt((vert(kx(1))%attr(VERTX))**2+(vert(kx(1))%attr(VERTY))**2) d2 = sqrt((vert(ky(1))%attr(VERTX))**2+(vert(ky(1))%attr(VERTY))**2) ! if ( d1 < d2 ) then vc = kx(1) else vc = ky(1) endif ! endif ! ! store first boundary vertex ! vcf = vc blistot(1) = vc nbpol = 1 firstvert = .true. vert(vc)%atti(BPOL) = nbpol ! nptemp = 0 icntfc = 0 ! ! algorithm start to store next subsequent boundary vertices in ascending order ! k = 1 iface = 1 ! faceloop: do ! if ( face(iface)%atti(FMARKER) == 1 ) then ! if ( firstvert ) then ! icell = face(iface)%atti(FACEC1) ! ! identify the vertices and faces of the current cell ! v(1) = cell(icell)%atti(CELLV1) v(2) = cell(icell)%atti(CELLV2) v(3) = cell(icell)%atti(CELLV3) ! fc(1) = cell(icell)%face(1)%atti(FACEID) fc(2) = cell(icell)%face(2)%atti(FACEID) fc(3) = cell(icell)%face(3)%atti(FACEID) ! ! pick up next vertex (counterclockwise counting of vertices is assumed) ! vn = 0 do j = 1, cell(icell)%nov if ( v(j) == vc ) then vn = v(mod(j,cell(icell)%nov)+1) exit endif enddo ! if ( vn == 0 ) goto 10 if ( vert(vn)%atti(VMARKER) /= 1 ) goto 10 ! ! prevent algorithm from skipping sections of boundary by identifying where a bridge element is encountered 41.39 ! if ( j == 1 .and. face(fc(1))%atti(FMARKER) /= 1 ) then if ( ITEST >= 30 .or. idebug == 1 ) write(PRTEST,*) 'bridge element encountered on face 1' goto 10 endif if ( j == 2 .and. face(fc(2))%atti(FMARKER) /= 1 ) then if ( ITEST >= 30 .or. idebug == 1 ) write(PRTEST,*) 'bridge element encountered on face 2' goto 10 endif if ( j == 3 .and. face(fc(3))%atti(FMARKER) /= 1 ) then if ( ITEST >= 30 .or. idebug == 1 ) write(PRTEST,*) 'bridge element encountered on face 3' goto 10 endif ! if ( ITEST >= 30 .or. idebug == 1 ) write(PRTEST,*) 'next vertex found = ', vn ! firstvert = .false. ! endif ! ! we have found a correct boundary face and we continue to store subsequent boundary vertices ! v1 = face(iface)%atti(FACEV1) v2 = face(iface)%atti(FACEV2) ! if ( v1 == vc ) then ! if ( v2 == vcf .and. vcf /= blistot(k-1) ) vc = vcf if ( any( v2 == blistot ) ) goto 10 ! k = k + 1 blistot(k) = v2 vert(v2)%atti(BPOL) = nbpol vc = v2 ! icntfc = 0 ! reset counting of number of faces ! elseif ( v2 == vc ) then ! if ( v1 == vcf .and. vcf /= blistot(k-1) ) vc = vcf if ( any( v1 == blistot ) ) goto 10 ! k = k + 1 blistot(k) = v1 vert(v1)%atti(BPOL) = nbpol vc = v1 ! icntfc = 0 ! reset counting of number of faces ! elseif ( vc == vcf ) then ! end of considered boundary polygon is found ! ! prevent search algorithm from doubling back on itself 41.39 ! if ( vcf == blistot(k) ) then ! if ( ITEST >= 30 .or. idebug == 1 ) write(PRTEST,*) 'REPEAT!' k = k + 1 blistot(k) = vn vert(vn)%atti(BPOL) = nbpol vc = vn ! icntfc = 0 ! reset counting of number of faces ! goto 10 ! endif ! if ( any( v1 == blistot ) .and. any( v2 == blistot ) ) goto 10 ! ! store number of boundary vertices for present polygon ! nbpt(nbpol) = k - nptemp nptemp = k ! ! some diagnostics if ( ITEST >= 30 .or. idebug == 1 ) then write(PRTEST,*) 'END OF POLYGON' write(PRTEST,*) 'v1 = ', v1, ' v2 = ', v2 endif ! ! take first vertex of next boundary polygon ! vc = v1 vcf = vc k = k + 1 blistot(k) = vc nbpol = nbpol + 1 firstvert = .true. vert(vc)%atti(BPOL) = nbpol ! ! give error if more than 10000 boundary polygons are found ! if ( nbpol > 10000 ) call msgerr ( 2, ' More than 10000 boundary polygons are found in grid' ) ! icntfc = 0 ! reset counting of number of faces ! endif ! if ( k == nbptot ) exit faceloop ! endif ! 10 continue iface = iface + 1 if ( iface > nfaces ) iface = 1 ! icntfc = icntfc + 1 ! ! print diagnostics and stop program if count of number of faces suggests an endless loop 41.39 ! if ( icntfc > 4*nfaces ) then ! call msgerr ( 4, 'SwanBpntlist: list of boundary vertices could not be completed ' ) ! if ( ITEST >= 30 .or. idebug == 1 ) then write(PRTEST,*) 'error in SwanBpntlist: vertex vc not found' write(PRTEST,*) 'error in SwanBpntlist: vc = ', vc write(PRTEST,*) 'error in SwanBpntlist: vcf = ', vcf write(PRTEST,*) 'error in SwanBpntlist: k = ', k write(PRTEST,*) 'error in SwanBpntlist: nbptot = ', nbptot write(PRTEST,*) 'error in SwanBpntlist: blistot = ', blistot endif return ! endif ! enddo faceloop ! ! store number of boundary vertices for last polygon ! nbpt(nbpol) = nbptot - nptemp ! ! check if list contains boundary vertices only ! do j = 1, nbptot ! vc = blistot(j) if (vert(vc)%atti(VMARKER) /= 1) then write (msgstr, '(a,i4,a)') ' Vertex with index ',vc,' in boundary list is not a valid boundary point' call msgerr( 2, trim(msgstr) ) endif ! enddo ! ! determine maximum number of boundary vertices in set of polygons and allocate blist ! maxnbp = maxval(nbpt) ! istat = 0 if(.not.allocated(blist)) allocate (blist(maxnbp,nbpol), stat = istat) if ( istat /= 0 ) then call msgerr ( 4, 'Allocation problem in SwanBpntlist: array blist ' ) return endif blist = 0 ! ! fill blist in appropriate manner ! k = 0 ! do j = 1, nbpol ! do m = 1, nbpt(j) vc = blistot(k+m) blist(m,j) = vc vert(vc)%atti(BINDX) = m enddo ! k = k + nbpt(j) ! enddo ! deallocate(blistot) ! ! ! Add output curve corresponding to boundary 41.14 ! ALLOCATE(OPSTMP) OPSTMP%PSTYPE = 'C' MIP = nbpt(1) OPSTMP%MIP = MIP ALLOCATE(OPSTMP%XP(MIP)) ALLOCATE(OPSTMP%YP(MIP)) OPSTMP%PSNAME = 'BOUNDARY' do m = 1, MIP vc = blist(m,1) OPSTMP%XP(m) = vert(vc)%attr(VERTX) OPSTMP%YP(m) = vert(vc)%attr(VERTY) enddo IF (ITEST.GE.10) WRITE (PRTEST, 104) 'BOUNDARY', MIP 104 format (' Generated output curve ', A8, ' with ', I4, ' vertices.') ! ***** store number of points of the curve ***** NULLIFY(OPSTMP%NEXTOPS) IF ( .NOT.LOPS ) THEN FOPS = OPSTMP COPS => FOPS LOPS = .TRUE. ELSE COPS%NEXTOPS => OPSTMP COPS => OPSTMP END IF ! ! Determine highst value of VM ! VMMAX = 0 DO JBG = 1, nbpol DO IP = 1, nbpt(JBG) IX = blist(IP,JBG) VMMAX = MAX(VMMAX, vmark(IX)) ENDDO ENDDO !TEST write (prtest, *) 'test VMMAX ', VMMAX, nbpol ! ALLOCATE(IARR1(SUM(nbpt))) DO VM=1, VMMAX MIP = 0 JJ = 0 DO JBG = 1, nbpol ! ! first boundary polygon is assumed an outer one ! (sea/mainland boundary) and hence, content of blist ! is ordered in counterclockwise manner ! DO IP = 1, nbpt(JBG) IX = blist(IP,JBG) IF ( vmark(IX) == VM ) THEN MIP = MIP+1 IARR1(MIP) = IP if (JJ==0) then JJ=JBG elseif (JJ/=JBG) then !call msgerr (1, 'Side is part of 2 boundaries') endif ENDIF ENDDO ENDDO ! IF ( MIP/=0 ) THEN ! ALLOCATE(IARR2(MIP)) IARR2(1:MIP) = IARR1(1:MIP) ISH = 0 DO IPP = 2, MIP IF ( IARR2(IPP)/=IARR2(IPP-1)+1 ) THEN ISH = IPP-1 EXIT ENDIF ENDDO IARR2 = CSHIFT(IARR2,ISH) !TEST write (prtest, *) 'Shift ', ISH, MIP, IARR2(1) ! ALLOCATE(OPSTMP) OPSTMP%PSTYPE = 'C' OPSTMP%MIP = MIP ALLOCATE(OPSTMP%XP(MIP)) ALLOCATE(OPSTMP%YP(MIP)) write (PSNAME, 101) VM 101 format ('BOUND_',I2.2) OPSTMP%PSNAME = PSNAME DO IPP = 1, MIP IP = IARR2(IPP) IX = blist(IP,JJ) OPSTMP%XP(IPP) = vert(IX)%attr(VERTX) OPSTMP%YP(IPP) = vert(IX)%attr(VERTY) ENDDO DEALLOCATE(IARR2) IF (ITEST.GE.10) WRITE (PRTEST, 104) PSNAME, MIP NULLIFY(OPSTMP%NEXTOPS) IF ( .NOT.LOPS ) THEN FOPS = OPSTMP COPS => FOPS LOPS = .TRUE. ELSE COPS%NEXTOPS => OPSTMP COPS => OPSTMP END IF ! ENDIF ENDDO DEALLOCATE(IARR1) ! #endif end subroutine SwanBpntlist