module hpfprocessors

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

end module hpfprocessors

module gridfunctions
use hpfprocessors
real, dimension(2),parameter :: dr =(/ 1.0,1.0 /)
real, parameter              :: dt = 1.0

!grid functions            
integer, parameter           :: gnx=100,gny=100,gnz=1
integer, parameter           :: maxlev=2
integer, parameter           :: maxtimlevels=2
real, dimension(maxlev,maxtimlevels,gnx,gny,gnz)  :: g

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

integer, dimension(maxlev) :: nnx,nny,nnz

end module gridfunctions

program driver
use gridfunctions

integer icstep,ncsteps
real t

REAL time1(2), time2(2)

CALL dtime(time1)
ncsteps=10
g=0
t=0
nnx(1)=gnx
nny(1)=gny
nnz(1)=gnz

if(maxlev > 1) then
! pseduo-regridding
nnx(2)=100
nny(2)=100
nnz(2)=1
endif

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

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

stop
end

subroutine advance(t)
use gridfunctions
integer istep,i2step
integer, parameter :: nsteps=2
integer levels,curlev

do istep=1,nsteps
curlev=1
call update(curlev,nnx(curlev),nny(curlev),nnz(curlev))
call boundaries(t,curlev)
t=t+dt
write(6,*) t,'::',g(1,2,5,1:3,1), '*'
if(maxlev > 1) then
  curlev=2
  do i2step=1,nsteps
  call update(curlev,nnx(curlev),nny(curlev),nnz(curlev))
!  call boundaries(t,curlev)
  enddo
  call inject(curlev)
endif
enddo

return
end

subroutine update(curlev,nx,ny,nz)
use gridfunctions
integer :: curlev,nx,ny,nz

real, parameter :: cx=1.0,cy=1.0
real  sx,sy,sx2,sy2,sxy,sx1,sy1,sxy1

 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) g(curlev,1,i,j,1)=exp(sin(g(curlev,1,i,j,1)**2))
FORALL(i=1:nx,j=1:ny) g(curlev,1,i,j,1)=g(curlev,1,i,j,1)**3
enddo

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

g(curlev,2,:,:,:)=g(curlev,1,:,:,:)

return
end


subroutine boundaries(t,curlev)
use gridfunctions
integer curlev

!---------------------------------------------------
! set boundary conditions after grid update
!---------------------------------------------------

real                   :: t

g(curlev,2,1,:,1)=g(curlev,2,2,:,1)
g(curlev,2,nnx(curlev),:,1)=g(curlev,2,2,:,1)

! the following condition is imposed only when the edge of the
! grid g on a given level of refinement coincide with 
! edges of the base grid

if(curlev == 1) g(curlev,2,:,1,1)=sin(0.15*t)
return
end

subroutine inject(curlev)
integer curlev

return
end



