subroutine boundaries(t,curlev,nx,ibound,id,p,pp,pass)
use hpfprocessors
use param
real t
integer curlev,nx,id,pass
integer, dimension(3,2) :: ibound
real, dimension(:,:,:,:) :: p
real, dimension(:,:,:,:) :: pp

!HPF$ DISTRIBUTE p(*,BLOCK,*,*)
!HPF$ DISTRIBUTE pp(*,BLOCK,*,*)

type(ADAPTIVE_MESH),pointer :: grid
integer offsetx,offsety,offsetz

!---------------------------------------------------
! set boundary conditions after grid update
!---------------------------------------------------

write(6,*) '---- boundaries1, pass',pass,' ------'

! if ibound==1 then physical boundaries
if(ibound(1,1)==1) then
  p(1,1,:,1)=p(1,2,:,1)
else
 write(6,*) 'interpolation needed (1)'
 stop
endif
 
if(ibound(1,2)==1) then
  p(1,nx,:,1)=p(1,2,:,1)
else
 write(6,*) 'interpolation needed (2)'
 stop
endif


if(ibound(2,1)==1) then 
p(1,:,1,1)=sin(0.15*t)
else

if(pass==1) then

write(6,*) 'interpolation'
!here I assumed that there are two steps on fine grid;
!first step : boundary conditions are just interpolated when regridding
!             that is, this routine is called after first update
!second step: start with values that are average between retarded
!             and advanced solution on level-1
!             this is first call
!second call: who cares, the result will be overwritten by regridding
!             anyway. What about transfer????    
!
!write(6,*) 'insert interpolation here!!!'

!first transfer averaged data from parent to child

!compute offsets

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

  nnx=grid%nr(1)
  ny=grid%nr(2)
  nz=grid%nr(3)

if(nnx <> nx) then
  write(6,*) 'BOUNDARIES: something wrong with nx'
  stop
endif
  
  offsetx=(grid%rmin(1)-grid%parent%rmin(1))/lp(curlev-1)%dr+1
  offsety=(grid%rmin(2)-grid%parent%rmin(2))/lp(curlev-1)%dr+1
  offsetz=(grid%rmin(3)-grid%parent%rmin(3))/lp(curlev-1)%dr+1

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

!transfer data points from coarse grid to fine one
!for all x and z, and only y=1 (offsety)

  p(1,1:nx:rho,1,1:nz:rho)=  &
  0.5*( pp(1,offsetx:nxc,offsety,offsetz:nzc)+ & 
        pp(2,offsetx:nxc,offsety,offsetz:nzc))

!now interpolate

  if(rho == 2) then
    if(mod(nx,2) <> 1 .or. mod(ny,2) <> 1 .or. mod(nz,2) <> 1) then
     write(6,*) 'it does not work for even size'
     stop
    endif

!    interpolate each k plane (along x axis only)
  do k=1,nz,2
    forall(i=2:nx-1:2,j=1:ny:2)  &
     p(1,i,j,k)=0.5*(p(1,i-1,j,k)+p(1,i+1,j,k))
    
!    forall(i=1:nx,j=2:ny:2) &
!     p(:,i,j,k)=0.5*(p(:,i,j-1,k)+p(:,i,j+1,k))
  enddo

!   interpolate between k planes
     if(nz > 1) then
!      forall(i=1:nx,j=1:ny,k=2:nz:2) &
!      p(:,i,j,k)=0.5*(p(:,i,j,k-1)+p(:,i,j,k+1)) 
      forall(i=1:nx,k=2:nz:2) &
      p(1,i,1,k)=0.5*(p(1,i,1,k-1)+p(1,i,1,k+1)) 
      endif
  else
   write(6,*) 'must supply interpolation routine'
   stop 
  endif

endif
endif

!do not care about other boundaries

return
end



