interp_from_grid_linear Subroutine

  subroutine interp_from_grid_linear(igrid,xp,wp)
    ! get point values from given grid via linear interpolation
    ! multi-D method
    
    integer :: igrid
    double precision :: xp(1:ndim),wp(1:nw)

    double precision :: dxb1,dxb2,dxb3,xbmin1,xbmin2,xbmin3,xbmax1,xbmax2,&
       xbmax3
    integer :: inblock,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,j
    integer :: ixb1,ixb2,ixb3,ixi1,ixi2,ixi3
    double precision :: xc(0:1,0:1,0:1,1:ndim),wc(0:1,0:1,0:1,nw)

    ! block/grid boundaries
    xbmin1=rnode(rpxmin1_,igrid)
    xbmin2=rnode(rpxmin2_,igrid)
    xbmin3=rnode(rpxmin3_,igrid)
    xbmax1=rnode(rpxmax1_,igrid)
    xbmax2=rnode(rpxmax2_,igrid)
    xbmax3=rnode(rpxmax3_,igrid)

    ! whether or not next point is in this block/grid
    inblock=0
    if (xp(1)>=xbmin1 .and. xp(1)<xbmax1) inblock=inblock+1
    if (xp(2)>=xbmin2 .and. xp(2)<xbmax2) inblock=inblock+1
    if (xp(3)>=xbmin3 .and. xp(3)<xbmax3) inblock=inblock+1
    if (inblock/=ndim) then
      call mpistop('Interpolation error: given point is not in given grid')
    endif

    ! cell indexes for the point
    dxb1=rnode(rpdx1_,igrid)
    dxb2=rnode(rpdx2_,igrid)
    dxb3=rnode(rpdx3_,igrid)
    ixOmin1=ixmlo1
    ixOmin2=ixmlo2
    ixOmin3=ixmlo3
    ixb1=floor((xp(1)-ps(igrid)%x(ixOmin1,ixOmin2,ixOmin3,1))/dxb1)+ixOmin1
    ixb2=floor((xp(2)-ps(igrid)%x(ixOmin1,ixOmin2,ixOmin3,2))/dxb2)+ixOmin2
    ixb3=floor((xp(3)-ps(igrid)%x(ixOmin1,ixOmin2,ixOmin3,3))/dxb3)+ixOmin3

    ! nearby cells for interpolation
    do ixi1=0,1
    do ixi2=0,1
    do ixi3=0,1
      xc(ixi1,ixi2,ixi3,:)=ps(igrid)%x(ixb1+ixi1,ixb2+ixi2,ixb3+ixi3,:)
      wc(ixi1,ixi2,ixi3,:)=ps(igrid)%w(ixb1+ixi1,ixb2+ixi2,ixb3+ixi3,:)
    enddo
    enddo
    enddo

    call interp_linear_multiD(xc,xp,wc,wp,nw)

  end subroutine interp_from_grid_linear