subroutine readpoly (iunit, np, nseg, npts, list, lat, lon, ier)
integer iunit, np, nseg, npts, ier
integer :: list(1:np)
real :: lat(1:np), lon(1:np)    

!***********************************************************
!   This subroutine reads in a file containing the line segments of the 
! bounding polygon for the interpolation region. The file format is the 
! standard ASCII output file from GRASS. The line segments are analyzed
! for connectivity and placed in order in the polygon lat, lon arrays.
! The LIST array is also initialized into a default point order.

! On input:
!       IUNIT = unit number of file opened for reading

! Input parameters are not altered by this function.

! On output:
!        NP = number of points in polygon array


!       IER = Error indicator:
!             IER = 0 if no errors were encountered.
!             IER = 1 file read error
!             IER = 2 polygon not closed
!***********************************************************

! Local variables
! idx - index variable
! nseg - number of segments
! npts - total number of points
! ns - number of points in a single segment
! tidx - index tracking total number of points read in
! sidx - index tracking total number of segments read in
! pidx - index tracking total number of points placed into polygon array
! seglat - array of latitude points for all segments
! seglon - array of longitude points for all segments

integer idx, ns, tidx, sidx, pidx, segleft, lastsegleft, ios, nidx
integer, allocatable :: isegbeg(:), isegend(:), isegnpt(:)
logical, allocatable :: segavail(:)
real, allocatable, target :: seglat(:), seglon(:)
character*256 line

! allocate index arrays for pointing to segment arrays
allocate ( isegbeg(1:nseg), isegend(1:nseg), isegnpt(1:nseg), segavail(1:nseg) )
! allocate arrays for points for all segments
allocate ( seglat(1:npts), seglon(1:npts) )

! read header lines until beginning of polygon definition found
ios = 0
tidx = 0
sidx = 0
read (iunit,'(a)',iostat=ios) line
do while(ios .eq. 0)
    if( (line(1:1) .eq. 'B') .and. (line(2:2) .eq. ' ') ) then
        read(line(2:),*) ns
        sidx = sidx + 1
        isegnpt(sidx) = ns
        ! read in segment points
        do idx=1,ns
            read (iunit,'(a)',iostat=ios) line
            if ( ios .ne. 0 ) then
                ier = 1
                return
            end if
            tidx = tidx + 1
            read( line, *) seglon(tidx), seglat(tidx)
            ! set indexes for beginning and end of segment
            if( idx .eq. 1 ) then
                isegbeg(sidx) = tidx
            else if( idx .eq. ns ) then
                isegend(sidx) = tidx
            end if
        end do
    end if
    read (iunit,'(a)',iostat=ios) line
end do

! initialize switch array showing which segments are available.
do idx = 1, nseg
    segavail(idx) = .true.
end do

! Move first segment into polygon array
segavail(1) = .false.
nidx = 0
do idx = 1, isegnpt(1)
    lat(idx) = seglat(idx)
    lon(idx) = seglon(idx)
    nidx = nidx + 1
end do

! find next segment to be attached to previous
segleft = nseg - 1
do while(segleft .gt. 0)
    lastsegleft = segleft
    do sidx=2,nseg
        if( segavail(sidx) ) then
            if( (lat(nidx) .eq. seglat(isegbeg(sidx))) .and. (lon(nidx) .eq. seglon(isegbeg(sidx))) ) then
                ! end of polygon segment matches beginning of this segment, add it
                do idx = isegbeg(sidx)+1, isegend(sidx)
                    nidx = nidx + 1
                    lat(nidx) = seglat(idx)
                    lon(nidx) = seglon(idx)
                end do
                segavail(sidx) = .false.
                segleft = segleft - 1
            else if( (lat(nidx) .eq. seglat(isegend(sidx))) .and. (lon(nidx) .eq. seglon(isegend(sidx))) ) then
                ! end of polygon segment matches end of this segment, add it backwards
                do idx = isegend(sidx)-1, isegbeg(sidx), -1
                    nidx = nidx + 1
                    lat(nidx) = seglat(idx)
                    lon(nidx) = seglon(idx)
                end do
                segavail(sidx) = .false.
                segleft = segleft - 1
            end if
        end if
    end do
    if( lastsegleft .eq. segleft ) then
        ! this means no matching segment was found even though there are segments left
        if( (lat(nidx) .eq. seglat(1)) .and. (lon(nidx) .eq. lon(1)) ) then
            ! polygon is closed, there are extra segments
            np = nidx
            exit
        end if
        ! polygon is not closed and no segment matches
        ier = 2
        return
    end if
end do

! create default list array
do idx = 1,np
    list(idx) = idx
end do

ier = 0

return
end

subroutine sizepoly (iunit, np, nseg, npts, ier)
integer iunit, np, nseg, npts, ier

!***********************************************************
!   This subroutine reads in a file containing the line segments of the 
! bounding polygon for the interpolation region. The file format is the 
! standard ASCII output file from GRASS. The numbers of line segments and 
! points in the file are counted and the number of points required to
! define the complete single line segment polygon returned.

! On input:
!       IUNIT = unit number of file opened for reading

! Input parameters are not altered by this function.

! On output:
!        NP = number of points in polygon array

!       IER = Error indicator:
!             IER = 0 if no errors were encountered.
!             IER = 1 file read error
!***********************************************************

! Local variables
! ios - returns status of file read
! idx - index variable
! nseg - number of segments
! npts - total number of points
! line - character array holding individual line contents

integer ios, idx
character*256 line

! Local parameters:

! read entire file, counting the number of segments available, and total points
ios = 0
nseg = 0
npts = 0
read (iunit,'(a)',iostat=ios) line
do while(ios .eq. 0)
    if( (line(1:1) .eq. 'B') .and. (line(2:2) .eq. ' ') ) then
        nseg = nseg + 1
        read(line(2:),*) idx
        npts = npts + idx
    end if
    read (iunit,'(a)',iostat=ios) line
end do

! this assumes that identical points will be dropped
np = npts - nseg + 1

rewind(iunit)

return
end