subroutine regrid(base,base_grid,g,level,finest_level,lp)
use adaptive
use parameters
type(ADAPTIVE_MESH),pointer    :: base
real,dimension(:,:,:)          :: base_grid
real,dimension(:,:,:,:,:),target :: g
integer level,finest_level
type(LEVEL_INFO), dimension(maxlevel)   :: lp

real,dimension(size(base_grid,dim=2),size(base_grid,dim=3)) :: estre

!HPF$ PROCESSORS P(2,2)
!HPF$ DISTRIBUTE *(*,BLOCK,BLOCK) onto *P :: base_grid
!HPF$ DISTRIBUTE *(*,*,*,BLOCK,BLOCK) onto *P :: g
!HPF$ ALIGN estre(:,:) with base_grid(*,:,:)

! local variables
type(ADAPTIVE_MESH),pointer    :: grid,grid1,grid2
real, dimension(:,:,:),pointer :: gridpointer      
real, dimension(2)             :: rmin
real, dimension(2)             :: rmax
integer, dimension(2)          :: nr
real, dimension(0:2)           :: dr

!HPF$ DISTRIBUTE (*,BLOCK,BLOCK) :: gridpointer
INTERFACE

  subroutine initialize_child(parent,child,grid)
  use adaptive
  real,dimension(:,:,:) :: parent,child
  type(ADAPTIVE_MESH),pointer :: grid
!HPF$ PROCESSORS P(2,2)
!HPF$ DISTRIBUTE *(*,BLOCK,BLOCK) ONTO *P :: parent,child
  end subroutine initialize_child

END INTERFACE

!First, deallocate old fine structure
!We will take care about transfering data from the old structure to the
!new one later

call deallocate_FINE_LEVELS(level,finest_level,lp)

!look at all meshes in descending order of refinement

!write(6,*) 'start regridding from',finest_level,' to',level
LEVELS: DO ilev=finest_level,level,-1
!write(6,*) 'current level',ilev


! array ESTRE is supposed to store truncation error estimate
! now, "a poor people" solution is implemented instead
   
estre(:,:)=0.0
where(base_grid(1,:,:)>lp(ilev)%threshold) estre=1.0

! identify clusters, i.e. contiguous rectangular regions where
! estre = 1

! the following will not work for levels 3 and finer
! since there is no coordinate system yet introduced

if(ANY(estre>0.5)) then

 nx=size(estre,dim=1)
 ny=size(estre,dim=2)

! do i=1,nx
! do j=1,ny
! if(estre(i,j) > 0.5) write(6,*) i,j,estre(i,j),base_grid(1,i,j)
! enddo
! enddo


! because of particular symmetry of this problem I am cheatting here
! to gain performance

! do i=1,nx
!   if(ANY(estre(i,:) > 0.5)) exit
! enddo
! rmin(1)=float(i)
! do j=i+1,nx
!   if(ALL(estre(i,:) < 0.5)) exit
! enddo
! rmax(1)=float(j-1)

rmin(1)=1.0
rmax(1)=100.00

! write(6,*) rmin(1),rmax(1)
 do i=1,ny
   if(ANY(estre(:,i) > 0.5)) exit
 enddo
 rmin(2)=float(i)
 do j=i+1,ny
   if(ALL(estre(:,j) < 0.5)) exit
 enddo
 rmax(2)=float(j-1)
! write(6,*) rmin(2),rmax(2)


nr=rmax-rmin+1

! do not include boundaries into the region to be regridded

if(abs(rmin(1)-1.0)<1E-8) then
  rmin(1)=rmin(1)+1.0
  nr(1)=nr(1)-1
endif
if(abs(rmax(1)-100.0)<1E-8) then
  rmax(1)=rmax(1)-1.0
  nr(1)=nr(1)-1
endif
if(nr(1)<1) return

if(abs(rmin(2)-1.0)<1E-8) then
  rmin(2)=rmin(2)+1.0
  nr(2)=nr(2)-1
endif
if(abs(rmax(2)-100.0)<1E-8) then
  rmax(2)=rmax(2)-1.0
  nr(2)=nr(2)-1
endif
if(nr(2)<1) then
 write(6,*) 'regrid region is empty!'
 return
endif 

!to be able to interpolate, the size of the region must be at least 3
! in each direction

if(nr(1)< 3) then
   iadd=3-nr(1)
   nr(1)=nr(1)+iadd
   if(rmax(1)+iadd>float(nx)) then
     rmin(1)=rmin(1)-iadd
   else
      rmax(1)=rmax(1)+iadd
   endif
endif

if(nr(2)< 3) then
   iadd=3-nr(2)
   nr(1)=nr(2)+iadd
   if(rmax(2)+iadd>float(nx)) then
     rmin(2)=rmin(2)-iadd
   else
      rmax(2)=rmax(2)+iadd
   endif
endif

write(6,*) 'corrected:',rmin(1),rmax(1)
write(6,*) 'corrected:',rmin(2),rmax(2)

   dr=lp(ilev+1)%dr

!at this point, nr contains number of grid points with resolution
!of the current grid; now increase the resolution:

   nr=nr*rho-1

! because of limitation of my "psudo-interpolation", nr cannot be
! even; this will work fine with rho=2


!does the new grid fits allocated memory for g1,g2,...?

   if(nr(1).gt.mxf) then
     write(6,*) 'fine mesh along x is too large',nr(1)
     stop
   endif
   if(nr(2).gt.myf) then
     write(6,*) 'fine mesh along y is too large',nr(2)
     stop
   endif

!who is my parent?
   grid2=>lp(ilev)%head
   id=0
   DO WHILE (associated(grid2))
      grid1=>grid2
      grid2=>grid2%sibling
      if(grid1%rmin(1)<=rmin(1).and.grid1%rmax(1)>=rmax(1)      &
      .and. grid1%rmin(2)<=rmin(2).and.grid1%rmax(2)>=rmax(2)) then
           id=grid1%grid_id
           nullify(grid2)
      endif
   ENDDO

   grid2=>lp(ilev)%head
   if(id==0) then
     write(6,*) 'parent not found',ilev
     write(6,*) 'kid:',rmin,rmax
     DO WHILE (associated(grid2))
      write(6,*) grid2%rmin,grid2%rmax
      grid2=>grid2%sibling
     ENDDO
     stop   
  endif

   base=>grid1
   t=base%t

   call allocate_ADAPTIVE_MESH(lp,ilev+1,base,grid,t,rmin,rmax,nr)

   gridpointer=>g(ilev+1,grid%grid_id,:,:,:)

!  initialize data
   call initialize_child(base_grid,gridpointer,grid)
endif

ENDDO LEVELS

do ilev=level,maxlevel
if(lp(ilev)%curmax>0) finest_level=ilev
enddo

!write(6,*) 'end of regridding: finest level is',finest_level
return
end

