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=400,gny=400,gnz=1
integer, parameter           :: maxlev=1
integer, parameter           :: mtl=2 !max time levels
integer, parameter           :: mgd=3 !max grids on the same level

real, dimension(mtl,gnx,gny,gnz),target     :: g1 !base grid
real, dimension(:,:,:,:),allocatable,target :: g21,g22,g23 !3 grids at level 2
real, dimension(:,:,:,:),allocatable,target :: g31,g32,g33 !3 grids at level 3
real, dimension(:,:,:,:),allocatable        :: ctrest


!!  HPF$ TEMPLATE TMPL(gnx,gny)
!!  HPF$ DISTRIBUTE TMPL(BLOCK,BLOCK) onto PROC
!!  HPF$ ALIGN (*,:,:,*) WITH TMPL(:,:) :: g1,g21,g22,g23,g31,g32,g33
!!  HPF$ ALIGN ctrest (*,:,:,*) WITH TMPL(:,:)
!HPF$ DISTRIBUTE (*,BLOCK,BLOCK,*) onto PROC :: g1,g21,g22,g23,g31,g32,g33
!HPF$ DISTRIBUTE ctrest (*,BLOCK,BLOCK,*) onto PROC

real, dimension(:,:,:,:), pointer :: p,pp,pc

!!coordinates
real, dimension(maxlev,mgd,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

!regridding
real, dimension(maxlev) :: threshold

type ADAPTIVE_MESH
 integer                        :: grid_id
 integer, dimension(3)          :: nr
 real, dimension(3)             :: rmin
 real, dimension(3)             :: rmax
 integer, dimension(3,2)        :: ibound !3d,min,max =1 -physical boundaries
 type (ADAPTIVE_MESH), pointer  :: parent
 type (ADAPTIVE_MESH), pointer  :: neighbor
 type (ADAPTIVE_MESH), pointer  :: child  
end type ADAPTIVE_MESH

type LEVEL_INFO
 integer                        :: ngd
 type (ADAPTIVE_MESH), pointer  :: grid 
end type LEVEL_INFO

type (LEVEL_INFO), dimension(maxlev) :: lp

end module gridfunctions

program driver
use gridfunctions
   
integer icstep,ncsteps,finlev
real t

REAL time1(2), time2(2)

CALL dtime(time1)
ncsteps=45
t=0
g=0

do i=1,maxlev
threshold(i)=0.08*i
enddo

do i=1,maxlev
 lp(i)%ngd=0
 nullify(lp(i)%grid)
enddo

call init_base_grid
finlev=1

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

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

stop
end

subroutine init_base_grid
use gridfunctions

type (ADAPTIVE_MESH), pointer  :: basegrid

allocate(basegrid)
basegrid%grid_id=1
basegrid%nr(1)=gnx
basegrid%nr(2)=gny
basegrid%nr(3)=gnz
basegrid%rmin=1.0
basegrid%rmax=basegrid%nr
basegrid%ibound=1
nullify(basegrid%parent)
nullify(basegrid%child)
nullify(basegrid%neighbor)

lp(1)%ngd=1
lp(1)%grid=>basegrid


!coordinates
!rmin(1,1,:)=1
!rmax(1,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 values on the grid function
g1(:,:,:,:)=0.0
return
end

subroutine advance(t,finlev,icstep)
use gridfunctions
real t
integer istep,i2step,icstep
integer, parameter :: nsteps=2
integer levels,curlev,finlev
integer ipoint
real t2
type(ADAPTIVE_MESH),pointer :: grid

if(mod(icstep,2) == 1) call regrid(t,finlev)

do istep=1,nsteps

curlev=1
p=>g1
nx=lp(curlev)%grid%nr(1)
ny=lp(curlev)%grid%nr(2)
nz=lp(curlev)%grid%nr(3)
!write(6,*) 'update level',curlev,nx,ny,nz
call update(curlev,nx,ny,nz)
t=t+dt
id=1
call boundaries(t,curlev,nx,lp(1)%grid%ibound,id)
pp=>g1
if(finlev > 1) then
  curlev=2
  t2=t-dt
  do i2step=1,nsteps
   t2=t2+dt/rho
   if(allocated(g21)) then
    p=>g21
    nx=lp(curlev)%grid%nr(1)
    ny=lp(curlev)%grid%nr(2)
    nz=lp(curlev)%grid%nr(3)
    call update(curlev,nx,ny,nz)
    id=1
    call boundaries(t2,curlev,nx,lp(2)%grid%ibound,id) 
   endif
   if(allocated(g22)) then
     p=>g22
    grid=>lp(curlev)%grid
    grid=>grid%neighbor 
    nx=grid%nr(1)
    ny=grid%nr(2)
    nz=grid%nr(3)
    call update(curlev,nx,ny,nz)
    id=2
    call boundaries(t2,curlev,nx,grid%ibound,id) 
   endif
   if(allocated(g23)) then
     p=>g23
    grid=>lp(curlev)%grid
    grid=>grid%neighbor 
    nx=grid%nr(1)
    ny=grid%nr(2)
    nz=grid%nr(3)
    call update(curlev,nx,ny,nz)
    id=3
    call boundaries(t2,curlev,nx,grid%ibound,id) 
   endif

  if(allocated(g21)) g21(2,:,:,:)=g21(1,:,:,:)
  if(allocated(g22)) g22(2,:,:,:)=g22(1,:,:,:)
  if(allocated(g23)) g23(2,:,:,:)=g23(1,:,:,:)
  enddo

   if(allocated(g21)) then
    p=>g21
    pp=>g1
    id=1
    call inject(curlev,id)
    endif

   if(allocated(g22)) then
    p=>g22
    pp=>g1
    id=2
    call inject(curlev,id)
    endif

   if(allocated(g23)) then
    p=>g23
    pp=>g1
    id=3
    call inject(curlev,id)
    endif

endif

g1(2,:,:,:)=g1(1,:,:,:)
!call dumpdata(finlev)
write(6,*) g1(2,1,1:4,1)

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

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

!p(2,:,:,:)=p(1,:,:,:)

return
end


subroutine boundaries(t,curlev,nx,ibound,id)
use gridfunctions
real t
integer curlev,nx,id
integer, dimension(3,2) :: ibound
type(ADAPTIVE_MESH),pointer :: grid

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

! if ibound==1 then physical boundaries

if(ibound(1,1)==1) then
  p(1,1,:,1)=p(1,2,:,1)
else
 write(6,*) 'interpolation needed (1)'
 stop
endif
 
if(ibound(1,2)==1) then
  p(1,nx,:,1)=p(1,2,:,1)
else
 write(6,*) 'interpolation needed (2)'
 stop
endif

if(ibound(2,1)==1) then 
p(1,:,1,1)=sin(0.15*t)

else
!here I assumed that there are two steps on fine grid;
!first step : boundary conditions are just interpolated when regridding
!             that is, this routine is called after first update
!second step: start with values that are average between retarded
!             and advanced solution on level-1
!             this is first call
!second call: who cares, the result will be overwritten by regridding
!             anyway. What about transfer????    
!
!write(6,*) 'insert interpolation here!!!'

!first transfer averaged data from parent to child

!compute offsets

  grid=>lp(curlev)%grid
  if(id <> 1) then
    do while (grid%grid_id <> id)
     grid=>grid%neighbor
    enddo
  endif

  nnx=grid%nr(1)
  ny=grid%nr(2)
  nz=grid%nr(3)

if(nnx <> nx) then
  write(6,*) 'BOUNDARIES: something wrong with nx'
  stop
endif
  
  offsetx=(grid%rmin(1)-grid%parent%rmin(1))+1
  offsety=(grid%rmin(2)-grid%parent%rmin(2))+1
  offsetz=(grid%rmin(3)-grid%parent%rmin(3))+1

  nxc=offsetx+(nx+1)/rho-1
  nyc=offsety+(ny+1)/rho-1
  nzc=offsetz+(nz+1)/rho-1

!transfer data points from coarse grid to fine one
!for all x and z, and only y=1 (offsety)

!write(6,*) p(1,5,1,1),'*'
!write(6,*) pp(1,5,offsety,1),'*'
!write(6,*) pp(2,5,offsety,1),'*'


  p(1,1:nx:rho,1,1:nz:rho)=  &
  0.5*( pp(1,offsetx:nxc,offsety,offsetz:nzc)+ & 
        pp(2,offsetx:nxc,offsety,offsetz:nzc))

!write(6,*) p(1,5,1,1),'**'

!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 (along x axis only)
  do k=1,nz,2
    forall(i=2:nx-1:2,j=1:ny:2)  &
     p(1,i,j,k)=0.5*(p(1,i-1,j,k)+p(1,i+1,j,k))
    
!    forall(i=1:nx,j=2:ny:2) &
!     p(:,i,j,k)=0.5*(p(:,i,j-1,k)+p(:,i,j+1,k))
  enddo

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

endif

!do not care about other boundaries
return
end

subroutine inject(level,id)
use gridfunctions
integer level,id
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
type(ADAPTIVE_MESH), pointer :: grid

  grid=>lp(level)%grid
  if(id <> 1) then
    do while (grid%grid_id <> id)
     grid=>grid%neighbor
    enddo
  endif

  nx=grid%nr(1)
  ny=grid%nr(2)
  nz=grid%nr(3)
  
  offsetx=(grid%rmin(1)-grid%parent%rmin(1))+1
  offsety=(grid%rmin(2)-grid%parent%rmin(2))+1
  offsetz=(grid%rmin(3)-grid%parent%rmin(3))+1

  nxc=offsetx+(nx+1)/rho-1
  nyc=offsety+(ny+1)/rho-1 
  nzc=offsetz+(nz+1)/rho-1

!transfer data points from fine grid to coarse one

 pp(:,offsetx:nxc,offsety:nyc,offsetz:nzc)= & 
 p(:,1:nx:rho,1:ny:rho,1:nz:rho)

return
end










