#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