module hpfprocessors

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

end module hpfprocessors

module gridfunctions
use hpfprocessors
real, parameter              :: dr_base = 1.0
integer, parameter           :: rho = 2 ! regridding factor
real, parameter              :: dtbydr = 1.0

!grid functions            
integer, parameter           :: gnx=100,gny=100,gnz=1
integer, parameter           :: maxlev=3
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
 real                           :: dr
 real                           :: dt
 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
!ncsteps=6
t=0
g=0

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

lp(1)%dr=dr_base
do i=2,maxlev
lp(i)%dr=lp(i-1)%dr/rho
enddo

do i=1,maxlev
 lp(i)%ngd=0
 lp(i)%dt=lp(i)%dr*dtbydr
 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 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*dtbydr
 sy=cy*dtbydr
 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,2
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










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

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+lp(curlev)%dt
id=1
call boundaries(t,curlev,nx,lp(1)%grid%ibound,id)
pp=>g1
if(finlev > 1) then
  t2=t-lp(curlev)%dt
  curlev=2
  do i2step=1,nsteps
   t2=t2+lp(curlev)%dt
   if(allocated(g21)) then
    grid=>lp(curlev)%grid 
    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)
    if(associated(grid%child)) then
       gridbck=>grid
       grid=>grid%child
       t3=t2-lp(curlev)%dt
       curlev=curlev+1      
       do i3step=1,nsteps
           t3=t3+lp(curlev)%dt
           if(allocated(g31)) then
              p=>g31
              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(t3,curlev,nx,lp(3)%grid%ibound,id)
           endif
       if(allocated(g31)) g31(2,:,:,:)=g31(1,:,:,:)
       enddo  ! i3step   
       call boundaries(t3,curlev,nx,grid%ibound,id) 
       grid=>gridbck
       curlev=curlev-1  
     endif  !associated(grid%child)
   endif  !g21 
   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) 
    if(associated(grid%child)) then
       gridbck=>grid
       grid=>grid%child
       t3=t2-lp(curlev)%dt
       curlev=curlev+1      
       do i3step=1,nsteps
           t3=t3+lp(curlev)%dt
           if(allocated(g32)) then
              p=>g32
              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(t3,curlev,nx,lp(3)%grid%ibound,id)
           endif
       if(allocated(g32)) g32(2,:,:,:)=g32(1,:,:,:)
       enddo  ! i3step   
       call boundaries(t3,curlev,nx,grid%ibound,id) 
       grid=>gridbck
       curlev=curlev-1  
     endif  !associated(grid%child)
   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) 
    if(associated(grid%child)) then
       gridbck=>grid
       grid=>grid%child
       t3=t2-lp(curlev)%dt
       curlev=curlev+1      
       do i3step=1,nsteps
           t3=t3+lp(curlev)%dt
           if(allocated(g33)) then
              p=>g33
              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(t3,curlev,nx,lp(3)%grid%ibound,id)
           endif
       if(allocated(g33)) g33(2,:,:,:)=g33(1,:,:,:)
       enddo  ! i3step   
       call boundaries(t3,curlev,nx,grid%ibound,id) 
       grid=>gridbck
       curlev=curlev-1  
     endif  !associated(grid%child)
   endif

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

   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 !finlev>1

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

enddo ! istep

return
end
subroutine dumpdata(finlev)
use gridfunctions

type(ADAPTIVE_MESH),pointer :: grid
integer gridno,finlev
integer lev
integer, SAVE :: timestep=0

gridno=0

!first the base grid

lev=1
grid=>lp(1)%grid
write(72,*) lev
write(72,*) grid%nr(1:2)
write(72,*) grid%rmin(1:2)
write(72,*) grid%rmax(1:2)

gridno=gridno+1
do i=1,grid%nr(1)
 do j=1,grid%nr(2)
   write(72,*) g1(1,i,j,1)
 enddo
enddo

!now level 2
if(finlev == 2) then
  lev=2
  grid=>lp(2)%grid
  write(72,*) lev
  write(72,*) grid%nr(1:2)
  write(72,*) grid%rmin(1:2)
  write(72,*) grid%rmax(1:2)

  gridno=gridno+1
  do i=1,grid%nr(1)
   do j=1,grid%nr(2)
    write(72,*) g21(1,i,j,1)
    enddo
  enddo

  if(allocated(g22)) then 
   grid=>grid%neighbor
   write(72,*) lev
   write(72,*) grid%nr(1:2)
   write(72,*) grid%rmin(1:2)
   write(72,*) grid%rmax(1:2)

   gridno=gridno+1
   do i=1,grid%nr(1)
   do j=1,grid%nr(2)
    write(72,*) g22(1,i,j,1)
    enddo
   enddo
  endif   
endif

write(73,*) timestep,gridno
timestep=timestep+1

return
end




subroutine find_minmax(tmin,tmax,lx,ly,lz,nr,level)
use gridfunctions
real, dimension(:,:)    :: tmin,tmax !3 grids; 3d
integer lx,ly,lz
integer, dimension(3) :: nr
integer level
integer, dimension(1) :: id
real, dimension(3)    :: gmin,gmax

logical, dimension(size(ctrest,dim=2)) :: ledge_x
logical, dimension(size(ctrest,dim=3)) :: ledge_y
logical, dimension(size(ctrest,dim=4)) :: ledge_z
integer, dimension(size(ctrest,dim=2)) :: edge_x
integer, dimension(size(ctrest,dim=3)) :: edge_y
integer, dimension(size(ctrest,dim=4)) :: edge_z

!! HPF$ TEMPLATE TMPL(size(ledge_x),size(ledge_y))
!! HPF$ DISTRIBUTE TMPL(BLOCK,BLOCK)
!! HPF$ ALIGN (:) WITH TMPL(:,*) :: ledge_x,edge_x
!! HPF$ ALIGN (:) WITH TMPL(*,:) :: ledge_y,edge_y

gmin(:)=tmin(1,:)
gmax(:)=tmax(1,:)


tmin=1.0
do i=1,size(tmax,dim=1)
tmax(:,i)=real(nr(i))
enddo

!find xmin and xmax


ledge_x=any(ctrest(1,:,:,1) > 0, dim=2)
edge_x=0.5
where(ledge_x) edge_x=1
forall(i=1:nr(1)-2) edge_x(i)=edge_x(i)-edge_x(i+1)

lx=1
index=1000

do while(index>1)
    id=minloc(edge_x)
    index=id(1)
    tmin(lx,1)=real(index+1) 
    edge_y(index)=0
    id=maxloc(edge_x)
    index=id(1)
    tmax(lx,1)=real(index)
    edge_x(index)=0

    if(index>1) then
       if(lx>3) then
         write(6,*) 'Oops! more than 3 sibling grids!'
         stop
       endif
      else
         exit
         write(6,*) 'exit'         
endif 
lx=lx+1
if(tmin(lx,1)<1.1 .AND. tmax(lx,1) > gmax(1)-0.1) index=0
enddo
lx=lx-1
if(lx==0) then
  tmin(1,1)=1.0
  tmax(1,1)=real(nr(1))
  lx=1
endif

! now ymin, ymax

ledge_y=any(ctrest(1,:,:,1) > 0, dim=1)
edge_y=0.5
where(ledge_y) edge_y=1
! we do not want maximum at i=1
forall(i=1:nr(2)-2) edge_y(i)=edge_y(i)-edge_y(i+1)

tmin=1.0
ly=1
index=1000

do while(index>1)
    id=minloc(edge_y)
    index=id(1)
    tmin(ly,2)=real(index+1) 
    edge_y(index)=0
    id=maxloc(edge_y)
    index=id(1)
    tmax(ly,2)=real(index)
    edge_y(index)=0
    if(ly==1 .AND. tmin(ly,2)>tmax(ly,2) .AND. index==1) index=1000

    if(index>1) then
       if(ly>3) then
         write(6,*) 'Oops! more than 3 sibling grids! ly=',ly
         stop
       endif
      else
         exit
endif 
ly=ly+1
enddo
ly=ly-1

!finally zmin and zmax
lz=1
tmin(:,3)=1.0
tmax(:,3)=1.0

!now sort things up

!anything to regrid?

if(lx < 1 .OR. ly < 1) return
if(lx > 3) then
  write(6,*) 'Oops! Cannot handle lx > 3:',lx
  stop
endif
if(ly > 3) then
  write(6,*) 'Oops! Cannot handle ly > 3:',ly
  stop
endif

if(lx == 3 .AND. tmin(1,1) > tmin(3,1)) then
   temp=tmin(1,1)
   tmin(1,1)=tmin(3,1)
   tmin(3,1)=temp
endif
if(lx == 3 .AND. tmax(1,1) > tmax(3,1)) then
   temp=tmax(1,1)
   tmax(1,1)=tmax(3,1)
   tmax(3,1)=temp
endif

if(ly == 3 .AND. tmin(1,2) > tmin(3,2)) then
   temp=tmin(1,2)
   tmin(1,2)=tmin(3,2)
   tmin(3,2)=temp
endif
if(ly == 3 .AND. tmax(1,2) > tmax(3,2)) then
   temp=tmax(1,2)
   tmax(1,2)=tmax(3,2)
   tmax(3,2)=temp
endif

if(lx > 1 .AND. tmin(1,1) > tmin(2,1)) then
   temp=tmin(1,1)
   tmin(1,1)=tmin(2,1)
   tmin(2,1)=temp
endif
if(lx > 1 .AND. tmax(1,1) > tmax(2,1)) then
   temp=tmax(1,1)
   tmax(1,1)=tmax(2,1)
   tmax(2,1)=temp
endif

if(ly > 1 .AND. tmin(1,2) > tmin(2,2)) then
   temp=tmin(1,2)
   tmin(1,2)=tmin(2,2)
   tmin(2,2)=temp
endif
if(ly > 1 .AND. tmax(1,2) > tmax(2,2)) then
   temp=tmax(1,2)
   tmax(1,2)=tmax(2,2)
   tmax(2,2)=temp
endif

if(lx == 3 .AND. tmin(2,1) > tmin(3,1)) then
   temp=tmin(2,1)
   tmin(2,1)=tmin(3,1)
   tmin(3,1)=temp
endif
if(lx == 3 .AND. tmax(2,1) > tmax(3,1)) then
   temp=tmax(2,1)
   tmax(2,1)=tmax(3,1)
   tmax(3,1)=temp
endif

if(ly == 3 .AND. tmin(2,2) > tmin(3,2)) then
   temp=tmin(2,2)
   tmin(2,2)=tmin(3,2)
   tmin(3,2)=temp
endif
if(ly == 3 .AND. tmax(2,2) > tmax(3,2)) then
   temp=tmax(2,2)
   tmax(2,2)=tmax(3,2)
   tmax(3,2)=temp
endif

! convert indices to coordinates

do i=1,lx
   tmin(i,1)=gmin(1)+tmin(i,1)*lp(level)%dr
   tmax(i,1)=gmin(1)+tmax(i,1)*lp(level)%dr
enddo 
do i=1,ly
   tmin(i,2)=gmin(2)+tmin(i,2)*lp(level)%dr
   tmax(i,2)=gmin(2)+tmax(i,2)*lp(level)%dr
!   tmax(i,2)=gmin(2)+tmax(i,2)
enddo 

!if rmin is too close to the physical edge, extend the grid
do i=1,lx
if(tmin(i,1) < gmin(1)+1.5*lp(level)%dr) tmin(i,1)=gmin(1)
enddo
do i=1,ly
if(tmin(i,2) < gmin(2)+1.5*lp(level)%dr) tmin(i,2)=gmin(2)
enddo

!if rmax is too close to the physical edge, extend the grid
do i=1,lx
if(tmax(i,1) > gmax(1)-1.5*lp(level)%dr) tmax(i,1)=gmax(1)
enddo
do i=1,ly
if(tmax(i,2) > gmax(2)-1.5*lp(level)%dr) tmax(i,2)=gmax(1)
enddo

!for transport, since I do not handle ymax boundary conditions
!proprely, extend ymax
!also make regions larger to decrease regridding frequency 

do i=1,ly
tmax(i,2)=tmax(i,2)+10.0
tmax(i,2)=min(tmax(i,2),gmax(2)-1)
enddo

!lp(level)%ngd=lx*ly
if(lx*ly > mgd) then
  write(6,*) 'Oops! maximum of',mgd, 'grids on level',level,' exceeded:',lx*ly
  stop
endif

if(lx*ly == 0) then
  write(6,*) 'Is it OK? No new grids on level', level
endif

return
end

subroutine allocate_fine_grid(level,minx,maxx,miny,maxy,minz,maxz,finlev)
use gridfunctions
real minx,maxx,miny,maxy,minz,maxz
integer finlev
type (ADAPTIVE_MESH),pointer :: grid,newgrid

n1=(maxx-minx)/lp(level)%dr+1
n2=(maxy-miny)/lp(level)%dr+1
n3=(maxz-minz)/lp(level)%dr+1


!reject too small clusters
if(abs(n1) < 5.0) return
if(abs(n2) < 5.0) return

!check nesting

!add buffers

allocate(newgrid) 
newgrid%nr(1)=n1
newgrid%nr(2)=n2
newgrid%nr(3)=n3
newgrid%rmin(1)=minx
newgrid%rmin(2)=miny
newgrid%rmin(3)=minz
newgrid%rmax(1)=maxx
newgrid%rmax(2)=maxy
newgrid%rmax(3)=maxz
nullify(newgrid%child)
nullify(newgrid%neighbor)
newgrid%parent=>lp(level-1)%grid

newgrid%ibound=0
if(abs(minx-newgrid%parent%rmin(1)) < 0.001) newgrid%ibound(1,1)=1
if(abs(miny-newgrid%parent%rmin(2)) < 0.001) newgrid%ibound(2,1)=1
if(abs(minz-newgrid%parent%rmin(3)) < 0.001) newgrid%ibound(3,1)=1
if(abs(maxx-newgrid%parent%rmax(1)) < 0.001) newgrid%ibound(1,2)=1
if(abs(maxy-newgrid%parent%rmax(2)) < 0.001) newgrid%ibound(2,2)=1
if(abs(maxz-newgrid%parent%rmax(3)) < 0.001) newgrid%ibound(3,2)=1



pp=>g1

if(lp(level)%ngd == 0) then
  lp(level)%grid=>newgrid
else
  grid=>lp(level)%grid
  do while(associated(grid%neighbor))
     grid=>grid%neighbor
  enddo
  grid%neighbor=>newgrid
endif
lp(level)%ngd=lp(level)%ngd+1
newgrid%grid_id=lp(level)%ngd

if(level==2) then
  if(newgrid%grid_id == 1) then
     if(allocated(g21)) deallocate(g21) 
     allocate(g21(mtl,n1,n2,n3))
     pc=>g21
     call init_fine_grid(level,newgrid%grid_id,finlev)
  endif
  if(newgrid%grid_id == 2) then 
     if(allocated(g22)) deallocate(g22) 
     allocate(g22(mtl,n1,n2,n3))
     pc=>g22
     call init_fine_grid(level,newgrid%grid_id,finlev)
   endif
   if(newgrid%grid_id == 3) then 
     if(allocated(g23)) deallocate(g23) 
     allocate(g23(mtl,n1,n2,n3))
     pc=>g23
     call init_fine_grid(level,newgrid%grid_id,finlev)
  endif
  finlev=2
endif  

return
end

subroutine deallocate_fine_grids
use gridfunctions
type (ADAPTIVE_MESH),pointer :: grid1,grid2

do i=2,maxlev
   if(lp(i)%ngd > 0) then
     grid1=>lp(i)%grid
     do j=1,lp(i)%ngd
        grid2=>grid1%neighbor 
        deallocate(grid1)
        grid1=>grid2
     enddo
   endif
   lp(i)%ngd=0
enddo

if(allocated(g21)) deallocate(g21)
if(allocated(g22)) deallocate(g22)
if(allocated(g23)) deallocate(g23)
if(allocated(g31)) deallocate(g31)
if(allocated(g32)) deallocate(g31)
if(allocated(g33)) deallocate(g31)

return
end

subroutine trunc_error(nr,level)
use gridfunctions
integer, dimension(3) :: nr

allocate(ctrest(1,nr(1),nr(2),nr(3)))

ctrest(1,:,:,:) = 0 
where(p(2,:,:,:) > threshold(level)) ctrest(1,:,:,:) = 1.0

return
endsubroutine init_fine_grid(level,id,finlev)
use gridfunctions
integer level,id,finlev
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

!pp is the parent grid
!pc is the child grid


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

!the following is for clarity only

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

!coordinates

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

!initial values on the grid function

!compute offsets
  
  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

  pc=0

  pc(:,1:nx:rho,1:ny:rho,1:nz:rho)=  &
  pp(:,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)  &
     pc(:,i,j,k)=0.5*(pc(:,i-1,j,k)+pc(:,i+1,j,k))
    
    forall(i=1:nx,j=2:ny:2) &
     pc(:,i,j,k)=0.5*(pc(:,i,j-1,k)+pc(:,i,j+1,k))
  enddo

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

finlev=level
return
end



















subroutine regrid(t,finlev)
use gridfunctions
real t
integer level,finlev
integer ngd,igrid

integer, dimension(1) :: id
real, dimension(3,3)    :: tmin,tmax !3 grids; 3d
integer lx,ly,lz
integer nx,ny,nz,finlev_old
type (ADAPTIVE_MESH), pointer :: grid

INTERFACE
subroutine find_minmax(tmin,tmax,lx,ly,lz,nr,level)
use gridfunctions
real, dimension(:,:)    :: tmin,tmax !3 grids; 3d
integer lx,ly,lz
integer, dimension(3) :: nr
integer level
end subroutine find_minmax
END INTERFACE

finlev_old=finlev
if(finlev_old > maxlev-1) finlev_old=maxlev-1

call deallocate_fine_grids

!always start from the base grid


do level=1,finlev_old

write(6,*) 'start level',level,' of',finlev_old

if(.not.associated(lp(level)%grid)) cycle
ngd=lp(level)%ngd

grid=>lp(level)%grid

if(ngd < 1) then
  write(6,*) 'find_minmax: something wrong, ngd=',ngd
  stop
endif

do igrid=1,ngd

SELECT CASE (level)

CASE (1)

      p=>g1

CASE (2)

      SELECT CASE (igrid)
      CASE (1)
        p=>g21 
      CASE (2)  
        p=>g22
      CASE (3)
        p=>g23
      CASE DEFAULT
        write(6,*) 'find_minmax: cannot handle level',level,' grid',igrid
      END SELECT

CASE (3)

      SELECT CASE (igrid)
      CASE (1)
        p=>g31 
      CASE (2)  
        p=>g32
      CASE (3)
        p=>g33
      CASE DEFAULT
        write(6,*) 'find_minmax: cannot handle level',level,' grid',igrid
      END SELECT

CASE (4)
      
      write(6,*) 'find_minmax: cannot handle level',level
      stop

END SELECT

if(allocated(ctrest)) deallocate(ctrest)
call trunc_error(grid%nr,level)

if(all(ctrest(1,:,:,:) < 0.5)) then
 write(6,*) 'all ctrest < 0.5'
 deallocate(ctrest)
 return
endif


do i=1,3
tmin(:,i)=grid%rmin(i)
tmax(:,i)=grid%rmax(i)
enddo


call find_minmax(tmin,tmax,lx,ly,lz,grid%nr,level)

!write(6,*) 'tmin'
!write(6,*) tmin
!write(6,*) 'tmax'
!write(6,*) tmax
!read(5,*) inp
!if(inp==0) stop

if(level < maxlev) then 

write(6,*) 'call allocate_fine_grid',lx,ly,lz
do i=1,lx
do j=1,ly
do k=1,lz
call allocate_fine_grid(level+1,tmin(i,1),tmax(i,1),tmin(j,2),tmax(j,2), &
                        tmin(k,3),tmax(k,3),finlev)
enddo
enddo
enddo
write(6,*) 'new finlev',finlev

endif

if(ly > 0) then
do i=1,ly
write(6,*) 'grid', i
write(6,*) tmin(i,:) 
write(6,*) tmax(i,:)
enddo
else
write(6,*) 'nothing to regrid'
endif 

grid=>grid%neighbor 

if(.not.associated(grid)) then
   write(6,*) 'no more grids at this level',level
   write(6,*) igrid,' of',ngd
endif

enddo !igrid
enddo !level

if(allocated(ctrest)) deallocate(ctrest)
return
end







