      subroutine propagate(lp,crsstp,finlev)
      use amr
      type(LEVEL_INFO),dimension(maxlev) :: lp     
      integer crsstp
    
      integer, parameter :: stacksize = 1000
      integer steps,curlev,finlev
      integer, array(stacksize) :: stack
      real*8 xmaxim
      integer, array(3):: imaxim  
 

      INTERFACE 
       subroutine timestep(lp,level)
       use amr
       type(LEVEL_INFO),dimension(maxlev) :: lp     
       integer level     
       end subroutine timestep

       subroutine regrid(lp,curlev,finlev,first)
       use amr
       type(LEVEL_INFO), dimension(maxlev) :: lp
       integer curlev,finlev
       logical first
       end subroutine regrid
       
      END INTERFACE

      finlev =1
      maxsteps=1

      curlev=crslev
      stack(1)=curlev
      stack(2)=crsstp
      top=2

      DO WHILE(top.gt.0)
       curlev=stack(top-1)
       steps =stack(top)
c       write(6,*) 'steps to be executed',steps,' on level',curlev
       
       if(steps.gt.0) then

         IF(curlev.lt.maxlev. AND. lp(curlev)%lsteps > 0) THEN
            IF(MOD(lp(curlev)%lsteps,lp(curlev)%lrgint)==0) THEN
               IF(curlev==crslev) THEN
c         write(6,*) ' ========== REGRID TIME (BASE GRID) ========='
c         write(6,*) 'step@level',lp(curlev)%lsteps,curlev
c         write(6,*) 'PROPAGATE: current status' 
c         do iii=1,maxlev
c         if(associated(lp(iii)%ml)) then
c           write(6,*) 'level',iii,' exists'
c         else
c           write(6,*) 'level',iii,' does not'
c         endif
c         enddo       
                  call regrid(lp,curlev,finlev,.false.)
c         write(6,*) 'PROPAGATE: current status AFTER' 
c         do iii=1,maxlev
c         if(associated(lp(iii)%ml)) then
c           write(6,*) 'level',iii,' exists'
c         else
c           write(6,*) 'level',iii,' does not'
c         endif
c         enddo       
c         write(6,*) ' =========== END REGRIDDING ================'
               ELSE
                 IF(MOD(lp(curlev-1)%lsteps-1,lp(curlev-1)%lrgint) <> 0
     &              .OR.      
     &            int(abs(lp(curlev-1)%ml%grid%gfcn(1)%t(1) - 
     &            lp(curlev-1)%ldt -  
     &            lp(curlev)%ml%grid%gfcn(1)%t(1))/lp(curlev)%ldt) <> 0)
     &              THEN
c         write(6,*) ' ========== REGRID TIME (FINE GRID) ========='
c         write(6,*) 'step@level',lp(curlev)%lsteps,curlev
c          write(6,*) 'PROPAGATE: current status' 
c         do iii=1,maxlev
c         if(associated(lp(iii)%ml)) then
c           write(6,*) 'level',iii,' exists'
c         else
c           write(6,*) 'level',iii,' does not'
c         endif
c         enddo       
                  call regrid(lp,curlev,finlev,.false.)
c         write(6,*) ' =========== END REGRIDDING ================'
                  ENDIF 
               ENDIF  
             ENDIF      

         ENDIF 

         stack(top-1)=curlev
         stack(top)=-steps

         call timestep(lp,curlev)

c        xmaxim=maxval(lp(1)%ml%grid%gfcn(PIN)%tlev(1)%data(:,1,1))
c        imaxim=
c     &  maxloc(lp(1)%ml%grid%gfcn(PIN)%tlev(1)%data)
c        write(6,*) 'maxvalue@i:',xmaxim,imaxim(1)+1

        SELECT CASE(curlev)
          CASE(1) 
          write(6,*) '***', 
     &    crsstp-steps+1
          CASE(2) 
          write(6,*) '   ***', 
     &    int(lp(curlev)%rcxshp)-steps+1
          CASE(3) 
          write(6,*) '      ***', 
     &    int(lp(curlev)%rcxshp)-steps+1
          CASE(4)
          write(6,*) '         ***', 
     &    int(lp(curlev)%rcxshp)-steps+1
          CASE(5) 
          write(6,*) '             ***', 
     &    int(lp(curlev)%rcxshp)-steps+1
          CASE DEFAULT 
          write(6,*) 'even finer', 
     &    int(lp(curlev)%rcxshp)-steps+1
        END SELECT
       
         if(curlev.lt.finlev) then
           if(top+2.gt.stacksize) then
            write(6,*) 'stacksize exceeded: emergency stop'
            stop
           endif
           stack(top+1)=curlev+1
          stack(top+2)=int(lp(curlev+1)%rcxshp)
          top=top+2
         endif 

       else
 
         steps=abs(steps)-1       
         IF(steps.gt.0) then
             stack(top-1)=curlev
             stack(top)  =steps   
         ELSE
             top=top-2
             IF (curlev.gt.crslev) THEN
c              write(6,*) 'INJECT'
              call linject(lp(curlev)%ml,curlev,int(lp(curlev)%rcxshp))
             ENDIF
         ENDIF
       endif
      ENDDO

      return
      end

      subroutine timestep(lp,level)
      use amr
      type(LEVEL_INFO),dimension(maxlev) :: lp     
      integer level     

      type(ADAPTIVE_MESH), pointer :: p 


      INTERFACE
      subroutine grid_step(grid,dt,dr,nr,gflb,gfrb)
      use user_param
      use u_interface
      use mesh
      type(GF3TL), pointer :: grid
      real*8               :: dt
      real*8, dimension(3) :: dr
      integer, dimension(3):: nr,gflb,gfrb
      end subroutine grid_step
      END INTERFACE


!     for each grid at this level perform one temporal step

      p=>lp(level)%ml
      if(.not.associated(p)) then
         write(6,*) 'timestep should not be called, level',level
         stop
      endif

      DO WHILE(associated(p))
       call grid_step(p%grid,lp(level)%ldt,lp(level)%ldr,p%gnr,p%gflb,
     & p%gfrb)
       p=>p%neighbor
      END DO

      lp(level)%lsteps=lp(level)%lsteps+1 

      return
      end


      subroutine grid_step(grid,dt,dr,nr,gflb,gfrb)
      use user_param
      use u_interface
      use mesh
      type(GF3TL), pointer :: grid
      real*8               :: dt
      real*8, dimension(3) :: dr
      integer, dimension(3):: nr,gflb,gfrb

      integer igf,it 
      real*8  newt 

      DO igf=1,grid%n_gfcn
        IF(grid%gfcn(igf)%n_t.eq.1) cycle
        DO it=grid%gfcn(igf)%n_t,2,-1
           grid%gfcn(igf)%tlev(it)%data=grid%gfcn(igf)%tlev(it-1)%data
           grid%gfcn(igf)%t(it)=grid%gfcn(igf)%t(it-1)
        ENDDO
      ENDDO

C     call set_boundaries 

       
! 1-d evolution of Choptiuk     

      nx=nr(1)
c
c        Inner (r = 0) boundary condition (regularity) ...
c
       if(gflb(1).eq.flbex) grid%gfcn(PPN)%tlev(1)%data(1,1,1)=0.d0    
c
c        Interior evolution ...
c
       grid%gfcn(W1)%tlev(1)%data=grid%gfcn(PPN)%tlev(1)%data

       call elfppi(
     &   grid%gfcn(PPN)%tlev(1)%data(:,1,1),            !PPNP1
     &   grid%gfcn(PPN)%tlev(2)%data(:,1,1),            !PPN
     &   grid%gfcn(W1)%tlev(1)%data(:,1,1),             !W1
     &   grid%gfcn(PIN)%tlev(2)%data(:,1,1),            !PIN
     &   grid%gfcn(1)%tlev(1)%c(1)%coord,               !R 
     &   grid%gfcn(R2)%tlev(1)%data(:,1,1),             !R**2
     &   grid%gfcn(R3)%tlev(1)%data(:,1,1),             !R**3
     &   nx,dr(1),dt,2,nx-1,1)

c
c        Outer (r = R) boundary condition (outgoing radiation) ...
c
       if(gfrb(1).eq.frbex) 
     &   grid%gfcn(PPN)%tlev(1)%data(nx,1,1) =
     &   ((4.0d0 * grid%gfcn(PPN)%tlev(2)%data(nx,1,1)
     &      - grid%gfcn(W1)%tlev(1)%data(nx,1,1)) / dt +
     &   (4.0d0 * grid%gfcn(PPN)%tlev(1)%data(nx-1,1,1) 
     &      -  grid%gfcn(PPN)%tlev(1)%data(nx-2,1,1) / dr(1)) /
     &   (3.0d0 / dt + 3.0d0 / dr(1) + 2.0d0 / 
     &   grid%gfcn(PPN)%tlev(1)%c(1)%coord(nx)))  

        call adddis(
     &   grid%gfcn(PPN)%tlev(1)%data(:,1,1),  !PPNP1
     &   grid%gfcn(W1)%tlev(1)%data(:,1,1),   !W1
     &   epsdis,3,nx-6)


!      ----------------------------------------------------------
!       ADVANCE PIN
!      ----------------------------------------------------------

c
c      interior evolution
c
       grid%gfcn(W1)%tlev(1)%data=grid%gfcn(PIN)%tlev(1)%data

       call elfpii(
     &   grid%gfcn(PIN)%tlev(1)%data(:,1,1),            !PINP1
     &   grid%gfcn(PIN)%tlev(2)%data(:,1,1),            !PIN
     &   grid%gfcn(W1)%tlev(1)%data(:,1,1),             !W1
     &   grid%gfcn(PPN)%tlev(2)%data(:,1,1),            !PPN
     &   grid%gfcn(1)%tlev(1)%c(1)%coord,               !R 
     &   grid%gfcn(R2)%tlev(1)%data(:,1,1),             !R**2
     &   grid%gfcn(R3)%tlev(1)%data(:,1,1),             !R**3
     &   nx,dr(1),dt,2,nx-1,1)

c
c        Inner (r = 0) boundary condition (regularity) ...
c
       if(gflb(1).eq.flbex) grid%gfcn(PIN)%tlev(1)%data(1,1,1)=
     &   (4.0d0 * grid%gfcn(PIN)%tlev(1)%data(2,1,1) -    
     &    grid%gfcn(PIN)%tlev(1)%data(3,1,1))/3.0d0

c
c        Outer (r = R) boundary condition (outgoing radiation) ...
c
       if(gfrb(1).eq.frbex) grid%gfcn(PIN)%tlev(1)%data(nx,1,1)=
     &   ((4.0d0 * grid%gfcn(PIN)%tlev(2)%data(nx,1,1)
     &      - grid%gfcn(W1)%tlev(1)%data(nx,1,1))  / dt +
     &     (4.0d0 * grid%gfcn(PIN)%tlev(1)%data(nx-1,1,1) 
     &      -  grid%gfcn(PIN)%tlev(1)%data(nx-2,1,1) / dr(1)) /
     &   (3.0d0 / dt + 3.0d0 / dr(1) + 2.0d0 / 
     &   grid%gfcn(PPN)%tlev(1)%c(1)%coord(nx)))  

        call adddis(
     &   grid%gfcn(PIN)%tlev(1)%data(:,1,1),  !PINP1
     &   grid%gfcn(W1)%tlev(1)%data(:,1,1), !W1
     &   epsdis,3,nx-6)


!      ----------------------------------------------------------
!       ADVANCE PHN
!      ----------------------------------------------------------

       call elfph(
     &  grid%gfcn(PHN)%tlev(1)%data(:,1,1),
     &  grid%gfcn(PHN)%tlev(1)%data(:,1,1),
     &  grid%gfcn(PIN)%tlev(2)%data(:,1,1),
     &  nx,dt
     &  )

       call clcmas(
     &  grid%gfcn(MASS)%tlev(1)%data,
     &  grid%gfcn(MOVR)%tlev(1)%data,
     &  grid%gfcn(DMDR)%tlev(1)%data,
     &  grid%gfcn(PPN)%tlev(1)%data,
     &  grid%gfcn(PIN)%tlev(1)%data,
     &  grid%gfcn(1)%tlev(1)%c(1)%coord,nx
     &  )

!     end of Choptiuk's code

C     call user_function

      newt=grid%gfcn(1)%t(1)+dt

      DO igf=1,grid%n_gfcn
        grid%gfcn(igf)%t(1)=newt
      ENDDO


       return
       end

