Below is a function that interpolates an array using a Piecewise Cubic Hermitian Interpolator. The required external methods are provided by the PCHIP module.
function interpolate(x, y, p) result(r)
!! This function constructs a piecewise cubic Hermitian interpolation of an array y(x) based on discrete numerical data,
!! and evaluates the interpolation at points p. Note that the mesh spacing of x does not necessarily have to be uniform.
real(wp), intent(in) :: x(:) !! Variable x
real(wp), intent(in) :: y(size(x)) !! Function y(x)
real(wp), intent(in) :: p(:) !! Interpolation domain p
real(wp) :: r(size(p)) !! Interpolation result y(p)
real(wp) :: d(size(x))
integer :: err
! Create a PCHIP interpolation of the input data
call dpchez(size(x), x, y, d, .false., 0, 0, err)
! Extract the interpolated data at provided points
call dpchfe(size(x), x, y, d, 1, .false., size(p), p, r, err)
end function