      subroutine interpolate(childgrid,order,dim,refinerat,xx,yy,zz)
      use types
      implicit none
     
      real(SP), dimension(:,:,:) :: childgrid
      integer(I4B)               :: order, dim, refinerat,no_of_planes
      integer(I4B)               :: i,j,k,no_of_cols,no_of_rows, err
      real(SP), dimension(:)     :: xx,yy,zz     
      real(SP), dimension(:), allocatable :: xtmp, ytmp
      real(SP)::dy

       
      INTERFACE 

      subroutine polint(xa,ya,x,y,dy)
      use types
      implicit none
      real(SP), dimension(:), intent(in) :: xa,ya
      real(SP), intent(in)::x
      real(SP), intent(out) ::y,dy
      end subroutine polint

      END INTERFACE
      
      if (order.eq.3) then
      allocate(xtmp(3),ytmp(3), stat=err)
      else 
      allocate(xtmp(2),ytmp(2), stat=err)
      endif
      if(err.ne.0) then
      write(6,*) 'interpolate:Allocation failed'
      stop
      endif

 


      no_of_planes=size(childgrid,dim=3)

      do k=1,no_of_planes,2

! interpolate along all the rows first
      no_of_rows=size(childgrid,dim=2)
      no_of_cols=size(childgrid,dim=1)
      if(order.eq.3) then
      no_of_cols=no_of_cols-2
      endif 

      do i=1,no_of_rows,2
         do j=1,no_of_cols-1,2
 
! Initialise the x and y values for 1d interpolation
! and call the 1d interpolation routine

          xtmp(1)=xx(j)
          xtmp(2)=xx(j+2)
          ytmp(1)=childgrid(j,i,k)
          ytmp(2)=childgrid(j+2,i,k)
          if(order.eq.3) then
          xtmp(3)=xx(j+4)
          ytmp(3)=childgrid(j+4,i,k)
          endif

   
          call polint(xtmp,ytmp,xx(j+1),childgrid(j+1,i,k),dy)
         
         enddo
      enddo

! now interpolate along the columns

      if(order.eq.3) then
      no_of_rows=no_of_rows-2
      endif 


      do i=1,no_of_cols
         do j=1,no_of_rows-1,2
 
! Initialise the x and y values for 1d interpolation
! and call the 1d interpolation routine

          xtmp(1)=yy(j)
          xtmp(2)=yy(j+2)
          ytmp(1)=childgrid(i,j,k)
          ytmp(2)=childgrid(i,j+2,k)
          if(order.eq.3) then
          xtmp(3)=yy(j+4)
          ytmp(3)=childgrid(i,j+4,k)
          endif

   
          call polint(xtmp,ytmp,yy(j+1),childgrid(i,j+1,k),dy)
         
         enddo
      enddo 

      enddo ! do k=1,n_of_planes,2

!     Now that we have all the alternate planes available
!     interpolate along the planes to find the planes in between
 
       if(order.eq.3) then
       no_of_planes=no_of_planes-2
       endif

      do i=1, no_of_rows
         do j=1, no_of_cols
           do k=1, no_of_planes-1,2

          xtmp(1)=zz(k)
          xtmp(2)=zz(k+2)
          ytmp(1)=childgrid(i,j,k)
          ytmp(2)=childgrid(i,j,k+2)
          if(order.eq.3) then
          xtmp(3)=zz(k+4)
          ytmp(3)=childgrid(i,j,k+4)
          endif

   
          call polint(xtmp,ytmp,zz(k+1),childgrid(i,j,k+1),dy)
          enddo
        enddo
       enddo

      deallocate(xtmp)
      deallocate(ytmp) 
      end subroutine interpolate






