	PROGRAM polint

        integer, parameter :: ngx=2000,ngy=200  ! size of the grid
        integer, parameter :: n=3     ! order of interpolation
        integer, parameter :: maxlx=ngx-2*(n-1)
        integer, parameter :: maxly=ngy-2*(n-1)


!HPF$ TEMPLATE TMPL(ngx)
!HPF$ DISTRIBUTE TMPL(BLOCK)

        REAL, DIMENSION(ngx,ngy) :: grid
!HPF$ ALIGN grid(:,*) WITH TMPL(:) 

        REAL, DIMENSION(ngx)     :: coordx 
!HPF$ ALIGN (:) WITH TMPL(:)         :: coordx



	REAL, DIMENSION(ngx,n) :: cx
	INTEGER :: m,i,k,nmax

!HPF$ ALIGN (:,*) WITH TMPL(:) :: cx


!       initialize coordinates and function values 
        forall(i=1:ngx) coordx(i)=i-1
        forall(i=1:ngx:2,j=1:ngy:2) &
           grid(i,j)=3*coordx(i)+1  !initalize every second
!                                          the rest is to be interpolated

!write(6,*) 'input'
!
!write(6,*) grid(5:8,5)
!write(6,*) grid(5:8,6)
!write(6,*) grid(5:8,7)
!write(6,*) grid(5:8,8)

!t1=secnds(0.0)

! =========================================================================
do j=1,ngy,2

!for each row interpolate columnwise

forall(i=1:ngx:2)   cx(i,1)=grid(i,j)
!forall(i=1:ngx-2:2) cx(i,2)=grid(i+2,j)
!forall(i=1:ngx-4:2) cx(i,3)=grid(i+4,j)
enddo !rows

	END program polint

