recursive subroutine propagate(level,base_grid,base,rmin,rmax, &
                                               nr,dr)
use adaptive
use parameters
real, dimension(:,:,:) :: base_grid
type(ADAPTIVE_MESH), pointer :: base
integer     :: level
real, dimension(2) :: rmin
real, dimension(2) :: rmax
integer, dimension(2) :: nr
real, dimension(0:2) :: dr

real, dimension(:,:,:,:), target, allocatable, save :: g1,g2,g3,g4
real, dimension(:,:,:), pointer :: gridpointer
type(ADAPTIVE_MESH), pointer :: grid,child,sibling

!HPF$ PROCESSORS P(2,2)
!HPF$ DISTRIBUTE (*,BLOCK,BLOCK) onto P :: base_grid
!HPF$ DISTRIBUTE (*,BLOCK,BLOCK) onto P :: gridpointer
!!HPF$ DISTRIBUTE (*,*,BLOCK,BLOCK) onto P :: g1,g2,g3,g4

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

write(6,*) 'PROPAGATE at level',level

if(.not.associated(base)) then
write(6,*) '???'
stop
endif


do istep=1,2
write(6,*) 'step',istep,' on level',level

! ================ regridding =====================================
!time for regridding?
if(mod(istep,2)==0 .and. level < maxlevel) then
   write(6,*) 'regrid'
   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)

   write(6,*) 'regrid: create grid', grid%grid_id
   select case(grid%grid_id)
   case (1)
            if(.not.allocated(g1)) &
            allocate(g1(maxlevel,2,nr(1),nr(2)))
            gridpointer=>g1(level+1,:,:,:)
   case (2)
            if(.not.allocated(g2)) &
            allocate(g2(maxlevel,2,nr(1),nr(2)))
            gridpointer=>g2(level+1,:,:,:)
   case (3)
            if(.not.allocated(g3)) &
            allocate(g3(maxlevel,2,nr(1),nr(2)))
            gridpointer=>g3(level+1,:,:,:)
   case (4)
            if(.not.allocated(g4)) &
            allocate(g4(maxlevel,2,nr(1),nr(2)))
            gridpointer=>g4(level+1,:,:,:)
   case default
            write(6,*) 'no more grids allowed at this level'
            stop
   end select

!  initialize data
!   write(6,*) 'call initialize_child'
   call initialize_child(base_grid,gridpointer,grid)
!   write(6,*) gridpointer(1,1:4,2)
!   write(6,*) g1(2,1,1:4,2)
   write(6,*) 'regrid done'
endif
! ====================== end of regridding ===========================

! ====================================================================
! do one step on the current level 

!  write(6,*) 'before',base_grid(1,1:4,1)
  call update(base_grid)
!  write(6,*) 'after',base_grid(1,1:4,1)

! if this grid has childern, update child

  if(associated(base%child)) then
     level=level+1
     child=>base%child 
     write(6,*) 'goto child:',level,child%grid_id
     write(6,*) g1(level,1,1:4,1) 
      select case(child%grid_id)
      case (1)
      gridpointer=>g1(level,:,:,:)
      case (2)
      gridpointer=>g2(level,:,:,:)
      case (3)
      gridpointer=>g3(level,:,:,:)
      case (4)
      gridpointer=>g4(level,:,:,:)
      case default
      write(6,*) 'illegal grid'
      stop
      end select

!      write(6,*) 'call',gridpointer(1,1:4,1)
      call propagate(level,gridpointer,child,child%rmin,child%rmax, &
                    child%nr,child%dr)
      level=level-1 
    endif 

enddo   !two steps on the current level

DO WHILE(associated(base))

! process siblings if any

base=>base%sibling
if(associated(base)) then
      select case(base%grid_id)
      case (1)
      gridpointer=>g1(level,:,:,:)
      case (2)
      gridpointer=>g2(level,:,:,:)
      case (3)
      gridpointer=>g3(level,:,:,:)
      case (4)
      gridpointer=>g4(level,:,:,:)
      case default
      write(6,*) 'illegal grid'
      stop
      end select
endif

enddo !WHILE associated(base)

return
end






