the dimensions of p and y in the calling program can be different
from the dimensions in this subroutine. for instance,
real*10 p(4,2),y(7,3,3).
in this case the calling statement might look like
x=terpf(p(1,m),y(1,1,n))
positions 6-7 are ignored here, since this is the 10-point
interpolator, but they are used for the 14-point interpolator
NOTE: if these unused positions are filled with zeroes, then
a call to TERP14 should give exactly the same result