get_point_w_ingrid Subroutine

  subroutine get_point_w_ingrid(igrid,xp,wp,variable_type)
    double precision :: xp(1:ndim),wp(1:nw)
    character(*) :: variable_type
    integer, intent(in) :: igrid

    double precision :: dxb1,dxb2,dxb3,xd1,xd2,xd3,xbmin1,xbmin2,xbmin3,xbmax1,&
       xbmax2,xbmax3,temp
    integer :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,ixOmin1,ixOmin2,&
       ixOmin3,ixOmax1,ixOmax2,ixOmax3
    integer :: ixbl1,ixbl2,ixbl3,ix1,ix2,ix3,ixAmin1,ixAmin2,ixAmin3,ixAmax1,&
       ixAmax2,ixAmax3,j,ingrid
    double precision :: factor(0:1,0:1,0:1)

    ixImin1=ixGlo1;ixImin2=ixGlo2;ixImin3=ixGlo3;ixImax1=ixGhi1
    ixImax2=ixGhi2;ixImax3=ixGhi3;
    ixOmin1=ixMlo1;ixOmin2=ixMlo2;ixOmin3=ixMlo3;ixOmax1=ixMhi1
    ixOmax2=ixMhi2;ixOmax3=ixMhi3;
    xbmin1=rnode(rpxmin1_,igrid)-rnode(rpdx1_,igrid)
    xbmin2=rnode(rpxmin2_,igrid)-rnode(rpdx2_,igrid)
    xbmin3=rnode(rpxmin3_,igrid)-rnode(rpdx3_,igrid);
    xbmax1=rnode(rpxmax1_,igrid)+rnode(rpdx1_,igrid)
    xbmax2=rnode(rpxmax2_,igrid)+rnode(rpdx2_,igrid)
    xbmax3=rnode(rpxmax3_,igrid)+rnode(rpdx3_,igrid);
    ingrid=0
    if (xp(3)>xbmin3 .and. xp(3)<xbmax3) ingrid=ingrid+1
    if (xp(2)>xbmin2 .and. xp(2)<xbmax2) ingrid=ingrid+1
    if (xp(1)>xbmin1 .and. xp(1)<xbmax1) ingrid=ingrid+1

    if (ingrid==ndim) then
      dxb1=rnode(rpdx1_,igrid)
      dxb2=rnode(rpdx2_,igrid)
      dxb3=rnode(rpdx3_,igrid)
      ixbl1=floor((xp(1)-ps(igrid)%x(ixOmin1,ixOmin2,ixOmin3,1))/dxb1)+ixOmin1
      ixbl2=floor((xp(2)-ps(igrid)%x(ixOmin1,ixOmin2,ixOmin3,2))/dxb2)+ixOmin2
      ixbl3=floor((xp(3)-ps(igrid)%x(ixOmin1,ixOmin2,ixOmin3,3))/dxb3)+ixOmin3
      xd1=(xp(1)-ps(igrid)%x(ixbl1,ixbl2,ixbl3,1))/dxb1
      xd2=(xp(2)-ps(igrid)%x(ixbl1,ixbl2,ixbl3,2))/dxb2
      xd3=(xp(3)-ps(igrid)%x(ixbl1,ixbl2,ixbl3,3))/dxb3
      ixAmin1=ixbl1
      ixAmin2=ixbl2
      ixAmin3=ixbl3
      ixAmax1=ixbl1+1
      ixAmax2=ixbl2+1
      ixAmax3=ixbl3+1

      do ix1=0,1
      do ix2=0,1
      do ix3=0,1
        factor(ix1,ix2,ix3)=abs(1-ix1-xd1)*abs(1-ix2-xd2)*abs(1-ix3-xd3)
      enddo
      enddo
      enddo

      select case(variable_type)
        case('primitive')
          call phys_to_primitive(ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
             ixImax3,ixAmin1,ixAmin2,ixAmin3,ixAmax1,ixAmax2,ixAmax3,&
             ps(igrid)%w,ps(igrid)%x)
        case('conserved')

        case default
          call mpistop("get_point_w_ingrid: Unknown variable type!")
      end select

      wp=0.d0
      do ix1=0,1
      do ix2=0,1
      do ix3=0,1
        do j=1,nw
          wp(j)=wp(j)+factor(ix1,ix2,ix3)*ps(igrid)%w(ixbl1+ix1,ixbl2+ix2,&
             ixbl3+ix3,j)
        enddo
      enddo
      enddo
      enddo

      select case(variable_type)
        case('primitive')
        call phys_to_conserved(ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
           ixAmin1,ixAmin2,ixAmin3,ixAmax1,ixAmax2,ixAmax3,ps(igrid)%w,&
           ps(igrid)%x)
        case('conserved')

        case default
          call mpistop("get_point_w_ingrid: Unknown variable type!")
      end select

      if(physics_type=='mhd') then
        wp(iw_mag(1):iw_mag(ndir))=0.d0
        do ix1=0,1
        do ix2=0,1
        do ix3=0,1
          do j=1,ndir
            if (b0field) then
              temp=ps(igrid)%w(ixbl1+ix1,ixbl2+ix2,ixbl3+ix3,&
                 iw_mag(j))+ps(igrid)%B0(ixbl1+ix1,ixbl2+ix2,ixbl3+ix3,j,0)
              wp(iw_mag(j))=wp(iw_mag(j))+factor(ix1,ix2,ix3)*temp
            else
              temp=ps(igrid)%w(ixbl1+ix1,ixbl2+ix2,ixbl3+ix3,iw_mag(j))
              wp(iw_mag(j))=wp(iw_mag(j))+factor(ix1,ix2,ix3)*temp
            endif
          enddo
        enddo
        enddo
        enddo
      endif
    else
      call mpistop("get_point_w_ingrid: The point is not in given grid!")
    endif

  end subroutine get_point_w_ingrid