subroutine idfind(nst, p, n, x, y, z, seln, dist)
integer nst, n, seln(*)
real    p(3), x(n), y(n), z(n), dist(*)

!***********************************************************
!
!                                                   Fred Fox
!                                       USDA-ARS-CGAHR-EWERU
!                                               Manhattan KS
!                                      fred.fox@ars.usda.gov
!                                                  6/03/2011
!
!   This subroutine locates the NST points that are closest to
! point P from the set of N points supplied and returns the point
! index and the distances for each point selected.
!
! On input:
!
!       NST = Number of points to be returned
!
!       P = Array of length 3 containing the x, y, and z
!           coordinates (in that order) of the point P to be
!           located.
!
!       N = Number of points to be searched.  N .GE. NST.
!
!       X,Y,Z = Arrays of length N containing the Cartesian
!               coordinates of the triangulation nodes (unit
!               vectors).  (X(I),Y(I),Z(I)) defines node I
!               for I = 1 to N.
!
! Input parameters are not altered by this routine.
!
! On output:
!
!       SELN = Array of length NST of indexes of the nearest points
!
!       DIST = Array of length NST with distance from point P
!              to point reference by SELN.
!
! Modules required by IDFIND:  arc_len
!
! Intrinsic function called by TRFIND:  None
!
!***********************************************************

integer idx, jdx, imaxdist
real dist1, maxdist, dx(2), dy(2), dz(2)

! function declaration
real arc_len

! set up location point
dx(1) = p(1)
dy(1) = p(2)
dz(1) = p(3)

! set up initial distance array to large values
do idx = 1,nst
    dist(idx) = 10.0
    seln(idx) = 0
end do
! distances are in radians so should be max 2pi
! this value triggers replacement of originals
maxdist = 10.0
imaxdist = 1

! loop through all candidate points selecting the closest ones
do idx = 1, n

    ! insert point into distance array
    dx(2) = x(idx)
    dy(2) = y(idx)
    dz(2) = z(idx)

    dist1 = arc_len(dx, dy, dz)

    !write(*,*) 'imaxdist, maxdist, dist1:', imaxdist, maxdist, dist1
    if( dist1 .lt. maxdist ) then
        ! this point is closer than the max distance point, replace it.
        dist(imaxdist) = dist1
        seln(imaxdist) = idx
        ! find point in selected array that is new max distance
        do jdx = 1, nst
            if( dist(jdx) .gt. dist(imaxdist) ) then
                maxdist = dist(jdx)
                imaxdist = jdx
            end if
        end do
    end if

end do
return
end