      subroutine interpolate(childgrid,order,dim,refinerat,xx,yy,zz)
      use types
      use parameters
      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
      integer(I4B)               :: step, m
      real(SP), dimension(:)     :: xx,yy,zz     
      real(SP), dimension(:), allocatable     :: xtmp, ytmp
      real(SP)::dy

!!HPF$ PROCESSORS P(PX,PY,PZ)
!!HPF$ DISTRIBUTE *(BLOCK,BLOCK,BLOCK) onto *P :: childgrid
!!HPF$ ALIGN xx(:) with *childgrid(:,*,*)
!!HPF$ ALIGN yy(:) with *childgrid(*,:,*)
!!HPF$ ALIGN zz(:) with *childgrid(*,*,:)
!!!!HPF$ ALIGN xx(*) with xtmp(*)
!!!HPF$ ALIGN yy(*) with xtmp(*)
!!!HPF$ ALIGN zz(*) with xtmp(*)
!!!HPF$ ALIGN ytmp(*) with childgrid(*,*,*)
       
      INTERFACE 

      subroutine polint(xa,ya,x,y,dy)
      use types
      use parameters
      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

 

      step = refinerat
      no_of_planes=size(childgrid,dim=3)

      do k=1,no_of_planes,step

!      forall(k=1:no_of_planes:step)

c 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-step
      endif 

      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,k)
          ytmp(2)=childgrid(j+step,i,k)
          if(order.eq.3) then
          xtmp(3)=xx(j+2*step)
          ytmp(3)=childgrid(j+2*step,i,k)
          endif

          do m=1, refinerat-1   
!          forall(m=1:refinerat-1)
          call polint(xtmp,ytmp,xx(j+m),childgrid(j+m,i,k),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,k)
          ytmp(2)=childgrid(i,j+step,k)
          if(order.eq.3) then
          xtmp(3)=yy(j+2*step)
          ytmp(3)=childgrid(i,j+2*step,k)
          endif

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

      enddo ! do k=1,n_of_planes,2

!!!!!      end forall

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

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

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

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

          enddo
        enddo
       enddo

      deallocate(xtmp)
      deallocate(ytmp) 
      end subroutine interpolate












