!$Author$
!$Date$
!$Revision$
!$HeadURL$

!          INTERPOLATE:  Triangular interpolation on the unit sphere
!                        using STRIPACK
!                         05/28/09

!   Uses the software package STRIPACK to construct a Delaunay
! triangulation and Voronoi diagram of a set of nodes on the surface
! of the unit sphere, and then use that to interpolate a value from
! those nodes to a point.

! standard input consists of a list of any number of lat,lon pairs

! The nodes are read in from the windgen histogram database index file
! specified on the command line. These values are then used to derive:
!    (RLAT(I),RLON(I), I = 1,nval) = Nodal coordinates in
!                 degrees latitude (-90 to 90) and degrees
!                 longitude (-180 to 180)

! The nodes and points are screened for presence within a bounding polygon
! which is input from a file designated on the command line.
!
! The format of this input bounding polygon is a GRASS standard format
! ASCII output file defining points of type B (boundary), which is then
! transformed into a Counter Clock Wise ordered list of point pairs,
! LONGITUDE, then LATITUDE GIS convention, with the last point linking
! to the first point.

! Example input bounding polygon (remove ! from beginning of line)
! where: B indicates the beginning of a line segment
!        the number after B indicates the number of points in the segment
!        lines with longitude, latitude numbers for the points start with a space
!        The first and last points in each segment match to other segments
!        forming a closed polygon.
! note: this program does not use the header lines at all

!ORGANIZATION: 
!DIGIT DATE:   
!DIGIT NAME:   fredfox
!MAP NAME:     
!MAP DATE:     Tue May 19 11:00:04 2009
!MAP SCALE:    1
!OTHER INFO:   
!ZONE:         0
!MAP THRESH:   0.000000
!VERTI:
!B  7
! -104.23332673 32.96272671 
! -104.26025644 33.11496538 
! -104.27960716 33.24932663 
! -104.31144531 33.39107902 
! -104.32906042 33.47254686 
! -104.35506583 33.55777434 
! -104.41459214 33.65805705 
!B  20
! -104.00220319 40.54063452 
! -104.09150432 40.7046569  
! -104.11869136 40.87716633 
! ...

!
! Only points inside the polygon are interpolated and output. The polygon
! is also considered the boundary for the interpolation region, with
! points mirrored from inside to outside or just to the boundary to handle
! triangulation or nearest neighbor interpolation near the boundary.

! This program must be linked to STRIPACK.

! for command line processing
integer iargc
character*256 polyfilename, histfilename, argv, outputfilename
character*256 line
character*25 fmt_str_lat, fmt_str_lon
character*25 fmt_str_lat1, fmt_str_lon1
character*25 fmt_str_in1, fmt_str_in2
character*25 fmt_str_nb
integer :: inplon_len = 0, inplat_len = 0
integer :: inplon_dlen = 0, inplat_dlen = 0
real :: inplon, inplat

! for standard in processing
integer :: idx, idxm1, idx1, jdx, jdx1, kdx   ! loop counters
integer :: ntr  ! the count of numbers to be read
integer :: ios   ! iostat return value
!real :: temp, in(2), point(3), px(1), py(1), pz(1)
real :: temp, in(2), point(3), savpoint(3), bndpoint(3), refpoint(3)
integer :: in_len(2), in_dlen(2)
real :: mx(4), my(4), mz(4), px(4), py(4), pz(4), intx(1), inty(1), intz(1)
integer :: plist(4)
real, allocatable :: trx(:), try(:), trz(:)
real, allocatable :: t1(:), t2(:), t3(:), t4(:), t5(:), t6(:), trlat(:), trlon(:)
integer, allocatable :: s1(:)
real :: latitude, longitude
real :: latmin, latmax, lonmin, lonmax
real :: plat(1:4), plon(1:4), mlat(1:4), mlon(1:4)

integer :: nval  ! the count of wind histrogram nodes to be read
integer :: lat_deg, lat_min, lon_deg, lon_min, elev, num, stnum, idxnum
character :: lat_NS, lon_EW, stat1, stat2, stat3
character*8 :: date1, date2
integer, allocatable :: stanum(:)

! bounding polygon variables
integer :: nline   ! line of file where polygon data starts
integer :: nppoly  ! count of points in bounding polygon
integer, allocatable :: polylist(:)
real, allocatable :: polylat(:), polylon(:), polyx(:), polyy(:), polyz(:)
real, allocatable :: biglat(:), biglon(:), bigx(:), bigy(:), bigz(:)
real, allocatable :: tbiglat(:), tbiglon(:), tbigx(:), tbigy(:), tbigz(:)
real, allocatable :: litlat(:), litlon(:), litx(:), lity(:), litz(:)
real, allocatable :: tlitlat(:), tlitlon(:), tlitx(:), tlity(:), tlitz(:)
real :: maxdist, incdist, deltadist  ! distance polygons are expanded, shrunk to get mirrored points
real :: bigdist, litdist    ! distance the little polygon has been moved.
real :: rlitdist  ! remaining distance for little polygon to be shrunk to get mirrored points
integer ninc   ! number of steps to expand/shrink polygons

! box and reflection tracking variables, m - point to be mirrored, r - reflected point
integer mval, tval, eval
real, allocatable :: rptx(:), rpty(:), rptz(:), rptlat(:), rptlon(:)
logical, allocatable :: mir_flg(:)

! triangulation variables
integer ier, lwk, lnew, nrow, nt6, ntmx, vrmx
real    sc
parameter (nrow=9)

! Array storage for the triangulation, work space, and nodal coordinates.
integer, allocatable :: list(:), lptr(:), lend(:), iwk(:)
real, allocatable :: ds(:), rlat(:), rlon(:), x(:), y(:), z(:)

! Array storage for the triangle list.
integer, allocatable :: ltri(:,:)

! return values from trfind
real :: b1, b2, b3, nb1, nb2, nb3
integer :: i1, i2, i3
integer :: i

! return values for bnodes
integer, allocatable :: bounds(:)
integer :: nbnd, narcs, ntri

! return values for inside
logical tin, tin1, tin2, sel_point

! inverse distance point count
integer :: invdist_pts

! command line input flags
logical cligen, debug, voronoi, windgen, invdist
! command line input flag indicating value to follow argument
logical hist_f, out_f, nopoly_f, poly_f, lon_f, lat_f, invdist_f, bound_f

! function called
logical inside
logical check_cross
logical arc_longer

! working variables for CRLIST and VROUT (VORONOI diagram output)
integer, allocatable :: listc(:), lbtri(:,:)
real, allocatable :: xc(:), yc(:), zc(:), rc(:), vclon(:), vclat(:)

! variables for inverse distance lists
integer, allocatable :: seln(:)
real, allocatable :: dist(:), idw(:), nidw(:)

! set defaults for flags
cligen = .FALSE.
debug = .FALSE.
voronoi = .FALSE.
windgen = .TRUE. ! windgen index is the default nodal file format
hist_f = .FALSE.
nopoly_f = .FALSE.
out_f = .FALSE.
poly_f = .FALSE.
bound_f = .TRUE. ! mirror points close to the bounding polygon only to polygon boundary
lon_f = .FALSE.
lat_f = .FALSE.
invdist_f = .FALSE.
invdist = .FALSE.
invdist_pts = 0

! set default input file names
histfilename = './wind_gen_his_upper_US.idx'
polyfilename = './boundary.pol'
outputfilename = './interpolate.txt'
inplat = 1000.0; !Initial values used to check if correctly initialized (or read in)
inplon = 1000.0;

! set coordinate transformation constant
sc = atan(1.)/45.

! process command line
if (iargc() .gt. 0) then
  !write(*,*) 'Argc ', iargc()
  do idx = 1, iargc()
    call getarg(idx,argv)
    !write(*,*) 'idx ', idx
    !write(*,*) 'Argv ', argv
    ! if flag set, then read second argument to command line switch
    if( hist_f ) then			! switch indicates histogram index file name is in next argument
        histfilename = argv		!wind histogram index file name
        hist_f = .FALSE.
    else if( poly_f ) then		! switch indicates bounding polygon file name is in next argument
        polyfilename = argv		!bounding polygon file name
        poly_f = .FALSE.
    else if( out_f ) then		! switch indicates output file name is in next argument
        outputfilename = argv	!output file name
        open(6, file=outputfilename, iostat=ios)
        if( ios .ne. 0 ) then
            write(*,*) 'ERROR, unable to open output file named: ', outputfilename
            stop
        end if
        out_f = .FALSE.
    else if( lat_f ) then		! switch indicates lattitude is in next argument
        inplat_len = len_trim(trim(argv))
         !print *, 'index_of_period:', index(trim(argv),".")
         !print *, 'substring_with_period:', argv(index(trim(argv),"."):inplat_len)
         !print *, 'substring_without_period:', argv(index(trim(argv),".")+1:inplat_len)
        inplat_dlen = len_trim(argv(index(trim(argv),".")+1:inplat_len))
        read(argv, *) inplat
        !write(*, FMT='(A)', ADVANCE='NO') 'Lat: '
        !write(fmt_str_lat, '("(f", i0, ".", i0")")') inplat_len, inplat_dlen
        !write(*, fmt_str_lat, ADVANCE='NO') inplat
        !write(*, FMT='(A,i2,A,i2)', ADVANCE='YES') ' len=', inplat_len, ' dec len=', inplat_dlen
        lat_f = .FALSE.
    else if( lon_f ) then		! switch indicates lattitude is in next argument
        inplon_len = len_trim(trim(argv))
        inplon_dlen = len_trim(argv(index(trim(argv),".")+1:inplon_len))
        read(argv, *) inplon
        !write(*, FMT='(A)', ADVANCE='NO') 'Lon: '
        !write(fmt_str_lon, '("(f", i0, ".", i0")")') inplon_len, inplon_dlen
        !write(*, fmt_str_lon, ADVANCE='NO') inplon
        !write(*, FMT='(A,i2,A,i2)', ADVANCE='YES') ' len=', inplon_len, ' dec len=', inplon_dlen
        lon_f = .FALSE.
    else if( invdist_f ) then		! switch indicates inverse distance point count is in next argument
        read(argv, *) invdist_pts
        !write(*,*) 'invdist_pts ', invdist_pts
        if( invdist_pts .lt. 1 ) then
            write(*,*) 'Specified point count for inverse distance weighting was ', invdist_pts
            write(*,*) 'At least one point is required.'
            stop
        end if
        invdist_f = .FALSE.
    else

        !command line help prompt
        if ((argv(2:2).eq.'?').or.(argv(2:2).eq.'h')) then
          write(*,*) 'Valid command line options:'
          write(*,*) '-? or -h Display this help screen'
          write(*,*) '-b Boundary points reflected only to boundary, not beyond'
          write(*,*) '-c index file is in cligen format'
          write(*,*) '-d create debugging output files'
          write(*,*) '-f [filename] windgen or cligen index file name'
          write(*,*) '-i [#] use inverse distance weighting with [#] points'
          write(*,*) '-lat [#] Lattitude (+ north)'
          write(*,*) '-lon [#] Longitude (+ east)'
          write(*,*) '-m Boundary points mirrored to outside bounding polygon'
          write(*,*) '-n Bounding polygon NOT used'
          write(*,*) '-o [filename] output file name'
          write(*,*) '-p [filename] Bounding polygon file name'
          write(*,*) '-t use Delauney triangulation weighting'
          write(*,*) '-v [filename] create voronoi output file'
          write(*,*) '-w index file is in windgen format'
          stop
        else if (argv(2:2) == 'b') then
            bound_f = .TRUE.
        else if (argv(2:2) == 'c') then
            cligen = .TRUE.
            windgen = .FALSE.
            write(*,*) 'Entered command arg c'
        else if (argv(2:2) == 'd') then
            debug = .TRUE.
        else if (argv(2:2) == 'f') then
            hist_f = .TRUE.
        else if (argv(2:2) == 'i') then
            invdist = .TRUE.
            invdist_f = .TRUE.
        else if (argv(2:4) == 'lat') then
            lat_f = .TRUE.
        else if (argv(2:4) == 'lon') then
            lon_f = .TRUE.
        else if (argv(2:2) == 'm') then
            bound_f = .FALSE.
        else if (argv(2:2) == 'n') then
            nopoly_f = .TRUE.
        else if (argv(2:2) == 'o') then
            out_f = .TRUE.
        else if (argv(2:2) == 'p') then
            poly_f = .TRUE.
        else if (argv(2:2) == 't') then
            invdist = .FALSE.
        else if (argv(2:2) == 'v') then
            voronoi = .TRUE.
        else if (argv(2:2) == 'w') then
            windgen = .TRUE.
        end if
    end if
  end do
endif
write(*,*) '# Using ', histfilename(1:len_trim(histfilename)), ' for known point coordinates'

if( nopoly_f ) then
    write(*,*) '# Not Using a bounding polygon'
else
    write(*,*) '# Using ', polyfilename(1:len_trim(polyfilename)), ' for bounding polygon'

    ! read in bounding polygon
    ! count number of data values in file
    write(*,*) '# Opening ', polyfilename(1:len_trim(polyfilename))
    open(2,FILE=polyfilename(1:len_trim(polyfilename)),STATUS='OLD',iostat=ios)
    if( ios .ne. 0 ) then
        write(*,*) 'ERROR, unable to open bounding polygon file named: ', polyfilename(1:len_trim(polyfilename))
        stop
    end if
    ! read number of polygon values and input data counts from file, then rewind
    call sizepoly (2, nppoly, idx, idx1, ier)
    ! allocate arrays for polygon
    allocate ( polylist(1:nppoly), polylat(1:nppoly), polylon(1:nppoly) )
    ! read in and assemble polygon points
    call readpoly ( 2, nppoly, idx, idx1, polylist, polylat, polylon, ier )
    if( ier .eq. 2 ) then
        write(*,*) "# ERROR: polygon does not close properly"
    end if
    close(2)
    ! debug, output polygon read for gnuplot
    if( debug ) then
        open(3,FILE='windgen_int_poly_as_read.pol',STATUS='REPLACE')
        write(3,1000) 'B ', nppoly
        do idx = 1, nppoly
            write(3,1020) polylon(polylist(idx)), polylat(polylist(idx))
        end do
        close(3)
    end if

    ! allocate arrays for polygon in cartesian coordinates
    allocate ( polyx(1:nppoly), polyy(1:nppoly), polyz(1:nppoly) )

    ! set up arrays for bounding polygon checking
    ! Set X and Y to the values of RLON and RLAT, respectively,
    ! in radians.  (RLON and RLAT are saved for output)
    do idx = 1, nppoly
        polyx(idx) = sc * polylon(idx)
        polyy(idx) = sc * polylat(idx)
    end do

    ! Transform spherical coordinates X and Y to Cartesian
    ! coordinates (X,Y,Z) on the unit sphere (X**2 + Y**2 + Z**2 = 1).
    call trans(nppoly, polyy, polyx, polyx, polyy, polyz)

    ! check polygon for single hemisphere coverage (not done at this point)

    ! check clockwise or counter-clockwise order of points.
    ! select point know to be outside of polygon
    ! (for our North American Polygons we use the South Pacific over water)
    px(1) = sc * (-115.0) ! longitutde
    py(1) = sc * (-15.0) ! latitude
    call trans(1, py, px, px, py, pz)
    refpoint(1) = px(1)
    refpoint(2) = py(1)
    refpoint(3) = pz(1)
    if( inside(refpoint, nppoly, polyx, polyy, polyz, nppoly, polylist, ier) ) then
        ! Other hemisphere should be outside, reverse point order
        do idx = 1, nppoly
            polylist(idx) = nppoly - idx + 1
        end do
        write(*,*)'# Polygon List order reversed'
        ! debug, output polygon read for gnuplot
        if( debug ) then
            open(3,FILE='windgen_int_poly_as_reordered.pol',STATUS='REPLACE')
            write(3,1000) 'B ', nppoly
            do idx = 1, nppoly
                write(3,1020) polylon(polylist(idx)), polylat(polylist(idx))
            end do
            close(3)
        end if
    else
        write(*,*)'# Polygon List order unchanged'
    end if
end if

! read points that are to be interpolated
! read stdin, check if in bounding polygon, dynamically expanding array as needed
allocate (t1(1:1), t2(1:1), t3(1:1), t4(1:1), t5(1:1))
t1(1) = 0.0
t2(1) = 0.0
t3(1) = 0.0
t4(1) = 0.0
t5(1) = 0.0

! read lat/lon pairs to be interpolated
ntr = 0
if( inplon.gt.999.0 ) then
    if( inplat.gt.999.0 ) then
        !Not sure what is going on here
        !Will leave alone for now - LEW
        read(*,*, iostat=ios) in(1), in(2)

!Was attempting to read in data similar to inplat/inplon, but can't get program to reach the code here - LEw
!        read(*,*, iostat=ios) line
!      write(*,*) 'line: ', line
!        read(line,*, iostat=ios) in(1), in(2)
!      write(*,*) 'in(1):', in(1), ' in(2):', in(2)
!Not 100% working code here yet - LEW
!        in_len(1) = len_trim(trim(in(1)))
!        in_dlen(1) = len_trim(in(1)(index(trim(in(1)),".")+1:in_len(1)))
!        write(*, FMT='(A)', ADVANCE='NO') 'Lat: '
!        write(fmt_str_in1, '("(f", i0, ".", i0")")') in_len(1), in_dlen(1)
!        write(*, fmt_str_in1, ADVANCE='NO') in(1)
!        write(*, FMT='(A,i2,A,i2)', ADVANCE='YES') ' len=', in_len(1), ' dec len=', in_dlen(1)
!        in_len(2) = len_trim(trim(in(2)))
!        in_dlen(2) = len_trim(in(2)(index(trim(in(2)),".")+1:in_len(2)))
!        write(*, FMT='(A)', ADVANCE='NO') 'Lon: '
!        write(fmt_str_in2, '("(f", i0, ".", i0")")') in_len(2), in_dlen(2)
!        write(*, fmt_str_in2, ADVANCE='NO') in(2)
!        write(*, FMT='(A,i2,A,i2)', ADVANCE='YES') ' len=', in_len(2), ' dec len=', in_dlen(2)

    else
        write(*,*) 'ERROR: Longitude needed as command line argument (-lon)'
        stop
    end if
else
    if( inplat.gt.999.0 ) then
        write(*,*) 'ERROR: Latitude needed as command line argument (-lat)'
        stop
    end if
    in(1) = inplat
    in(2) = inplon
    ios = 0

    !Not sure why the above is done, but added the following so we can print out the lat/lon data to desired number of digits, etc. - LEW
    in_len(1) = inplat_len
    in_len(2) = inplon_len
    in_dlen(1) = inplat_dlen
    in_dlen(2) = inplon_dlen
    !write(*,*) '# lat/lon ', in(1), in(2)

    write(*,'(A)',ADVANCE='NO') ' # lat/lon '
    write(fmt_str_lat1, '("(f", i0, ".", i0")")') inplat_len, inplat_dlen
    write(fmt_str_lon1, '("(f", i0, ".", i0")")') inplon_len+1, inplon_dlen
    write(*, fmt_str_lat1, ADVANCE='NO') inplat
    write(*, fmt_str_lon1, ADVANCE='YES') inplon

endif
do while(ios .eq. 0)
    px(1) = sc * in(2)
    py(1) = sc * in(1)
    call trans(1, py, px, px, py, pz)
    point(1) = px(1)
    point(2) = py(1)
    point(3) = pz(1)
    ! check point location (inside or outside of polygon) or (polygon not used)
    if( nopoly_f ) then
        sel_point = .true.
    else
        sel_point = inside(point, nppoly, polyx, polyy, polyz, nppoly, polylist, ier)
    end if
    if( sel_point ) then
        ! point is inside or polygon not used so take all points
       if( ntr .gt. 0 ) then
            allocate (t1(1:ntr), t2(1:ntr), t3(1:ntr), t4(1:ntr), t5(1:ntr) )
            ! whole array assignment
            t1 = trlat
            t2 = trlon
            t3 = trx
            t4 = try
            t5 = trz
            deallocate (trlat, trlon, trx, try, trz)
        end if
        ! it is inside, add to list to be processed
        ntr = ntr + 1
        allocate (trlat(1:ntr), trlon(1:ntr), trx(1:ntr), try(1:ntr), trz(1:ntr))
        ! whole array assignment
        do idx = 1, ntr-1
            trlat(idx) = t1(idx)
            trlon(idx) = t2(idx)
            trx(idx) = t3(idx)
            try(idx) = t4(idx)
            trz(idx) = t5(idx)
        end do
        deallocate (t1, t2, t3, t4, t5)
        trlat(ntr) = in(1)
        trlon(ntr) = in(2)
        trx(ntr) = px(1)
        try(ntr) = py(1)
        trz(ntr) = pz(1)
    end if
    if (inplon.gt.999.0) then
        !Again, not sure what is going on here - LEW
        read(*,*, iostat=ios) in(1), in(2)
    else
        ios = 1
    endif
end do

if( ntr .eq. 0 ) then
    write(*,*) '# No points to be interpolated'
    stop
end if

! process input nodes from point index file (windgen, cligen)
! count number of data values in file
if( windgen ) then
    write(*,*) '#Using Windgen format index file'
else if( cligen ) then
    write(*,*) '#Using Cligen format index file'
end if

write(*,*) '#Opening ', histfilename(1:len_trim(histfilename))
open(2,FILE=histfilename(1:len_trim(histfilename)),STATUS='OLD',iostat=ios)
if( ios .ne. 0 ) then
    write(*,*) 'ERROR, unable to open input nodes file named: ', histfilename(1:len_trim(histfilename))
    stop
end if

! read nodal points, check if in bounding polygon, dynamically expanding array as needed
allocate (t1(1:1), t2(1:1), t3(1:1), t4(1:1), t5(1:1), s1(1:1))
t1(1) = 0.0
t2(1) = 0.0
t3(1) = 0.0
t4(1) = 0.0
t5(1) = 0.0
s1(1) = 0

nval = 0
if( windgen ) then
    read (2,*,iostat=ios) lat_deg, lat_min, lat_NS, lon_deg, lon_min, lon_EW, elev, date1, date2, stat1, stat2, stat3, num
else if( cligen ) then
    read (2,*,iostat=ios) stnum, num, idxnum, lat_deg, lat_min, lat_NS, lon_deg, lon_min, lon_EW, elev
end if
do while(ios .eq. 0)
    latitude = lat_deg + lat_min/60.0
    if (lat_NS=='S') latitude = -1*latitude
    longitude = lon_deg + lon_min/60.0
    if (lon_EW=='W') longitude = -1*longitude
    px(1) = sc * longitude
    py(1) = sc * latitude
    call trans(1, py, px, px, py, pz)
    point(1) = px(1)
    point(2) = py(1)
    point(3) = pz(1)
    ! check if point is (inside bounding polygon) or (no polygon is used)
    if( nopoly_f ) then
        sel_point = .true.
    else
        sel_point = inside(point, nppoly, polyx, polyy, polyz, nppoly, polylist, ier)
    end if
    if( sel_point ) then
        ! point is inside or polygon not used so take all points
        if( nval .gt. 0 ) then
            allocate (t1(1:nval), t2(1:nval), t3(1:nval), t4(1:nval), t5(1:nval), s1(1:nval) )
            ! whole array assignment
            t1 = rlat
            t2 = rlon
            t3 = x
            t4 = y
            t5 = z
            s1 = stanum
            deallocate (rlat, rlon, stanum, x, y, z)
        end if
        ! it is inside, add to list to be processed
        nval = nval + 1
        allocate (rlat(1:nval), rlon(1:nval),stanum(1:nval), x(1:nval), y(1:nval), z(1:nval))
        ! whole array assignment
        do idx = 1, nval-1
            rlat(idx) = t1(idx)
            rlon(idx) = t2(idx)
            x(idx) = t3(idx)
            y(idx) = t4(idx)
            z(idx) = t5(idx)
            stanum(idx) = s1(idx)
        end do
        deallocate (t1, t2, t3, t4, t5, s1)
        rlat(nval) = latitude
        rlon(nval) = longitude
        stanum(nval) = num
        x(nval) = px(1)
        y(nval) = py(1)
        z(nval) = pz(1)
    end if
    if( windgen ) then
        read (2,*,iostat=ios) lat_deg, lat_min, lat_NS, lon_deg, lon_min, lon_EW, elev, date1, date2, stat1, stat2, stat3, num
    else if( cligen ) then
        read (2,*,iostat=ios) stnum, num, idxnum, lat_deg, lat_min, lat_NS, lon_deg, lon_min, lon_EW, elev
    end if
end do

if (nval .lt. 3) then
    write(*,*) "#ERROR: need at least three points for triangulation"
    stop
end if
close(2)

! debug, output points within bounding polygon for gnuplot
if( debug ) then
    open(3,FILE='windgen_int_inbound.pts',STATUS='REPLACE')
    do idx = 1, nval
        write(3,*) rlon(idx), rlat(idx)
    end do
    close(3)
end if

if( nopoly_f .or. invdist ) then
    ! no bounding polygon, so triangulate with point set as given
    ! or
    ! inverse distance weighting, polygon point mirroring not wanted
    eval = nval
else
    ! create interior polygon for mirror station selection
    allocate (litlat(1:nppoly), litlon(1:nppoly), litx(1:nppoly), lity(1:nppoly), litz(1:nppoly))
    allocate (tlitlat(1:nppoly), tlitlon(1:nppoly), tlitx(1:nppoly), tlity(1:nppoly), tlitz(1:nppoly))
    ! create exterior polygon for mirror station placement
    allocate (biglat(1:nppoly), biglon(1:nppoly), bigx(1:nppoly), bigy(1:nppoly), bigz(1:nppoly))
    allocate (tbiglat(1:nppoly), tbiglon(1:nppoly), tbigx(1:nppoly), tbigy(1:nppoly), tbigz(1:nppoly))
    !set mir_flg to .TRUE. if station is mirrored, .FALSE. otherwise
    allocate ( mir_flg(1:nval) )
    ! allocate storage for reflected points.
    allocate( rptx(1:nval), rpty(1:nval), rptz(1:nval), rptlat(1:nval), rptlon(1:nval) )
    mval = 0 ! zero out conter for number of mirrored points (reflected points)
    ! set up point index for four point array
    do idx = 1,4
        plist(idx) = idx
    end do
    ! copy the polygon into a temporary polygon for stepping to smaller
    do idx = 1,nppoly
        tlitx(idx) = polyx(idx)
        tlity(idx) = polyy(idx)
        tlitz(idx) = polyz(idx)
    end do
    ! copy the polygon into a temporary polygon for stepping to larger
    do idx = 1,nppoly
        tbigx(idx) = polyx(idx)
        tbigy(idx) = polyy(idx)
        tbigz(idx) = polyz(idx)
    end do

    ! initialize flags for points to be mirrored
    do idx = 1,nval
        mir_flg(idx) = .false.
    end do

    if( debug ) then
        ! debug, open file for boxes
        write(*,*) "# Open windgen_int_boxes.pol"
        open(3,FILE='windgen_int_boxes.pol',STATUS='REPLACE')
    end if

    ! set distance used to mirrored points
    maxdist = (200.0/6378.0)
    ! move polygons in smaller increments, to get better placement of near boundary points.
    ninc = 10
    incdist = maxdist / ninc
    do while( ninc .gt. 0 )
        ninc = ninc - 1
        ! set starting distance the two polygons have been moved
        litdist = 0.0
        bigdist = 0.0
        rlitdist = -incdist
        do while( rlitdist .lt. 0.0 )
            ! makes a smaller polygon, returns negative ier if terminates before full distance reached
            ! this behavior built into polysizing to support early point mirroring on peninsulas that
            ! disappear in the smaller polygon
            call polysizing (rlitdist, nppoly, polylist, tlitx, tlity, tlitz, litx, lity, litz, ier)
            ! now move outer polygon the same distance
            deltadist = incdist + litdist + rlitdist ! made a positive number
            !write(*,*) "# incdist, litdist, deltadist, rlitdist: ", incdist, litdist, deltadist, rlitdist
            ! update bookmarks before deltadist gets zeroed
            bigdist = bigdist + deltadist
            litdist = -incdist - rlitdist
            ! move big polygon
            call polysizing (deltadist, nppoly, polylist, tbigx, tbigy, tbigz, bigx, bigy, bigz, ier)
            ! copy the smaller polygon into larger polygon so ready for next step
            do idx = 1,nppoly
                tlitx(idx) = litx(idx)
                tlity(idx) = lity(idx)
                tlitz(idx) = litz(idx)
            end do
            ! copy the larger polygon into smaller polygon so ready for next step
            do idx = 1,nppoly
                tbigx(idx) = bigx(idx)
                tbigy(idx) = bigy(idx)
                tbigz(idx) = bigz(idx)
            end do
            ! mirror points that can be mirrored at this step and mark as mirrored.
            do idx = 1,nval
                if( .not. mir_flg(idx) ) then
                    point(1) = x(idx)
                    point(2) = y(idx)
                    point(3) = z(idx)
                    tin = inside(point, nppoly, litx, lity, litz, nppoly, polylist, ier)
                    if( ier .gt. 0) then
                        write(*,*) '#',ier,'error in inside routine. idx=',idx, 'idx1=', idx1
                        cycle
                    else if( tin ) then
                        ! point is inside small polygon so it is rejected
                        mir_flg(idx) = .false.
                        rptx(idx) = point(1)
                        rpty(idx) = point(2)
                        rptz(idx) = point(3)
                    else
                        ! point is outside small polygon so it is accepted (it is already known to be inside boundary polygon)
                        mval = mval + 1
                        mir_flg(idx) = .true.
                        ! find box containing point to be mirrored
                        do jdx = 1,nppoly
                            ! set box points (counter clock wise box order), first segment bounding polygon line
                            if( jdx .eq. nppoly ) then
                                jdx1 = 1
                            else
                                jdx1 = jdx + 1
                            end if
                            px(1) = polyx(polylist(jdx))
                            py(1) = polyy(polylist(jdx))
                            pz(1) = polyz(polylist(jdx))
                            px(2) = polyx(polylist(jdx1))
                            py(2) = polyy(polylist(jdx1))
                            pz(2) = polyz(polylist(jdx1))
                            px(3) = litx(polylist(jdx1))
                            py(3) = lity(polylist(jdx1))
                            pz(3) = litz(polylist(jdx1))
                            px(4) = litx(polylist(jdx))
                            py(4) = lity(polylist(jdx))
                            pz(4) = litz(polylist(jdx))
                            ! check with reference point for a good box
                            tin1 = inside(refpoint, 4, px, py, pz, 4, plist, ier)
                            if( (ier .gt. 0) .or. tin1 ) then
                                !write(*,*) "# inside error ", ier, "index ", jdx
                                cycle
                            else
                                ! check for inside (point already set above)
                                tin2 = inside(point, 4, px, py, pz, 4, plist, ier)
                                if( ier .gt. 0) then
                                    !call rtrans(4, px, py, pz, plat, plon)
                                    write(*,*) '#',ier,'error in inside routine. idx=',idx, 'idx1=', idx1
                                    cycle
                                else if( tin2 ) then
                                    ! find location of reflected point
                                    ! set outer box for reflected point location
                                    mx(1) = polyx(polylist(jdx))
                                    my(1) = polyy(polylist(jdx))
                                    mz(1) = polyz(polylist(jdx))
                                    mx(2) = polyx(polylist(jdx1))
                                    my(2) = polyy(polylist(jdx1))
                                    mz(2) = polyz(polylist(jdx1))
                                    mx(3) = bigx(polylist(jdx1))
                                    my(3) = bigy(polylist(jdx1))
                                    mz(3) = bigz(polylist(jdx1))
                                    mx(4) = bigx(polylist(jdx))
                                    my(4) = bigy(polylist(jdx))
                                    mz(4) = bigz(polylist(jdx))
                                    ! save point to be mirrored or moved
                                    savpoint = point
                                    if( bound_f ) then
                                        ! find point on boundary instead of using reflected point
                                        call tobound_point (point, px,py,pz, ier)
                                    else
                                        ! mirror point
                                        call mirror_point (point, px,py,pz, mx,my,mz, ier)
                                    end if
                                    if( ier .gt. 0) then
                                        ! error on return
                                        write(*,*) "# mirror point error ", ier
                                    end if

                                    ! copy new point into array
                                    rptx(idx) = point(1)
                                    rpty(idx) = point(2)
                                    rptz(idx) = point(3)
                                    if( debug ) then
                                        ! debug, write out boxes used to move point
                                        write(3,1015) 'B  5'
                                        call rtrans(4, px, py, pz, plat, plon)
                                        do kdx = 1, 4
                                            write(3,1020) plon(kdx)/sc, plat(kdx)/sc
                                        end do
                                        write(3,1020) plon(1)/sc, plat(1)/sc, idx
                                        call rtrans(4, mx, my, mz, mlat, mlon)
                                        write(3,1015) 'B  5'
                                        do kdx = 1, 4
                                            write(3,1020) mlon(kdx)/sc, mlat(kdx)/sc
                                        end do
                                        write(3,1020) mlon(1)/sc, mlat(1)/sc, idx
                                    end if
                                end if
                            end if
                        end do
                    end if
                end if
            end do
        end do
    end do

    if( debug ) then
        ! close box output file
        close(3)
    end if

    if( debug ) then
        ! debug, output polygon created smaller
        write(*,*) "# Small polygon created, writing file: windgen_int_small.pol"
        call rtrans(nppoly, litx, lity, litz, litlat, litlon)
        open(3,FILE='windgen_int_small.pol',STATUS='REPLACE')
        write(3,1000) 'B ', nppoly
        do idx = 1, nppoly
            write(3,1020) litlon(polylist(idx))/sc, litlat(polylist(idx))/sc
        end do
        close(3)
        ! debug, output polygon created larger
        write(*,*) "# Large polygon created, writing file: windgen_int_large.pol"
        call rtrans(nppoly, bigx, bigy, bigz, biglat, biglon)
        open(3,FILE='windgen_int_large.pol',STATUS='REPLACE')
        write(3,1000) 'B ', nppoly
        do idx = 1, nppoly
            write(3,1020) biglon(polylist(idx))/sc, biglat(polylist(idx))/sc
        end do
        close(3)
    end if

    ! expand trangulation array to receive mirrored points
    allocate (t1(1:nval), t2(1:nval), t3(1:nval), t4(1:nval), t5(1:nval), s1(1:nval) )
    ! whole array assignment
    t1 = rlat
    t2 = rlon
    t3 = x
    t4 = y
    t5 = z
    s1 = stanum
    deallocate (rlat, rlon, stanum, x, y, z)
    ! add to list to be processed
    eval = nval + mval
    allocate (rlat(1:eval), rlon(1:eval), stanum(1:eval), x(1:eval), y(1:eval), z(1:eval))
    ! whole array assignment
    do idx = 1, nval
        rlat(idx) = t1(idx)
        rlon(idx) = t2(idx)
        x(idx) = t3(idx)
        y(idx) = t4(idx)
        z(idx) = t5(idx)
        stanum(idx) = s1(idx)
    end do
    deallocate (t1, t2, t3, t4, t5, s1)

    ! output debugging files
    if( debug ) then
        ! output arc of mirrored point pairs to file (grass import format).
        open(100,FILE='windgen_int_reflected.edg',STATUS='REPLACE')
        call rtrans(nval, rptx, rpty, rptz, rptlat, rptlon)
    end if

    ! add valid mirrored points to triangulation array
    mval = 0
    do idx = 1, nval
        if( mir_flg(idx) ) then
            ! a reflected value has been created
            mval = mval + 1
            rlat(nval+mval) = rptlat(idx)/sc
            rlon(nval+mval) = rptlon(idx)/sc
            stanum(nval+mval) = stanum(idx)
            x(nval+mval) = rptx(idx)
            y(nval+mval) = rpty(idx)
            z(nval+mval) = rptz(idx)
            ! output debugging files
            if( debug ) then
                ! show point mirrored and reflected point as arc for grass import
                write(100,1010) 'L 2'
                write(100,1020) rlon(idx), rlat(idx), idx
                write(100,1020) rlon(nval+mval), rlat(nval+mval)
            end if
        end if
    end do

    !add format statements to eliminate space at beginning of line
    1000 format(a2,i7)
    1010 format(a3)
    1015 format(a4)
    1020 format(f20.15,1x,f20.15,1x,i5)

    ! close debugging files
    if( debug ) then
        close(100)
    end if

end if
    
! selection of interpolation type
if( invdist ) then
    ! use inverse distance weighting

    ! output debugging files
    if( debug ) then
        ! output file for points
        open(3,FILE='windgen_int_internal.pts',STATUS='REPLACE')
    end if

    ! loop through all interpolation points
    do idx = 1, ntr
        ! find invdist_pts closest points to interpolation point

        ! allocate arrays for return indexes and distances
        allocate (seln(1:invdist_pts), dist(1:invdist_pts), idw(1:invdist_pts), nidw(1:invdist_pts))

        ! copy point into point array
        point(1) = trx(idx)
        point(2) = try(idx)
        point(3) = trz(idx)

        call idfind(invdist_pts, point, eval, x,y,z, seln, dist)

        if (seln(invdist_pts) .gt. 0 ) then
            ! A complete list of points was returned
            ! calculate normalized inverse distance coefficients
            temp = 0.0
            do jdx = 1, invdist_pts
                ! find inverse distance
                if( dist(jdx) .gt. TINY(dist) ) then
                    idw(jdx) = 1.0 / dist(jdx)
                else
                    idw(jdx) = HUGE(dist)
                end if
                ! sum distances
                temp = temp + idw(jdx)
            end do
            do jdx = 1, invdist_pts
                ! normalize distance values
                nidw(jdx) = idw(jdx) / temp
            end do
            ! output points
            do jdx = 1, invdist_pts
                write(*,FMT="(1x,i9,1x,f11.8)",ADVANCE="NO") stanum(seln(jdx)), nidw(jdx)
            end do
            ! output carriage return
            write(*,*)
            ! list points on different output device
            do jdx = 1, invdist_pts
                write(6,*) stanum(seln(jdx)), nidw(jdx)
            end do
            if( debug ) then
                ! write points to file
                write(3,*) trlon(idx), trlat(idx)
                write(3,*)
                do jdx = 1, invdist_pts
                    write(3,*) rlon(seln(jdx)), rlat(seln(jdx)), stanum(seln(jdx)), dist(jdx), idw(jdx), nidw(jdx)
                end do
                write(3,*) rlon(seln(1)), rlat(seln(1)), stanum(seln(1)), dist(1), idw(1), nidw(1)
                write(3,*)
            end if
        else
            ! not enough points to do averaging
            write (*,*) '# not enough points found'
            stop
        end if
    end do

    ! output debugging files
    if( debug ) then
        close(3)
    end if

else
    ! use triangulation
    ! set array dimensions for triangulation based on number of values
    ntmx=2*eval
    nt6=6*eval
    lwk=2*eval

    ! Array storage for the triangulation, work space, and nodal
    ! coordinates.
    allocate (list(1:nt6), lptr(1:nt6), lend(1:eval), iwk(1:lwk))
    allocate (ds(1:eval))

    ! Array storage for the triangle list.
    allocate (ltri(1:nrow,1:ntmx))

    ! array storage for list of boundary nodes if needed
    allocate (bounds(1:eval))

    ! Create the triangulation and test the error flag.
    call trmesh(eval, x, y, z, list, lptr, lend, lnew, iwk, iwk(eval+1), ds, ier)
    if (ier .eq. -2) then
        write (*,*) '# error in trmesh: the first three nodes are collinear'
        stop
    else if (ier .gt. 0) then
        write (*,*) '# error in trmesh: duplicate nodes encountered'
        stop
    endif

    ! output debugging files
    if( debug ) then
        ! output point pairs to file (grass import format).
        open(3,FILE='windgen_int_internal.edg',STATUS='REPLACE')
        call trarcs (3,eval,rlon,rlat,stanum,list,lptr,lend,ier)
        close(3)
    end if

    ! output debugging files
    if( debug ) then
        ! output file for points
        open(3,FILE='windgen_int_internal.pts',STATUS='REPLACE')
    end if
    ! for all the points
    do idx = 1, ntr
        ! place the point within the triangulation and find points of the enclosing triangle
        point(1) = trx(idx)
        point(2) = try(idx)
        point(3) = trz(idx)
        call trfind (1, point, eval, x,y,z, list,lptr,lend, b1,b2,b3, i1,i2,i3)

        if (i3 .gt. 0 ) then
            ! Point is within a triangle
            ! calculate normalized barycentric coefficients
            temp = b1+b2+b3
            nb1 = b1 / temp
            nb2 = b2 / temp
            nb3 = b3 / temp
            ! output points
            !write(*,*) stanum(i1), nb1, stanum(i2), nb2, stanum(i3), nb3
            !write(6,*) stanum(i1), nb1
            !write(6,*) stanum(i2), nb2
            !write(6,*) stanum(i3), nb3

            !Changed above write statements to display to a fixed precision for OS/compiler consistency - LEW
            write(fmt_str_nb, '("(ES", i0, ".", i0, "e2)")') precision(nb1)+6, precision(nb1)

            write(*, '(6x, i6, 1x)', ADVANCE='NO') stanum(i1)
            write(*, fmt_str_nb, ADVANCE='NO') nb1
            write(*, '(6x, i6, 1x)', ADVANCE='NO') stanum(i2)
            write(*, fmt_str_nb, ADVANCE='NO') nb2
            write(*, '(6x, i6, 1x)', ADVANCE='NO') stanum(i3)
            write(*, fmt_str_nb, ADVANCE='YES') nb3

            write(*, '(6x, i6, 1x)', ADVANCE='NO') stanum(i1)
            write(*, fmt_str_nb, ADVANCE='YES') nb1
            write(*, '(6x, i6, 1x)', ADVANCE='NO') stanum(i2)
            write(*, fmt_str_nb, ADVANCE='YES') nb2
            write(*, '(6x, i6, 1x)', ADVANCE='NO') stanum(i3)
            write(*, fmt_str_nb, ADVANCE='YES') nb3

            ! output debugging files
            if( debug ) then
                ! write points to file
                write(3,*) trlon(idx), trlat(idx)
                write(3,*)
                write(3,*) rlon(i1), rlat(i1), stanum(i1)
                write(3,*) rlon(i2), rlat(i2), stanum(i2)
                write(3,*) rlon(i3), rlat(i3), stanum(i3)
                write(3,*) rlon(i1), rlat(i1), stanum(i1)
                write(3,*)
            end if
        else if (i2 .gt. 0) then
            ! Point is outside of triangulation, search for closest boundary node
            ! and two neighbors and find weighting coefficients
            ! find boundary nodes
            call bnodes (eval,list,lptr,lend, bounds,nbnd,narcs,ntri)
            ! find closest boundary node and return indexes and proportional coefficients
            call outside_weights( point, eval, x,y,z, nbnd, bounds, b1, b2, i1, i2, ier, 4, mx, my, mz)
            temp = b1+b2
            nb1 = b1 / temp
            nb2 = b2 / temp
            ! output points
            !write(*,*) stanum(i1), nb1, stanum(i2), nb2
            !write(6,*) stanum(i1), nb1
            !write(6,*) stanum(i2), nb2

            !Changed above write statements to display to a fixed precision for OS/compiler consistency - LEW
            write(fmt_str_nb, '("(ES", i0, ".", i0, "e2)")') precision(nb1)+6, precision(nb1)

            write(*, '(6x, i6, 1x)', ADVANCE='NO') stanum(i1)
            write(*, fmt_str_nb, ADVANCE='NO') nb1
            write(*, '(6x, i6, 1x)', ADVANCE='NO') stanum(i2)
            write(*, fmt_str_nb, ADVANCE='YES') nb2

            write(*, '(6x, i6, 1x)', ADVANCE='NO') stanum(i1)
            write(*, fmt_str_nb, ADVANCE='YES') nb1
            write(*, '(6x, i6, 1x)', ADVANCE='NO') stanum(i2)
            write(*, fmt_str_nb, ADVANCE='YES') nb2

            ! output debugging files
            if( debug ) then
                call rtrans(1, mx, my, mz, mlat, mlon)
                ! write points to file
                write(3,*) trlon(idx), trlat(idx)
                write(3,*)
                write(3,*) rlon(i1), rlat(i1), stanum(i1)
                write(3,*) rlon(i2), rlat(i2), stanum(i2)
                write(3,*)
                write(3,*) mlon(1)/sc, mlat(1)/sc
                write(3,*)
            end if
        else
            ! all points lie on the same great circle (or are coplanar)
            ! no output
            write (*,*) '# no triangulation possible'
            stop
        end if
    end do

    ! output debugging files
    if( debug ) then
        close(3)
    end if

end if

! find voronoi lines and write to output file (grass import format)
if( voronoi ) then

    ! allocate memory for arrays needed for voronoi calculations
    vrmx = 2*eval-4
    allocate ( listc(1:nt6), lbtri(1:6,1:eval) )
    allocate (xc(1:ntmx), yc(1:vrmx), zc(1:vrmx), rc(1:vrmx), vclon(1:vrmx), vclat(1:vrmx) )

!    call trmesh(eval, x, y, z, list, lptr, lend, lnew, iwk, iwk(eval+1), ds, ier)

    ! Note that the triangulation data structure is altered if NB > 0.
    ! for that reason, this is put after the triangulation is competed
    call crlist (eval, eval, x, y, z, list, lend, lptr, lnew, lbtri, listc, nbnd, xc, yc, zc, rc, ier)
    if (ier .gt. 0) then
        write (*,*) '# error in crlist: ier = ', ier
        stop
    endif

    ! translate the voronoi centers into latitude, longitude list
    call rtrans(vrmx, xc, yc, zc, vclat, vclon)
    ! convert values from radians to degrees
    do idx = 1, vrmx
        vclat(idx) = vclat(idx)/sc
        vclon(idx) = vclon(idx)/sc
    end do

    ! opened file for writing
    open(3,FILE='voronoi_all.edg',STATUS='REPLACE')

    ! output the voronoi lines
    call vrout (3,eval,vrmx,listc,lptr,lend,vclon,vclat,ier)

    close(3)
end if

close(6)
end