subroutine polin2(x1a,x2a,ya,x1,x2,y,dy)
real(8), dimension(:), intent(in) :: x1a,x2a
real(8), dimension(:,:), intent(in) :: ya
real(8), intent(in) :: x1,x2
real(8), intent(out) :: y,dy
integer :: j,m,ndum
real(8), dimension(size(x1a)) :: ymtmp
real(8), dimension(size(x2a)) :: yntmp
m = size(x1a);if(m/=size(ya,1))stop "POLINT: wrong dimensions m"
ndum=size(x2a);if(ndum/=size(ya,2))stop "POLINT: wrong dimensions ndum"
do j=1,m
yntmp=ya(j,:)
call polint(x2a,yntmp,x2,ymtmp(j),dy)
end do
call polint(x1a,ymtmp,x1,y,dy)
end subroutine polin2