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,wgrid
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 > 45) stop
!if(mod(icstep,2) == 1) call regrid(t,finlev)
call regrid(t,finlev)

do istep=1,nsteps

! ---------------------- first advance base grid ----------------------------
curlev=1
nx=lp(curlev)%grid%nr(1)
ny=lp(curlev)%grid%nr(2)
nz=lp(curlev)%grid%nr(3)
call update(curlev,nx,ny,nz,g1)
t=t+lp(curlev)%dt
id=1
call base_bound(t,curlev,nx,lp(1)%grid%ibound,id)


! ---------------------- now advance grids on level 2 -----------------------
if(finlev > 1) then
  t2=t-lp(curlev)%dt
  curlev=2

! ----------------------- grid 21 -------------------------------------------
  do i2step=1,nsteps
   t2=t2+lp(curlev)%dt
   if(allocated(g21)) then
    grid=>lp(curlev)%grid 
    nx=lp(curlev)%grid%nr(1)
    ny=lp(curlev)%grid%nr(2)
    nz=lp(curlev)%grid%nr(3)
    call update(curlev,nx,ny,nz,g21)
    id=1
    call boundaries(t2,curlev,i2step,nx,lp(2)%grid%ibound,id,g21,g1)

! ----------------------- grid 31 on level3, if exisits --------------------

    if(associated(grid%child)) then
     if(grid%child%grid_id==1) then
       gridbck=>grid
       grid=>grid%child
       t3=t2-lp(curlev)%dt
       curlev=curlev+1      
           do i3step=1,nsteps
           t3=t3+lp(curlev)%dt
           if(allocated(g31)) then
              nx=lp(curlev)%grid%nr(1)
              ny=lp(curlev)%grid%nr(2)
              nz=lp(curlev)%grid%nr(3)
              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 

              call update(curlev,nx,ny,nz,g31)
              id=1
              call boundaries(t3,curlev,i3step,nx,grid%ibound,id,g31,g21)
           endif !g31
       if(allocated(g31)) g31(2,:,:,:)=g31(1,:,:,:)
       enddo  ! i3step   
       call inject(curlev,id,g31,g21)

! ---------------------- back to level 2 ------------------------------------
       grid=>gridbck
       curlev=curlev-1  
      endif !id=1 
     endif  !associated(grid%child)
   endif  !g21 

! --------------------   grid 22 --------------------------------------------
   if(allocated(g22)) then
    grid=>lp(curlev)%grid
    grid=>grid%neighbor
    nx=grid%nr(1)
    ny=grid%nr(2)
    nz=grid%nr(3)
    call update(curlev,nx,ny,nz,g22)
    id=grid%grid_id
    if(id /= 2) then
       write(6,*) 'id should be 2, and not',id
       stop
    endif
    call boundaries(t2,curlev,i2step,nx,grid%ibound,id,g22,g1) 

! --------------------- grid 32 on level 3 ---------------------------------

    if(associated(grid%child)) then
       gridbck=>grid
       nullify(grid)
       curlev=curlev+1  
       id=2
       wgrid=>lp(curlev)%grid
       DO WHILE(associated(wgrid))
        if(wgrid%grid_id == id) grid=>wgrid
        wgrid=>wgrid%neighbor
       enddo !while 
       nullify(wgrid)
       if(.not.associated(grid)) then
         write(6,*) 'Oops, grid 2 on level',curlev,' not found'
       else
       
       t3=t2-lp(curlev)%dt
       if(.not.associated(grid)) then
          write(6,*) 'Oops, grid not associated with g32',curlev
          stop
       endif
       id=grid%grid_id
       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
              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)
                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 

              call update(curlev,nx,ny,nz,g32)
              call boundaries(t3,curlev,i3step,nx,grid%ibound,id,g32,g22)
           endif !g32
       if(allocated(g32)) g32(2,:,:,:)=g32(1,:,:,:)
       enddo  ! i3step  
        call inject(curlev,id,g32,g22)
       endif !grid found 
! --------- back to level 2 ---------------------------------------------
       grid=>gridbck
       curlev=curlev-1  
     endif  !associated(grid%child)
   endif !g22

! ----------------- grid 23 ----------------------------------------------
   if(allocated(g23)) then
    grid=>lp(curlev)%grid
    grid=>grid%neighbor%neighbor
              nx=grid%nr(1)
              ny=grid%nr(2)
              nz=grid%nr(3)
    call update(curlev,nx,ny,nz,g23)
    id=grid%grid_id
    if(id /= 3) then
      write(6,*) 'id should be 3, and not,id'
      stop
    endif
    call boundaries(t2,curlev,i2step,nx,grid%ibound,id,g23,g1) 

! ----------------- grid 33 on level 3 ------------------------------------
    if(associated(grid%child)) then
       gridbck=>grid
       nullify(grid)
       curlev=curlev+1      
       id=3
       wgrid=>lp(curlev)%grid
       DO WHILE(associated(wgrid))
        if(wgrid%grid_id == id) grid=>wgrid
        wgrid=>wgrid%neighbor
       enddo !while 
       nullify(wgrid)
       if(.not.associated(grid)) then
         write(6,*) 'Oops, grid 2 on level',curlev,' not found'
       else
       t3=t2-lp(curlev)%dt
       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
              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)
                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 
              call update(curlev,nx,ny,nz,g33)
              call boundaries(t3,curlev,i3step,nx,grid%ibound,id,g33,g23)
           endif !g33
       if(allocated(g33)) g33(2,:,:,:)=g33(1,:,:,:)
       enddo  ! i3step   
        call inject(curlev,id,g33,g23)
        endif !grid found
! -------------- back to level 2 ------------------------------------------

       grid=>gridbck
       curlev=curlev-1  
     endif  !associated(grid%child)
   endif ! g23

  if(allocated(g21)) g21(2,:,:,:)=g21(1,:,:,:)
  if(allocated(g22)) g22(2,:,:,:)=g22(1,:,:,:)
  if(allocated(g23)) g23(2,:,:,:)=g23(1,:,:,:)
  enddo !i2step

   if(allocated(g21)) then
    id=1
    call inject(curlev,id,g21,g1)
    endif

   if(allocated(g22)) then
    id=2
    call inject(curlev,id,g22,g1)
    endif

   if(allocated(g23)) then
    id=3
    call inject(curlev,id,g23,g1)
    endif

endif !finlev>1

! ----------------- back to level 1 ------------------------------------

g1(2,:,:,:)=g1(1,:,:,:)

enddo ! istep

return
end
