program main
use adaptive

type(ADAPTIVE_MESH), pointer :: grid,parent,base
real, dimension(2) :: rmin 
real, dimension(2) :: rmax 
integer, dimension(2) :: nr 
real, dimension(0:2) :: dr 

real, dimension(:,:), target, allocatable :: base_grid
real, dimension(:,:), target, allocatable :: g21,g22,g23,g24
real, dimension(:,:), target, allocatable :: g31,g32,g33,g34
real, dimension(:,:), target, allocatable :: g41,g42,g43,g44
real, dimension(:,:), pointer :: gridpointer

!HPF$ PROCESSORS P(2,2)
!HPF$ DISTRIBUTE base_grid(BLOCK,BLOCK) onto P
!HPF$ DISTRIBUTE (BLOCK,BLOCK) onto P :: g21,g22,g23,g24
!HPF$ DISTRIBUTE (BLOCK,BLOCK) onto P :: gridpointer

INTERFACE
  subroutine update(a)
  real, dimension(:,:) :: a
!HPF$ PROCESSORS P(2,2)
!HPF$ DISTRIBUTE *(BLOCK,BLOCK) onto *P :: a
  end subroutine update

  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

nullify(base)
!allocate base_grid
nullify(parent)
rmin=(/1.0,1.0/)
rmax=(/100.0,100.0/)
nr=(/100,100/)
dr=(/1.0,1.0,1.0/)
call allocate_ADAPTIVE_MESH(1,parent,base,rmin,rmax,nr,dr)
allocate(base_grid(nr(1),nr(2)))
! initialize base grid
forall(i=1:nr(1),j=1:nr(2)) base_grid(i,j)=float(i)
maxlev=1

do istep=1,10
!time for regridding?
if(mod(istep,2) == 0) then
  level=2
  rmin=(/istep*10,istep*10/)
  rmax=rmin+10
  nr  =(/11,11/)
  dr  =(/0.5,0.5,0.5/)
  call allocate_ADAPTIVE_MESH(level,base,grid,rmin,rmax,nr,dr)

  select case (grid%grid_id)
   case (1)  
            allocate(g21(nr(1),nr(2)))
            gridpointer=>g21
   case (2) 
            allocate(g22(nr(1),nr(2)))
            gridpointer=>g22
   case (3) 
            allocate(g23(nr(1),nr(2)))
            gridpointer=>g23
   case (4) 
            allocate(g24(nr(1),nr(2)))
            gridpointer=>g24
   case default 
            write(6,*) 'no more grids allowed at this level'
            stop
  end select   

!  initialize data

  call initialize_child(base_grid,gridpointer,grid)
!  maxlev=maxlev+1
   maxlev=2
endif
call update(base_grid)
if(maxlev == 2 ) then
! two fine steps
  grid=>base%child
  if(.not.associated(grid)) then
    write(6,*) 'child not associated'
    stop
  endif
 do ifine=1,2
  write(6,*) 'there are', grid%curmax, 'fine grids at this level'
  do igrid=1,grid%curmax 
  select case( igrid)
    case (1)
     gridpointer=>g21
    case (2)
     gridpointer=>g22
    case (3)
     gridpointer=>g23
    case (4)
     gridpointer=>g24
    case default
     write(6,*) 'illegal grid'
     stop
  end select
  write(6,*) 'update grid',igrid
  call update(gridpointer)
  enddo
 enddo  
endif
enddo

stop
end





