module hpfprocessors

integer, parameter :: NPX=2,NPY=2
!HPF$ PROCESSORS PROC(NPX,NPY)

end module hpfprocessors

module gridfunctions
use hpfprocessors
integer, parameter :: nx=100,ny=100,nz=100

!base grid
real, dimension(100,100,100), save :: g1

!HPF$ DISTRIBUTE g1 (BLOCK,BLOCK,*)

end module gridfunctions

program driver
use gridfunctions

integer icstep,ncsteps
REAL time1(2), time2(2)

CALL dtime(time1)
ncsteps=10
g1=0
g1(:,1,1)=5.0

do icstep=1,ncsteps
write(6,*) 'advance',icstep
call advance
enddo

CALL dtime(time2)
PRINT *, 'Time taken = ', time2(1) - time1(1)

stop
end

subroutine advance
use gridfunctions
integer istep
integer, parameter :: nsteps=2

do istep=1,nsteps
call update
enddo

return
end

subroutine update
use gridfunctions


!!!!!  real, dimension(:,:,:)   :: anew,aold
real                     :: t,dt
real, dimension(2)       :: dr
real, parameter :: cx=1.0,cy=1.0
real  sx,sy,sx2,sy2,sxy,sx1,sy1,sxy1

dr=1.0
dt=1.0

 sx=cx*dt/dr(1)
 sy=cy*dt/dr(2)
 sx2=sx**2
 sy2=sy**2
 sxy=sx*sy
 sx1=1.0/(1.0+sx)
 sy1=1.0/(1.0+sy)
 sxy1=1.0/(1.0+sx+sy)


do i=1,100
FORALL(i=1:nx,j=1:ny) g1(i,j,1)=exp(sin(g1(i,j,1)**2))
FORALL(i=1:nx,j=1:ny) g1(i,j,1)=g1(i,j,1)**3
enddo

Forall(i=2:nx-1,j=2:ny-1) g1(i,j,1)=                     &
       g1(i,j,1)                                         &
       -(g1(i+1,j,1)-g1(i-1,j,1))*sx/2.0-(g1(i,j+1,1)    &
       -g1(i,j-1,1))*sy/2.0                              &
       +(g1(i+1,j,1)-2.0*g1(i,j,1)+g1(i-1,j,1))*sx2/2.0  &
       +(g1(i,j+1,1)-2.0*g1(i,j,1)+g1(i,j-1,1))*sy2/2.0  &
       +(g1(i+1,j+1,1)-g1(i+1,j-1,1)-g1(i-1,j+1,1)       &
       +g1(i-1,j-1,1))*sxy/4.0 

!write(6,*) g1(5,1:10,1), '*'

return
end









