subroutine regrid(t,finlev)
use gridfunctions
real t
integer level,finlev
logical, dimension(gnx) :: ledge_x
logical, dimension(gny) :: ledge_y
logical, dimension(gnz) :: ledge_z
integer, dimension(gnx) :: edge_x
integer, dimension(gny) :: edge_y
integer, dimension(gnz) :: edge_z

!HPF$ ALIGN (:) WITH TMPL(:,*) :: ledge_x,edge_x
!HPF$ ALIGN (:) WITH TMPL(*,:) :: ledge_y,edge_y

integer, dimension(1) :: id
real, dimension(2,3)    :: tmin,tmax !two grids; 3d
integer lx,ly,lz

call deallocate_fine_grids

!always start from the base grid

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

if(all(pc(1,:,:,:) < 0.5)) return

!find xmin and xmax

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

tmin=1.0
tmax=gnx
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>2) then
         write(6,*) 'Oops! more than 2 sibling grids!'
         stop
       endif
      else
         exit
endif 
lx=lx+1
if(tmin(lx,1)<1.1 .AND. tmax(lx,1) > gnx-0.1) index=0
enddo
lx=lx-1

! now ymin, ymax

ledge_y=any(p(2,:,:,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:gny-2) edge_y(i)=edge_y(i)-edge_y(i+1)

tmin=1.0
tmax=gny
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>2) then
         write(6,*) 'Oops! more than 2 sibling grids!'
         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 > 2) then
  write(6,*) 'Oops! Cannot handle lx > 2:',lx
  stop
endif
if(ly > 2) then
  write(6,*) 'Oops! Cannot handle ly > 2:',ly
  stop
endif

if(lx == 2 .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 == 2 .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 == 2 .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 == 2 .AND. tmax(1,2) > tmax(2,2)) then
   temp=tmax(1,2)
   tmax(1,1)=tmax(2,2)
   tmax(2,2)=temp
endif

!if rmin is too close to the physical edge, extend the grid
do i=1,lx
if(tmin(i,1) < 2.5) tmin(i,1)=1.0
enddo
do i=1,ly
if(tmin(i,2) < 2.5) tmin(i,2)=1.0
enddo

!for transport, since I do not handle ymax boundary conditions
!proprely, extend ymax

do i=1,ly
if(tmax(i,2) < gny-4.0) then
   tmax(i,2)=tmax(i,2)+4.0
elseif (tmax(i,2) < gny-2.0) then
   tmax(i,2)=tmax(i,2)+2.0
endif
enddo

ngd(level)=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

level=2

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

write(6,*) lp(2)%ngd
if(lp(2)%ngd > 0) write(6,*) lp(2)%grid%nr

!initialize fine grid




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 

!allocate grids
!call init_fine_grid(level)


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

dr1=1.0/(rho**(level-1))

n1=(maxx-minx)/dr1+1
n2=(maxy-miny)/dr1+1
n3=(maxz-minz)/dr1+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
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)

return
end


