subroutine advance(t,finlev,icstep)
use gridfunctions
real t
integer istep,i2step,i3step,icstep
integer, parameter :: nsteps=2
integer levels,curlev,finlev
integer ipoint
real t2,t3
type(ADAPTIVE_MESH),pointer :: grid,gridbck
integer istop

INTERFACE
subroutine update(curlev,nx,ny,nz,p)
use hpfprocessors
use param
integer :: curlev,nx,ny,nz
real, dimension(:,:,:,:) :: p

!HPF$ DISTRIBUTE p(*,BLOCK,*,*) onto PROC
end subroutine update

subroutine boundaries(t,curlev,istep,nx,ibound,id,p,pp)
use hpfprocessors
use param
real t
integer curlev,nx,id,istep
integer, dimension(3,2) :: ibound
real, dimension(:,:,:,:) :: p
real, dimension(:,:,:,:) :: pp

!HPF$ DISTRIBUTE p(*,BLOCK,*,*)
!HPF$ DISTRIBUTE pp(*,BLOCK,*,*)
end subroutine boundaries

subroutine inject(level,id,p,pp)
use hpfprocessors
use param
integer level,id
real, dimension(:,:,:,:) :: p
real, dimension(:,:,:,:) :: pp
!HPF$ DISTRIBUTE p(*,BLOCK,*,*)
!HPF$ DISTRIBUTE pp(*,BLOCK,*,*)
end subroutine inject

END INTERFACE

istop=143
!if(icstep > 40) stop
!if(mod(icstep,2) == 1) call regrid(t,finlev)
call regrid(t,finlev)

write(6,*) 'g21',allocated(g21)
write(6,*) 'g22',allocated(g22)
write(6,*) 'g23',allocated(g23)
write(6,*) 'g31',allocated(g31)
write(6,*) 'g32',allocated(g32)
write(6,*) 'g33',allocated(g33)



do istep=1,nsteps


curlev=1
!p=>g1
nx=lp(curlev)%grid%nr(1)
ny=lp(curlev)%grid%nr(2)
nz=lp(curlev)%grid%nr(3)
!write(6,*) 'update level',curlev,nx,ny,nz
call update(curlev,nx,ny,nz,g1)
if(icstep>istop) write(6,*) g1(1,3,1:4,1), 'base grid (11)'
t=t+lp(curlev)%dt
id=1
call base_bound(t,curlev,nx,lp(1)%grid%ibound,id)
if(icstep>istop) write(6,*) g1(1,3,1:4,1), 'with bond (11)'
!pp=>g1
if(finlev > 1) then
  t2=t-lp(curlev)%dt
  curlev=2
if(icstep>istop) write(6,*) 'start processing level',curlev
  do i2step=1,nsteps
if(icstep>istop) write(6,*) 'i2step=',i2step,' on level',curlev   
   t2=t2+lp(curlev)%dt
   if(allocated(g21)) then
    grid=>lp(curlev)%grid 
!    p=>g21
    nx=lp(curlev)%grid%nr(1)
    ny=lp(curlev)%grid%nr(2)
    nz=lp(curlev)%grid%nr(3)
    write(6,*) nx,ny,nz,curlev,grid%grid_id,':21'

it=(lp(2)%grid%rmin(2)-lp(1)%grid%rmin(2))+1

!if(icstep>istop) write(6,*) g1(2,3,it:it+3,1), 'shifted (11)'
!if(icstep>istop) write(6,*) g21(1,3,1:8:2,1), '21: init'
    call update(curlev,nx,ny,nz,g21)
!if(icstep>istop) write(6,*) g21(1,3,1:8:2,1), '21 updated'
    id=1
    call boundaries(t2,curlev,i2step,nx,lp(2)%grid%ibound,id,g21,g1)
it=(grid%rmin(2)-grid%parent%rmin(2))*rho+1
!if(icstep>istop) write(6,*) g21(1,3,1:8:2,1), '21 bound'

!do kk=1,4
!k1=(kk-1)*2+1
!k2=k1+1
!if(icstep>istop) write(6,*) g21(1,3,k1:k2,1)
!enddo

if(icstep>istop) write(6,*) 'offset used',it
    if(associated(grid%child)) then
       gridbck=>grid
       grid=>grid%child
       t3=t2-lp(curlev)%dt
       curlev=curlev+1      
!if(icstep>istop) write(6,*) 'start processing level',curlev
           do i3step=1,nsteps
!if(icstep>istop) write(6,*) 'i3step=',i3step   
           t3=t3+lp(curlev)%dt
           if(allocated(g31)) then
!              p=>g31
              nx=lp(curlev)%grid%nr(1)
              ny=lp(curlev)%grid%nr(2)
              nz=lp(curlev)%grid%nr(3)
write(6,*) '>>> process grid 31'
              if(nx /= size(g31,dim=2)) then
                write(6,*) 'Oops, x-size of g31 is wrong',nx,size(g31,dim=2)
                stop
              endif 
              if(ny /= size(g31,dim=3)) then
                write(6,*) 'Oops, y-size of g31 is wrong',ny,size(g31,dim=3)
                stop
              endif 
              if(nz /= size(g31,dim=4)) then
                write(6,*) 'Oops, z-size of g31 is wrong',nz,size(g31,dim=4)
                stop
              endif 
if(icstep>istop) write(6,*) g31(1,3,1:16:4,1), '31: init'

!do kk=1,4
!k1=(kk-1)*4+1
!k2=k1+3
!if(icstep>istop) write(6,*) g31(1,3,k1:k2,1)
!enddo

              call update(curlev,nx,ny,nz,g31)
              id=1
if(icstep>istop) write(6,*) g31(1,3,1:16:4,1), '31: updated'
              call boundaries(t3,curlev,i3step,nx,lp(3)%grid%ibound,id,g31,g21)
if(icstep>istop) write(6,*) g31(1,3,1:16:4,1), '31: bound'


!do kk=1,4
!k1=(kk-1)*4+1
!k2=k1+3
!if(icstep>istop) write(6,*) g31(1,3,k1:k2,1)
!enddo
           endif

       if(allocated(g31)) g31(2,:,:,:)=g31(1,:,:,:)
       enddo  ! i3step   
!if(icstep>istop) write(6,*) 'i3step do exited'
it=(grid%rmin(2)-grid%parent%rmin(2))*rho+1
if(icstep>istop) write(6,*) g31(1,3,1:16:4,1), '31: inject'
        call inject(curlev,id,g31,g21)
if(icstep>istop) write(6,*) g21(1,3,it:it+6:2,1), '21injected'
       grid=>gridbck
       curlev=curlev-1  
     endif  !associated(grid%child)
   endif  !g21 
   if(allocated(g22)) then
!     p=>g22
    grid=>lp(curlev)%grid
    grid=>grid%neighbor
    nx=grid%nr(1)
    ny=grid%nr(2)
    nz=grid%nr(3)
write(6,*) 'process grid 22'
    write(6,*) nx,ny,nz,curlev,grid%grid_id,':22'
!if(icstep>istop) write(6,*) g22(1,3,1:8:2,1), '22: init'
    call update(curlev,nx,ny,nz,g22)
!if(icstep>istop) write(6,*) g22(1,3,1:8:2,1), '22: updated'
    id=2
    call boundaries(t2,curlev,i2step,nx,grid%ibound,id,g22,g1) 
    if(associated(grid%child)) then
       gridbck=>grid
       t3=t2-lp(curlev)%dt
       curlev=curlev+1      
       grid=>lp(curlev)%grid%neighbor     !warning g32 specific!
       id=grid%grid_id
       if(.not.associated(grid)) then
          write(6,*) 'Oops, grid not associated with g32',curlev
          stop
       endif
       if(id /= 2) then
          write(6,*) 'Oops, id should be 2, and not',id
          stop
       endif
       do i3step=1,nsteps
           t3=t3+lp(curlev)%dt
           if(allocated(g32)) then
!              p=>g32
write(6,*) '>>> process grid 32'
              nx=grid%nr(1)
              ny=grid%nr(2)
              nz=grid%nr(3)
              if(nx /= size(g32,dim=2)) then
                write(6,*) 'Oops, x-size of g32 is wrong',nx,size(g32,dim=2)
                stop
              endif 
              if(ny /= size(g32,dim=3)) then
                write(6,*) 'Oops, y-size of g32 is wrong',ny,size(g32,dim=3)
grid=>lp(curlev)%grid
write(6,*) lp(curlev)%ngd
do i=1,lp(curlev)%ngd
write(6,*) grid%grid_id
write(6,*) grid%nr
write(6,*) grid%rmin
write(6,*) grid%rmax
write(6,*) associated(grid%neighbor)
grid=>grid%neighbor
enddo                
                stop
              endif 
              if(nz /= size(g32,dim=4)) then
                write(6,*) 'Oops, z-size of g32 is wrong',nz,size(g32,dim=4)
                stop
              endif 
if(icstep>istop) write(6,*) g32(1,3,1:8:2,1), '32: init'
              call update(curlev,nx,ny,nz,g32)
if(icstep>istop) write(6,*) g32(1,3,1:8:2,1), '32: updated'
              call boundaries(t3,curlev,i3step,nx,lp(3)%grid%ibound,id,g32,g22)
if(icstep>istop) write(6,*) g32(1,3,1:8:2,1), '32: bound'
           endif
       if(allocated(g32)) g32(2,:,:,:)=g32(1,:,:,:)
       enddo  ! i3step  
        id=2  
        call inject(curlev,id,g32,g22)
       grid=>gridbck
       curlev=curlev-1  
     endif  !associated(grid%child)
   endif 
   if(allocated(g23)) then
write(6,*) 'wwwwwwwwwwwww GRID 23 wwwwwwwwwwwwwww'
!     p=>g23
    grid=>lp(curlev)%grid
    grid=>grid%neighbor%neighbor
              nx=grid%nr(1)
              ny=grid%nr(2)
              nz=grid%nr(3)
write(6,*) 'process grid 23'
!if(icstep>istop) write(6,*) g23(1,3,1:8:2,1), '23: init'
!if(icstep>istop) write(6,*) g23(1,3,:,1), '23: init'
    write(6,*) nx,ny,nz,curlev,grid%grid_id,':23'
    call update(curlev,nx,ny,nz,g23)
!if(icstep>istop) write(6,*) g23(1,3,1:8:2,1), '23: updated'
    id=3
    call boundaries(t2,curlev,i2step,nx,grid%ibound,id,g23,g1) 
    if(associated(grid%child)) then
write(6,*) 'wwwwwwwwwwwwww GRID 33 wwwwwwwwwwww'
       gridbck=>grid
       t3=t2-lp(curlev)%dt
       curlev=curlev+1      
       grid=>lp(curlev)%grid%neighbor%neighbor     !warning g33 specific!
       if(.not.associated(grid)) then
          write(6,*) 'Oops, grid not associated with g33',curlev
          stop
       endif
       id=grid%grid_id
       if(id /= 3) then
          write(6,*) 'Oops, id should be 3, and not',id
          stop
       endif
       do i3step=1,nsteps
           t3=t3+lp(curlev)%dt
           if(allocated(g33)) then
!              p=>g33   
              nx=grid%nr(1)
              ny=grid%nr(2)
              nz=grid%nr(3)
              if(nx /= size(g33,dim=2)) then
                write(6,*) 'Oops, x-size of g33 is wrong',nx,size(g33,dim=2)
                stop
              endif 
              if(ny /= size(g33,dim=3)) then
                write(6,*) 'Oops, y-size of g33 is wrong',ny,size(g33,dim=3)
grid=>lp(curlev)%grid
write(6,*) lp(curlev)%ngd
do i=1,lp(curlev)%ngd
write(6,*) grid%grid_id
write(6,*) grid%nr
write(6,*) grid%rmin
write(6,*) grid%rmax
write(6,*) associated(grid%neighbor)
grid=>grid%neighbor
enddo                
                stop
              endif 
              if(nz /= size(g33,dim=4)) then
                write(6,*) 'Oops, z-size of g33 is wrong',nz,size(g33,dim=4)
                stop
              endif 
!if(icstep>istop) write(6,*) g33(1,3,1:8:2,1), '33: init'
              call update(curlev,nx,ny,nz,g33)
              call boundaries(t3,curlev,i3step,nx,lp(3)%grid%ibound,id,g33,g23)
           endif
       if(allocated(g33)) g33(2,:,:,:)=g33(1,:,:,:)
       enddo  ! i3step   
        call inject(curlev,id,g33,g23)
!stop   
       grid=>gridbck
       curlev=curlev-1  
     endif  !associated(grid%child)
   endif

  if(allocated(g21)) g21(2,:,:,:)=g21(1,:,:,:)
  if(allocated(g22)) g22(2,:,:,:)=g22(1,:,:,:)
  if(allocated(g23)) g23(2,:,:,:)=g23(1,:,:,:)
  enddo !i2step
!if(icstep>istop) write(6,*) 'loop i2step exited'

   if(allocated(g21)) then
!    p=>g21
!    pp=>g1
    id=1
it=(lp(2)%grid%rmin(2)-lp(1)%grid%rmin(2))*rho+1
!write(6,*) g1(1,3,it:it+3,1), 'TEST IT'
!write(6,*) g1(1,3,:,1), 'TEST IT'
    call inject(curlev,id,g21,g1)
!write(6,*) g1(1,3,it:it+3,1), '11 injected'
!write(6,*) g1(1,3,:,1), '11 injected'
    endif

   if(allocated(g22)) then
!    p=>g22
!    pp=>g1
    id=2
    call inject(curlev,id,g22,g1)
    endif

   if(allocated(g23)) then
!    p=>g23
!    pp=>g1
    id=3
    call inject(curlev,id,g23,g1)
    endif

endif !finlev>1
g1(2,:,:,:)=g1(1,:,:,:)
!call dumpdata(finlev)
if(icstep>istop) write(6,*) g1(2,1,1:4,1)
if(icstep>istop) write(6,*) 'end of istep=',istep,'  +++++++++++++++++++++++++++++++++++++'
enddo ! istep

return
end
