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=1
integer, parameter           :: maxtimlevels=2
integer, parameter           :: maxgrids = 4
real, dimension(maxlev,maxtimlevels,maxgrids,gnx,gny,gnz)  :: g
real, dimension(maxlev,1,maxgrids,gnx,gny,gnz)  :: ctrest

!HPF$ TEMPLATE TMPL(gnx,gny)
!HPF$ DISTRIBUTE TMPL(BLOCK,BLOCK)
!HPF$ ALIGN g (*,*,*,:,:,*) WITH TMPL(:,:)
!HPF$ ALIGN ctrest (*,*,*,:,:,*) WITH TMPL(:,:)

integer, dimension(maxlev,maxgrids,3) :: nr
integer, dimension(maxgrids)          :: ngrids

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

!regridding
real, dimension(maxlev) :: threshold

end module gridfunctions

program driver
use gridfunctions

integer icstep,ncsteps,finlev
real t

REAL time1(2), time2(2)

CALL dtime(time1)
ncsteps=36
t=0
g=0
ngrids=0

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

call init_base_grid
finlev=1

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

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

stop
end

subroutine init_base_grid
use gridfunctions

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

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

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

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

call regrid(t,finlev)

do istep=1,nsteps

curlev=1
!write(6,*) 'update level',curlev
call update(curlev)
t=t+dt
call boundaries(t,curlev)
if(finlev > 1) then
  curlev=2
  t2=t-dt
  do ig=1,ngrids(2)
  do i2step=1,nsteps
write(6,*) 'update level',curlev
  call update(curlev,ig)
  t2=t2+dt/rho
  call boundaries(t2,curlev,ig) 
  enddo

  call inject(curlev,ig)
  enddo
endif

enddo

return
end

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

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)


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

do isteps=1,2
FORALL(i=2:nx-1,j=2:ny-1) &
                  g(curlev,1,ig,i,j,1)=exp(sin(g(curlev,1,ig,i,j,1)**2))
FORALL(i=2:nx-1,j=2:ny-1) g(curlev,1,ig,i,j,1)=g(curlev,1,ig,i,j,1)**3
enddo

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

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

return
end


subroutine boundaries(t,curlev,ig)
use gridfunctions
real t
integer curlev,ig

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

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


if(curlev==1) then 
g(curlev,2,ig,:,1,1)=sin(0.15*t)
else

if(rmin(curlev,ig,2)<1.5) then
! 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
  g(curlev,2,ig,:,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!!!'
stop
!  g(curlev,2,:,1,1)=0.5*(g(curlev-1,1,:,:,:)+g(curlev-1,2,:,:,:))
endif

endif

return
end

subroutine inject(level,ig)
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
integer ig


!the following is for clarity only

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


!compute offsets

  offsetx=(rmin(level,ig,1)-rmin(level-1,ig,1))+1
  offsety=(rmin(level,ig,2)-rmin(level-1,ig,2))+1
  offsetz=(rmin(level,ig,3)-rmin(level-1,ig,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

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

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
integer ig

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

!the following is for clarity only

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


!coordinates

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

!initial values on the grid function

!compute offsets

  offsetx=(rmin(level,ig,1)-rmin(level-1,ig,1))+1
  offsety=(rmin(level,ig,2)-rmin(level-1,ig,2))+1
  offsetz=(rmin(level,ig,3)-rmin(level-1,ig,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

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

!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,:,ig,i,j,k)=0.5*(g(level,:,ig,i-1,j,k)+g(level,:,ig,i+1,j,k))
    
    forall(i=1:nx,j=2:ny:2) &
     g(level,:,ig,i,j,k)=0.5*(g(level,:,ig,i,j-1,k)+g(level,:,ig,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,:,ig,i,j,k)=0.5*(g(level,:,ig,i,j,k-1)+g(level,:,ig,i,j,k+1)) 
      endif
  else
   write(6,*) 'must supply interpolation routine'
   stop 
  endif

return
end

subroutine regrid(t,finlev)
use gridfunctions
real t
integer level,finlev
logical, dimension(gny) :: ledge_y
integer, dimension(gny) :: edge_y
integer, dimension(1) :: id
real, dimension(2)    :: tmin,tmax
integer ig

!do i=1,10
!write(6,*) g(1,2,:,i,1),'*'
!enddo
!write(6,*) '++++++++++'

level=1
ig=1
ctrest(level,1,ig,:,:,:) = 0
where(g(level,2,ig,:,:,:) > threshold(level)) ctrest(level,1,ig,:,:,:) = 1

ledge_y=any(g(level,2,ig,:,:,1) > 0, dim=1)
edge_y=0.5
where(ledge_y) edge_y=1
forall(i=1:gny-2) edge_y(i)=edge_y(i)-edge_y(i+1)

tmin=1000.0
tmax=-1000.0
l=1
index=1000
do while(index>1)
id=minloc(edge_y)
index=id(1)
tmin(l)=real(index+1)
edge_y(index)=0
id=maxloc(edge_y)
index=id(1)
tmax(l)=real(index)
edge_y(index)=0
if(index>1) then
 if(l>2) then
    write(6,*) 'Oops! more than 2 sibling grids!'
  stop
 endif
 if(l==1) then
  rmin(level+1,ig,1)=tmin(1)
  rmax(level+1,ig,1)=tmax(1)
 endif
 if(l==2) then
   if(tmin(2)<tmin(1)) then
      rmin(level+1,ig,1)=tmin(2)
      rmin(level+1,ig,2)=tmin(1)
   else
      rmin(level+1,ig,1)=tmin(1)
      rmin(level+1,ig,2)=tmin(2)
   endif
   if(tmax(2)<tmax(1)) then
      rmax(level+1,ig,1)=tmax(2)
      rmax(level+1,ig,2)=tmax(1)
   else
      rmax(level+1,ig,1)=tmax(1)
      rmax(level+1,ig,2)=tmax(2)
   endif
  endif
else
  exit
endif 
l=l+1
enddo

l=l-1


if(l==0) then
  write(6,*) 'nothing to regrid'
elseif(l==1) then
  write(6,*) 'min ',rmin(level+1,ig,1)
  write(6,*) 'max ',rmax(level+1,ig,1)
else
  write(6,*) 'min ',rmin(level+1,ig,1:2)
  write(6,*) 'max ',rmax(level+1,ig,1:2)
endif

return
end

!!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
!g(level,:,:,:,:)=0
!rmin(level,1)=26
!rmax(level,1)=rmin(level,1)+real(gnx)/rho-1.0


!if(t<49.1) then
! rmin(level,2)=1.0
!else
! rmin(level,2)=t-49.0
!endif
!rmax(level,2)=rmin(level,2)+real(gny)/rho-1.0
!if(rmax(level,2) > rmax(level-1,2)) rmax(level,2)=rmax(level-1,2)

!rmin(level,3)=1.0
!rmax(level,3)=1.0

!write(6,*) 'RMIN',rmin(level,:)
!write(6,*) 'RMAX',rmax(level,:)
!call init_fine_grid(level)

!write(6,*) t,'+',g(2,2,1:3,1,1), '*'
!write(6,*) t,'+',g(2,2,1:3,2,1), '*'
!write(6,*) t,'+',g(2,2,1:3,3,1), '*'

