!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!   INTERPOLATE
!!!
!!!   Simple linear interpolation.  The routine
!!!   coarse_to_fine MUST be called before calling this routine, as the
!!!   interpolated values from the parent grid are assumed to exist
!!!   within the input CHILDGRID.
!!!   Warning: this routine assumes refinerat is 2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      subroutine interpolate(grid,order,dim,refinerat,xx,yy,zz,nx,ny,nz)

      use types
      use bh02

      real(SP), dimension(:,:,:)  :: grid
      integer(I4B)                :: order,dim,refinerat
      real(SP), dimension(:)      :: xx,yy,zz
      integer(I4B)                :: nx,ny,nz
!     LOCAL VARIABLES
      integer(I4B)                :: i,j,k,rr

      rr = refinerat
      if(rr > 2) then
         write(*,*) 'Error: interpolate: refinerat does not equal 2! '
         stop
      endif

!  LEGS
      
      do i=2,nx-1,rr
      do j=1,ny,rr
      do k=1,nz,rr
         grid(i,j,k) = (grid(i+1,j,k) + grid(i-1,j,k))/2.0
      enddo
      enddo
      enddo

      do i=1,nx,rr
      do j=2,ny-1,rr
      do k=1,nz,rr
         grid(i,j,k) = (grid(i,j+1,k) + grid(i,j-1,k))/2.0
      enddo
      enddo
      enddo
      
      do i=1,nx,rr
      do j=1,ny,rr
      do k=2,nz-1,rr
         grid(i,j,k) = (grid(i,j,k+1) + grid(i,j,k-1))/2.0
      enddo
      enddo
      enddo

!  FACES

      do i=2,nx-1,rr
      do j=2,ny-1,rr
      do k=1,nz,rr
         grid(i,j,k) = (grid(i+1,j+1,k) + grid(i+1,j-1,k) + &
                        grid(i-1,j+1,k) + grid(i-1,j-1,k) )/4.0
      enddo
      enddo
      enddo

      do i=2,nx-1,rr
      do j=1,ny,rr
      do k=2,nz-1,rr
         grid(i,j,k) = (grid(i+1,j,k+1) + grid(i+1,j,k-1) + &
                        grid(i-1,j,k+1) + grid(i-1,j,k-1) )/4.0
      enddo
      enddo
      enddo

      do i=1,nx,rr
      do j=2,ny-1,rr
      do k=2,nz-1,rr
         grid(i,j,k) = (grid(i,j+1,k+1) + grid(i,j+1,k-1) + &
                        grid(i,j-1,k+1) + grid(i,j-1,k-1) )/4.0
      enddo
      enddo
      enddo

!  CUBES

      do i=2,nx-1,rr
      do j=2,ny-1,rr
      do k=2,nz-1,rr
         grid(i,j,k) = (grid(i+1,j+1,k+1) + grid(i+1,j+1,k-1) + &
                        grid(i+1,j-1,k+1) + grid(i+1,j-1,k-1) + &
                        grid(i-1,j+1,k+1) + grid(i-1,j+1,k-1) + &
                        grid(i-1,j-1,k+1) + grid(i-1,j-1,k-1) )/8.0

      enddo
      enddo
      enddo
      

!     write(6,*) '------- INTERPOLATION DONE ------------'
      return
      end subroutine interpolate


