      subroutine interpolate(childgrid,order,dim,refinerat,xx,yy,zz)
      use types
      implicit none
     
      real(SP), dimension(:,:) :: childgrid
      integer(I4B)               :: order, dim, refinerat, step
      integer(I4B)               :: i,j,k,no_of_cols,no_of_rows, err, m
      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

      write(6,*) 'within interpolate'
      
      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_cols=size(childgrid,dim=1)
      no_of_rows=size(childgrid,dim=2)

      step=refinerat  

c interpolate along all the rows first

      if(order.eq.3) then
      no_of_cols=no_of_cols-step
      endif 
      write(6,*) 'step=', step,'refinerat=', refinerat


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

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

          do m=1,refinerat-1
          call polint(xtmp,ytmp,xx(j+m),childgrid(j+m,i),dy)
          enddo
         
         enddo
      enddo

c now interpolate along the columns
      if(order.eq.3) then
      no_of_rows=no_of_rows-step
      endif 

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

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

          do m=1,refinerat-1
          call polint(xtmp,ytmp,yy(j+m),childgrid(i,j+m),dy)
          enddo
         
         enddo
      enddo 

      end subroutine interpolate



