      MODULE types
      
c    symbolic names for 4, 2, and 1 byte integers
	
	INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
	INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
	INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)

! and for single and double precision reals
	INTEGER, PARAMETER :: SP = KIND(1.0)
	INTEGER, PARAMETER :: DP = KIND(1.0D0)

      END MODULE types


      module parameters
      integer, parameter  :: PX=2,PY=2 !number of processors
      integer, parameter  :: nn = 99
      end module parameters

      program main
      use types
      use parameters
      real(SP), dimension(nn,nn,nn)::child
      real(SP), dimension(nn)    :: xx,yy,zz
      integer(I4B)              :: order, dim, refinerat
      real t1(PX*PY), t2(PX*PY)

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

!     interpolation
      real(SP) , dimension(nn,2) :: ho
      real(SP), dimension(nn)    ::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 child(*,:,*) :: ho
!HPF$ ALIGN (:) with child(*,:,*) :: dy1,den
!     end interpolation

      INTERFACE
      extrinsic(hpf_local) subroutine timer(return_time, initial_time)
      real, dimension(:), intent(in)  :: initial_time
      real, dimension(:), intent(out) :: return_time
!HPF$ DISTRIBUTE *(block) :: initial_time
!HPF$ ALIGN return_time with initial_time
      end subroutine timer

      END INTERFACE

      
      forall(i=1:nn) xx(i)=i
      forall(i=1:nn) yy(i)=i
      forall(i=1:nn) zz(i)=i

      order = 2
      dim = 3
      refinerat = 2
         
c      do k=1,13
c       do i=1,13
c        do j=1,13,3
c      child(i,j,k) = xx(i)*xx(i)*yy(j)*zz(k)
c         enddo
c        enddo
c      enddo

      forall(i=1:nn,j=1:nn,k=1:nn:2) child(i,j,k)=
     &                            xx(i)*xx(i)*yy(j)*zz(k)
       
      
      write(6,*) 'calling interpolate'

      call timer(t2,t1)

      write(6,*) 'within interpolate'


      step = refinerat
      no_of_planes=size(child,dim=3)
      no_of_rows=size(child,dim=2)
      no_of_cols=size(child,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) child(j+1,i,k) = 
     &                                    child(j,i,k)

      forall(i=1:no_of_rows:step) den(i)=(
     &                     child(j+step,i,k)-
     &                     child(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) child(j+1,i,k) = 
     &                                    child(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) child(i,j+1,k) = 
     &                                    child(i,j,k)
      forall(i=1:no_of_cols) den(i)=(child(i,j+step,k)-
     &                             child(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) child(i,j+1,k) = 
     &                                    child(i,j+1,k) + dy1(i)
c!      write(6,*) 'after 3rd forall'



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

c!      write(6,*) child

     
      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) child(i,j,k+1) = 
     &                                    child(i,j,k)
!      write(6,*) k,k+step
      forall(j=1:no_of_cols) den(j)=(child(i,j,k+step)-
     &                             child(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) child(i,j,k+1) = 
     &                                    child(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


      call timer(t1,t2)
      print *,'time',t1,' seconds'
  
 
      end program
      

      extrinsic(hpf_local) subroutine timer(return_time, initial_time)
      real, dimension(:), intent(in)  :: initial_time
      real, dimension(:), intent(out) :: return_time

      return_time(lbound(return_time, 1)) =     
     &        secnds(initial_time(lbound(initial_time, 1)))

      end



















