find_y_nonuniform_x Function

  function find_y_nonuniform_x(x0,x,y,n)

    !! Routine to find y0 such that y0 = y(x0) given a set of
    !! values x(1:n), y(1:n)
    !! author: P J Knight, CCFE, Culham Science Centre
    !! x0 : input real : x value at which we want to find y
    !! x(1:n) : input real array : monotonically de/increasing x values
    !! y(1:n) : input real array : y values at each x
    !! n : input integer : size of array
    !! This routine performs a simple linear interpolation method
    !! to find the y value at x = x0. If x0 lies outside the
    !! range [x(1),x(n)], the y value at the nearest 'end' of the data
    !! is returned.
    !! None
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    real(dp) :: find_y_nonuniform_x

    !  Arguments

    integer, intent(in) :: n
    real(dp), intent(in) :: x0
    real(dp), dimension(n), intent(in) :: x
    real(dp), dimension(n), intent(in) :: y

    !  Local variables
    integer :: i,j

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Step through arrays until x crosses the value of interest

    j = 0
    rough_search: do i = 1,n-1

       if (((x(i)-x0)*(x(i+1)-x0)) <= 0.0D0) then
          j = i
          exit rough_search
       end if

    end do rough_search

    if (j /= 0) then

       !  Simply do a linear interpolation between the two grid points
       !  spanning the point of interest

       find_y_nonuniform_x = y(j) + (y(j+1)-y(j))*(x0-x(j))/(x(j+1)-x(j))

    else  !  No points found, so return the 'nearest' y value

       if (x(n) > x(1)) then  !  values are monotonically increasing
          if (x0 > x(n)) then
             find_y_nonuniform_x = y(n)
          else
             find_y_nonuniform_x = y(1)
          end if

       else  !  values are monotonically decreasing
          if (x0 < x(n)) then
             find_y_nonuniform_x = y(n)
          else
             find_y_nonuniform_x = y(1)
          end if

       end if

    end if

  end function find_y_nonuniform_x