subroutine inject(level,id,p,pp)
use hpfprocessors
use param
integer level,id
real, dimension(:,:,:,:) :: p
real, dimension(:,:,:,:) :: pp
!HPF$ DISTRIBUTE p(*,BLOCK,*,*)
!HPF$ DISTRIBUTE pp(*,BLOCK,*,*)

integer offsetx,offsety,offsetz
integer nx,ny,nz    !number of grid points on this level
integer nxc,nyc,nzc !number of gird points on parent level
type(ADAPTIVE_MESH), pointer :: grid

  grid=>lp(level)%grid
  if(id <> 1) then
    do while (grid%grid_id <> id)
     grid=>grid%neighbor
    enddo
  endif

  nx=grid%nr(1)-1
  ny=grid%nr(2)-2
  nz=grid%nr(3)-1
  
  offsetx=(grid%rmin(1)-grid%parent%rmin(1))*rho+1
  offsety=(grid%rmin(2)-grid%parent%rmin(2))*rho+2
  offsetz=(grid%rmin(3)-grid%parent%rmin(3))*rho+1

  nxc=offsetx+(nx+1)/rho-1
  nyc=offsety+(ny+1)/rho-2
  nzc=offsetz+(nz+1)/rho-1


!write(6,*) 'offset used in inject',offsety,grid%rmin(2),grid%parent%rmin(2)

!transfer data points from fine grid to coarse one

!sanity check for the wave equation
if(any(abs(pp(:,offsetx:nxc,offsety:nyc,offsetz:nzc)- &
            p(:,1:nx:rho,1:ny:rho,1:nz:rho)) > 0.001)) then
   write(6,*) 'something fishy going on here!'
   stop
endif

 pp(:,offsetx:nxc,offsety:nyc,offsetz:nzc)= & 
 p(:,1:nx:rho,1:ny:rho,1:nz:rho)

i=max(offsety-2,1)
j=i+3
write(6,*) 'begin'
write(6,*) pp(1,5,i:j,1)
write(6,*) 'middle'
write(6,*) p(1,5,1:4,1)
write(6,*) 'end'
j=min(nyc+2,size(pp,dim=3))
i=j-3
write(6,*) pp(1,5,i:j,1)

return
end

