	PROGRAM polint
	USE types

        integer, parameter :: ngx=200,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)


        REAL(SP), DIMENSION(ngx,ngy) :: grid
!HPF$ DISTRIBUTE grid(BLOCK,BLOCK)

        REAL(SP), DIMENSION(ngx)     :: coordx 
        REAL(SP), DIMENSION(ngy)     :: coordy 
!HPF$ ALIGN (:) WITH grid(:,*) :: coordx
!HPF$ ALIGN (:) WITH grid(*,:) :: coordy


        REAL(SP), DIMENSION(ngy) :: y
!HPF$ ALIGN (:) WITH grid(*,:) :: y

        REAL(SP), DIMENSION(ngx) :: dgx
        REAL(SP), DIMENSION(ngy) :: dgy
        REAL(SP), DIMENSION(ngx,n) :: xa
	REAL(SP), DIMENSION(ngx,n) :: cx,dx,hox,denx
        REAL(SP), DIMENSION(ngx,n) :: ya,gya
	REAL(SP), DIMENSION(ngx,n) :: cy,dy,hoy,deny
	INTEGER(I4B) :: m,i,k,nmax

!  !HPF$ PROCESSORS PROC(4)
!  !HPF$ DISTRIBUTE grid(BLOCK) ONTO PROC


!HPF$ ALIGN (:)   WITH grid(:,*) :: dgx
!HPF$ ALIGN (:)   WITH grid(*,:) :: dgy
!HPF$ ALIGN (:,*) WITH grid(:,*) :: xa,cx,dx,hox,denx
!HPF$ ALIGN (:,*) WITH grid(*,:) :: ya,gya,cy,dy,hoy,deny
  

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

!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)


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

!for each row interpolate columnwise

do k=1,n
nmax=ngx-2*(k-1)
forall(i=1:ngx:2) cx(i,k)=grid(i+2*(k-1),j)
forall(i=1:ngx:2) dx(i,k)=cx(i,k)
forall(i=1:nmax:2) hox(i,k)=coordx(i+2*(k-1))-coordx(i+1)
enddo

forall(i=1:ngx-2*(n-1):2) grid(i+1,j)=grid(i,j)


do m=1,n-1
do k=1,n-m   
forall(i=1:maxlx:2) denx(i,k)=hox(i,k)-hox(i,k+m)
forall(i=1:maxlx:2) denx(i,k)=(cx(i,k+1)-dx(i,k))/denx(i,k)
forall(i=1:maxlx:2) dx(i,k)=hox(i,k+m)*denx(i,k)
forall(i=1:maxlx:2) cx(i,k)=hox(i,k)*denx(i,k)
enddo
!forall(i=1:maxlx:2) dgx(i)=cx(i,ns(i)+1)
forall(i=1:maxlx:2) dgx(i)=cx(i,1)
forall(i=1:maxlx:2) grid(i+1,j)=grid(i+1,j)+dgx(i)
enddo

enddo !rows

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


! ===========================================================================
forall(j=1:ngy-n:2) y(j)=coordy(j+1)

do i=1,ngx
!do i=1,3

!for each colume interpolate rowwise

do k=1,n
nmax=ngy-2*(k-1)
forall(j=1:nmax:2) ya(j,k)=coordy(j+2*(k-1))
forall(j=1:nmax:2) gya(j,k)=grid(i,j+2*(k-1))
forall(j=1:ngy:2) cy(j,k)=gya(j,k)
forall(j=1:ngy:2) dy(j,k)=gya(j,k)
forall(j=1:nmax:2) hoy(j,k)=ya(j,k)-y(j)
enddo

forall(j=1:ngy-2*(n-1):2) grid(i,j+1)=gya(j,1)

do m=1,n-1
do k=1,n-m   
forall(j=1:maxly:2) deny(j,k)=hoy(j,k)-hoy(j,k+m)
forall(j=1:maxly:2) deny(j,k)=(cy(j,k+1)-dy(j,k))/deny(j,k)
forall(j=1:maxly:2) dy(j,k)=hoy(j,k+m)*deny(j,k)
forall(j=1:maxly:2) cy(j,k)=hoy(j,k)*deny(j,k)
enddo
!forall(i=1:maxlx:2) dgx(i)=cx(i,ns(i)+1)
forall(j=1:maxly:2) dgy(j)=cy(j,1)
forall(j=1:maxly:2) grid(i,j+1)=grid(i,j+1)+dgy(j)
enddo

enddo

write(6,*) 'after rows'

write(6,*) grid(5:8,5)
write(6,*) grid(5:8,6)
write(6,*) grid(5:8,7)
write(6,*) grid(5:8,8)



	END program polint

