!subroutine regrid(base,base_grid,grid,g,level,finest_level,lp, &
!                    base_coordx,base_coordy,coordx,coordy)
subroutine regrid(base,base_grid,grid,g,level,finest_level,lp)

use adaptive
use parameters
type(ADAPTIVE_MESH),pointer    :: base,grid
real,dimension(:,:,:)          :: base_grid
!real, dimension(:) :: base_coordx,base_coordy
real,dimension(:,:,:,:,:),target :: g
!real, dimension(:,:,:) :: coordx,coordy
integer level,finest_level
type(LEVEL_INFO), dimension(maxlevel)   :: lp

real,dimension(:,:), allocatable :: estre

!HPF$ PROCESSORS P(PX,PY)
!HPF$ DISTRIBUTE *(*,BLOCK,BLOCK) onto *P :: base_grid
!HPF$ DISTRIBUTE *(*,*,*,BLOCK,BLOCK) onto *P :: g
!! HPF$ ALIGN base_coordx(:) with *base_grid(*,:,*)
!! HPF$ ALIGN base_coordy(:) with *base_grid(*,*,:)
!! HPF$ ALIGN coordx(*,*,:) with *base_grid(*,:,*)
!! HPF$ ALIGN coordy(*,*,:) with *base_grid(*,*,:)
!HPF$ ALIGN estre(:,:) with base_grid(*,:,:)

!local variables
type(ADAPTIVE_MESH),pointer    :: grid1,grid2,grid3,grid4
real, dimension(:,:,:),pointer :: gridpointer1,gridpointer2
real, dimension(2)             :: rmin
real, dimension(2)             :: rmax
integer, dimension(2)          :: nr
real, dimension(0:2)           :: dr
integer                        :: ssz,fsz

!HPF$ DISTRIBUTE (*,BLOCK,BLOCK) :: gridpointer1,gridpointer2

INTERFACE

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

END INTERFACE



!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,*) 'REGRID: current level',ilev

   call deallocate_FINE_LEVELS(ilev+1,lp)

! array ESTRE is supposed to store truncation error estimate
! now, "a poor people" solution is implemented instead


   fsz = size(base_grid,dim=2)    
   ssz = size(base_grid,dim=3)
   allocate(estre(fsz,ssz))
   estre(:,:) = 0.0
   where(base_grid(1,:,:)>lp(ilev)%threshold) estre=1.0
   
   if(ANY(estre>0.5)) then
    nx=base%nr(1)
    ny=base%nr(2)
    rmin(1)=1.0
    rmax(1)=100.00
    do i=1,ny
      if(estre(5,i) > lp(ilev)%threshold) exit
    enddo
    rmin(2)=float(i)
    do j=i+1,ny
      if(estre(5,j) < lp(ilev)%threshold) exit
    enddo
    rmax(2)=float(j-1)
    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)-mx)<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)-my)<1E-8) then
    if((real(my)-rmax(2))<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!'
      lp(ilev)%curmax=0
      if(ilev==2) lp(3)%curmax=0
      deallocate(estre) 
      exit
    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(2)=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:

    if(ilev==1) nr=nr*rho-1
    if(ilev==2) nr=nr*rho*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

!   grid=>lp(ilev+1)%head

   base=>grid1
   t=base%t

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

!   grid3=>grid
!   if(.not.associated(grid3)) then
!   write(6,*) 'grid not associated in level',ilev
!   stop
!   endif

!   DO WHILE(associated(grid3))
!     coordx(ilev+1,grid3%grid_id,1:(mxf/rho):rho)=base_coordx(1:mx)
!     coordy(ilev+1,grid3%grid_id,1:(myf/rho):rho)=base_coordy(1:my)
!     forall(i=2:(mxf/rho):rho) coordx(ilev+1,grid3%grid_id,i)= &
!                          0.5*(coordx(ilev+1,grid3%grid_id,i-1)+ &
!                             coordx(ilev+1,grid3%grid_id,i+1))
!     forall(j=2:(myf/rho):rho) coordy(ilev+1,grid3%grid_id,j)= &
!                          0.5*(coordy(ilev+1,grid3%grid_id,j-1)+ &
!                             coordy(ilev+1,grid3%grid_id,j+1))
!     grid3=>grid3%sibling
!   ENDDO
!   write(6,*) 'base_coordx(10)',base_coordx(10)
!   write(6,*) 'coordx(2,1,10)',coordx(2,1,10)

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

!  initialize data

!                                                          ret from alloc
!                                     parent(g)  child(g)  child(AM)
   if(ilev==1) call initialize_child(base_grid,gridpointer1,grid)
   if(ilev==2) call initialize_child(g(ilev,1,:,:,:),gridpointer1,grid)
 endif

deallocate(estre) 
ENDDO LEVELS

do ilev=1,maxlevel
write(6,*) lp(ilev)%curmax
if(lp(ilev)%curmax>0) finest_level=ilev
enddo

write(6,*) 'end of regridding',finest_level
return    
end




