      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
      real(SP), dimension(:)     :: xx,yy,zz     


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

      real(SP) , dimension(size(childgrid,dim=2),2) :: ho
      real(SP), dimension(size(childgrid, dim =2))::dy1, den
      real(SP) :: dy
      integer(I4B)               :: i,j,k,no_of_cols,no_of_rows, err,l
      integer(I4B)               :: step, m


!HPF$ ALIGN (:,*) with childgrid(*,:,*) :: ho
!HPF$ ALIGN (:) with childgrid(*,:,*) :: dy1,den


      write(6,*) 'within interpolate'


      step = refinerat
      no_of_planes=size(childgrid,dim=3)
      no_of_rows=size(childgrid,dim=2)
      no_of_cols=size(childgrid,dim=1)

      ho(:,1) = xx(1)-xx(2)
      ho(:,2) = xx(3)-xx(2)




      do k=1,no_of_planes,step

      ho(:,1) = xx(1)-xx(2)
      ho(:,2) = xx(3)-xx(2)

       do j=1,no_of_cols-1,step

      den(:) = ho(:,1) - ho(:,2)
c!      write(6,*) 'den=' ,den

      forall(i=1:no_of_rows:step) childgrid(j+1,i,k) = 
     &                                    childgrid(j,i,k)


      forall(i=1:no_of_rows:step) den(i)=(
     &                     childgrid(j+step,i,k)-
     &                     childgrid(j,i,k)) / den(i)   
c!      write(6,*) 'after 1st forall'

      forall(i=1:no_of_rows:step) dy1(i)=ho(i,1) * den(i)
c!      write(6,*) 'after 2nd forall'

      forall(i=1:no_of_rows:step) childgrid(j+1,i,k) = 
     &                                    childgrid(j+1,i,k) + dy1(i)
c!      write(6,*) 'after 3rd forall'


       enddo ! j

      ho(:,1) = yy(1)-yy(2)
      ho(:,2) = yy(3)-yy(2)


      do j=1, no_of_rows-1,step

      den(:) = ho(:,1) - ho(:,2)
c!      write(6,*) 'den=' ,den

      forall(i=1:no_of_cols) childgrid(i,j+1,k) = 
     &                                    childgrid(i,j,k)
      forall(i=1:no_of_cols) den(i)=(childgrid(i,j+step,k)-
     &                             childgrid(i,j,k)) / den(i)   
c!      write(6,*) 'after 1st forall'

      forall(i=1:no_of_cols) dy1(i)=ho(i,1) * den(i)
c!      write(6,*) 'after 2nd forall'

      forall(i=1:no_of_cols) childgrid(i,j+1,k) = 
     &                                    childgrid(i,j+1,k) + dy1(i)
c!      write(6,*) 'after 3rd forall'



        enddo ! j=1, no_of_rows-1
       enddo ! k

c!      write(6,*) childgrid

     
      do i = 1, no_of_rows
       ho(:,1) = zz(1)-zz(2)
       ho(:,2) = zz(3)-zz(2)
          
        do k = 1, no_of_planes-1, step

      den(:) = ho(:,1) - ho(:,2)
      forall(j=1:no_of_cols) childgrid(i,j,k+1) = 
     &                                    childgrid(i,j,k)

c!      write(6,*) 'den=' ,den
      write(6,*) k,k+step
      forall(j=1:no_of_cols) den(j)=(childgrid(i,j,k+step)-
     &                             childgrid(i,j,k)) / den(j)   
c!      write(6,*) 'after 1st forall'

      forall(j=1:no_of_cols) dy1(j)=ho(j,1) * den(j)
c!      write(6,*) 'after 2nd forall'

      forall(j=1:no_of_cols) childgrid(i,j,k+1) = 
     &                                    childgrid(i,j,k+1) + dy1(j)
c!      write(6,*) 'after 3rd forall'



        enddo ! i= 1, no_of_rows
      enddo ! k= 1, no_of_planes-1

      end subroutine interpolate























