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
integer, parameter           :: rho = 2 ! regridding factor

!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$ TEMPLATE TMPL(gnx,gny)
!HPF$ DISTRIBUTE TMPL(BLOCK,BLOCK)
!HPF$ ALIGN g (*,*,:,:,*) WITH TMPL(:,:)

integer, dimension(maxlev,3) :: nr

!coordinates
real, dimension(maxlev,3) :: rmin,rmax
real, dimension(maxlev,gnx) :: coordx
real, dimension(maxlev,gny) :: coordy
real, dimension(maxlev,gnz) :: coordz
!HPF$ ALIGN coordx (*,:) WITH TMPL(:,*)  
!HPF$ ALIGN coordy (*,:) WITH TMPL(*,:)
! z coordinate is to be replicated

end module gridfunctions

program driver
use gridfunctions

integer icstep,ncsteps
real t

REAL time1(2), time2(2)

CALL dtime(time1)
ncsteps=10
t=0
call init_base_grid

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 init_base_grid
use gridfunctions

!number of grid poins
nr(1,1)=gnx
nr(1,2)=gny
nr(1,3)=gnz

!coordinates
rmin(1,:)=1
rmax(1,:)=nr(1,:)
FORALL(i=1:nr(1,1)) coordx(1,i)=real(i)
FORALL(i=1:nr(1,2)) coordy(1,i)=real(i)
FORALL(i=1:nr(1,3)) coordz(1,i)=real(i)

!initial valuse on the grid function
g(1,:,:,:,:)=0.0
return
end

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

do istep=1,nsteps
call regrid(t)
curlev=1
call update(curlev,nr)
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)
!  call boundaries(t,curlev)
  enddo
  call inject(curlev)
endif
enddo

return
end

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

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

nx=nr(curlev,1)
ny=nr(curlev,2)
nz=nr(curlev,3)

 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
real t
integer curlev

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

g(curlev,2,1,:,1)=g(curlev,2,2,:,1)
g(curlev,2,nr(curlev,1),:,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

subroutine init_fine_grid(level)
use gridfunctions
integer level
integer offsetx,offsety,offsetz
integer nx,ny,nz    !number of grid points on this level
integer nxc,nyc,nzc !number of gird points on parent level


!number of grid poins
nr(level,:)=(rmax(level,:)-rmin(level,:))*rho+1

!the following is for clarity only

 nx=nr(level,1)
 ny=nr(level,2)
 nz=nr(level,3)

 nxc=nr(level-1,1)
 nyc=nr(level-1,2)
 nzc=nr(level-1,3)

!coordinates

FORALL(i=1:nx) coordx(level,i)=rmin(level,1)+real(i-1)/rho
FORALL(i=1:ny) coordy(level,i)=rmin(level,2)+real(i-1)/rho
FORALL(i=1:nz) coordz(level,i)=rmin(level,3)+real(i-1)/rho

!initial values on the grid function

!compute offsets

  offsetx=(rmin(level,1)-rmin(level-1,1))+1
  offsety=(rmin(level,2)-rmin(level-1,2))+1
  offsetz=(rmin(level,3)-rmin(level-1,3))+1

!transfer data points from coarse grid to fine one

  g(level,:,1:nx:rho,1:ny:rho,1:nz:rho)=  &
  g(level-1,:,offsetx:nxc,offsety:nyc,offsetz:nyz) 

!now interpolate

  if(rho == 2) then
    if(mod(nx,2) <> 1 .or. mod(ny,2) <> 1 .or. mod(nz,2) <> 1) then
     write(6,*) 'it does not work for even size'
     stop
    endif

!    interpolate each k plane
  do k=1,nz,2
    forall(i=2:nx-1:2,j=1:ny:2)  &
     g(level,:,i,j,k)=0.5*(g(level,:,i-1,j,k)+g(level,:,i+1,j,k))
    
    forall(i=1:nx,j=2:ny:2) &
     g(level,:,i,j,k)=0.5*(g(level,:,i,j-1,k)+g(level,:,i,j+1,k))
  enddo

!   interpolate between k planes
     if(nz > 1) then
      forall(i=1:nx,j=1:ny,k=2:nz:2) &
      g(level,:,i,j,k)=0.5*(g(level,:,i,j,k-1)+g(level,:,i,j,k+1)) 
      endif
  else
   write(6,*) 'must supply interpolation routine'
   stop 
  endif
return
end

subroutine regrid(t)
use gridfunctions
real t
integer level

!very simple stuff; kind of placebo
!regrids up to maxlev

if(maxlev <> 2) then
   write(6,*) 'can do only maxlev=2'
   stop
endif

level=2
rmin(level,1)=26
rmin(level,2)=t+1.0
rmin(level,3)=1.0
rmax(level,1)=rmin(level,1)+real(gnx)/rho-1.0
rmax(level,2)=t+1.0+real(gny)/rho-1.0
if(rmax(level,2) > rmax(level-1,2)) rmax(level,2)=rmax(level-1,2)
rmax(level,3)=1.0
write(6,*) 'RMIN',rmin(level,:)
write(6,*) 'RMAX',rmax(level,:)
call init_fine_grid(level)

return
end


