      subroutine propagate(lp,crsstp,finlev)
      use amr
      type(LEVEL_INFO),dimension(maxlev) :: lp     
      integer crsstp
    
      integer, parameter :: stacksize = 1000
      integer steps,curlev,finlev
      integer, dimension(stacksize) :: stack
      real   xmaxim
      integer, dimension(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

      if(AVS) then
      write(11,1001) lp(1)%ml%gnr
      write(11,1002) lp(1)%ml%grid%gfcn(RG)%tlev(1)%data
      endif

      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

!          write(6,*) curlev, lp(curlev)%lsteps,lp(curlev)%lrgint

            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       
          write(6,*) 'call regrid'
          write(6,*) lp%ldt
                  call regrid(lp,curlev,finlev,.false.)
         write(6,*) 'PROPAGATE: current status AFTER' 
         do iii=1,maxlev
         if(associated(lp(iii)%ml)) then
           write(6,*) 'level',iii,' exists'
         else
           write(6,*) 'level',iii,' does not'
         endif
         enddo       
         write(6,*) ' =========== END REGRIDDING ================'
         write(6,*) lp%ldt
               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       
c                  call regrid(lp,curlev,finlev,.false.)
c         write(6,*) ' =========== END REGRIDDING ================'
                  ENDIF 
               ENDIF  
             ENDIF      

         ENDIF 

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

         write(6,*) 'call timestep on level ',curlev
     &, lp(curlev)%ml%grid%gfcn(1)%t(1)
         call timestep(lp,curlev)
         write(6,*) 'timestep done'
     &, lp(curlev)%ml%grid%gfcn(1)%t(1)


        IF(AVS) then 
        write(11,1001) lp(1)%ml%gnr
        write(11,1002) lp(1)%ml%grid%gfcn(RG)%tlev(1)%data
        ENDIF

 1001 Format(3I5)
 1002 Format(E15.6)   


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,*) lp(curlev)%ml%grid%gfcn(1)%t,'***', 
     &    crsstp-steps+1
          CASE(2) 
          write(6,*) lp(curlev)%ml%grid%gfcn(1)%t,'   ***', 
     &    int(lp(curlev)%rcxshp)-steps+1
          CASE(3) 
          write(6,*) lp(curlev)%ml%grid%gfcn(1)%t,'      ***', 
     &    int(lp(curlev)%rcxshp)-steps+1
          CASE(4)
          write(6,*) lp(curlev)%ml%grid%gfcn(1)%t,'         ***', 
     &    int(lp(curlev)%rcxshp)-steps+1
          CASE(5) 
          write(6,*) lp(curlev)%ml%grid%gfcn(1)%t,'             ***', 
     &    int(lp(curlev)%rcxshp)-steps+1
          CASE DEFAULT 
          write(6,*) lp(curlev)%ml%grid%gfcn(1)%t,'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                dt
      real, 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
c      write(6,*) 'time_step entered'
      p=>lp(level)%ml
      if(.not.associated(p)) then
         write(6,*) 'timestep should not be called, level',level
         stop
      endif

      DO WHILE(associated(p))
       write(6,*) 'call grid_step at level ',level
       write(6,*) 'HERE ::::::::::::::::::::::::'
       write(6,*) 'DT at level 1&2',lp%ldt
       write(6,*) 'DX at level 1&2',lp%ldr(1)
       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                dt
      real, dimension(3) :: dr,gflb,gfrb
      integer, dimension(3):: nr

      integer igf,it 
      real  newt,dx,dy,dz,t

      
      write(6,*) 'grid_step entered: dt',dt

      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 

      write(6,*) 'GRID STEP nr', nr ,dr
      nx=nr(1)
      ny=nr(2)
      nz=nr(3)
      dx=dr(1)
      dy=dr(2)
      dz=dr(3)
      write(6,*) 'gridstep (dx,dt,nx)' ,dx,dt,nx

C     call user_function
      t=grid%gfcn(1)%t(1) 
      newt=grid%gfcn(1)%t(1)+dt
      write(6,*) 'gridstep: dt',dt,'t',t,'newt',newt 

      DO igf=1,grid%n_gfcn
        grid%gfcn(igf)%t(1)=newt
      ENDDO
     
      write(6,*) 'call Method'
      write(6,*) nx,ny,nz,dt,dx,dy,dz 

      do i=35,47
      write(6,*) i,grid%gfcn(i)%tlev(1)%data(53,53,53)
      enddo

!      ny=size(grid%gfcn(xx)%tlev(1)%data,dim=2) 
!      write(6,*) ny 

      write(6,*) grid%gfcn(rr)%tlev(1)%data(1,1,:) 
!      write(6,*) size(grid%gfcn(alp)%tlev(1)%data,dim=1) 

      call Method(
     &     nx,ny,nz,
     &     dt,dx,dy,dz,
     &     grid%gfcn(xx)%tlev(1)%data,
     &     grid%gfcn(yy)%tlev(1)%data,
     &     grid%gfcn(zz)%tlev(1)%data,
     &     grid%gfcn(rr)%tlev(1)%data,
     &     grid%gfcn(psi)%tlev(1)%data,
     &     grid%gfcn(alp)%tlev(1)%data,
     &     grid%gfcn(cux)%tlev(1)%data,
     &     grid%gfcn(cuy)%tlev(1)%data,
     &     grid%gfcn(cuz)%tlev(1)%data,
     &     grid%gfcn(rg)%tlev(1)%data,
     &     grid%gfcn(uxx)%tlev(1)%data,
     &     grid%gfcn(uxy)%tlev(1)%data,
     &     grid%gfcn(uxz)%tlev(1)%data,
     &     grid%gfcn(uyy)%tlev(1)%data,
     &     grid%gfcn(uyz)%tlev(1)%data,
     &     grid%gfcn(uzz)%tlev(1)%data,
     &     grid%gfcn(gxx)%tlev(1)%data,
     &     grid%gfcn(gxy)%tlev(1)%data,
     &     grid%gfcn(gxz)%tlev(1)%data,
     &     grid%gfcn(gyy)%tlev(1)%data,
     &     grid%gfcn(gyz)%tlev(1)%data,
     &     grid%gfcn(gzz)%tlev(1)%data,
     &     grid%gfcn(qxx)%tlev(1)%data,
     &     grid%gfcn(qxy)%tlev(1)%data,
     &     grid%gfcn(qxz)%tlev(1)%data,
     &     grid%gfcn(qyy)%tlev(1)%data,
     &     grid%gfcn(qyz)%tlev(1)%data,
     &     grid%gfcn(qzz)%tlev(1)%data,
     &     grid%gfcn(dxuxx)%tlev(1)%data,
     &     grid%gfcn(dxuxy)%tlev(1)%data,
     &     grid%gfcn(dxuxz)%tlev(1)%data,
     &     grid%gfcn(dxuyy)%tlev(1)%data,
     &     grid%gfcn(dxuyz)%tlev(1)%data,
     &     grid%gfcn(dxuzz)%tlev(1)%data,
     &     grid%gfcn(dyuxx)%tlev(1)%data,
     &     grid%gfcn(dyuxy)%tlev(1)%data,
     &     grid%gfcn(dyuxz)%tlev(1)%data,
     &     grid%gfcn(dyuyy)%tlev(1)%data,
     &     grid%gfcn(dyuyz)%tlev(1)%data,
     &     grid%gfcn(dyuzz)%tlev(1)%data,
     &     grid%gfcn(dzuxx)%tlev(1)%data,
     &     grid%gfcn(dzuxy)%tlev(1)%data,
     &     grid%gfcn(dzuxz)%tlev(1)%data,
     &     grid%gfcn(dzuyy)%tlev(1)%data,    
     &     grid%gfcn(dzuyz)%tlev(1)%data,
     &     grid%gfcn(dzuzz)%tlev(1)%data)

      write(6,*) 'call Boundaries'
      call Boundaries(
     &     nx,ny,nz,t,
     &     dt,dx,dy,dz,
     &     grid%gfcn(xx)%tlev(1)%data,
     &     grid%gfcn(yy)%tlev(1)%data,
     &     grid%gfcn(zz)%tlev(1)%data,
     &     grid%gfcn(rr)%tlev(1)%data,
     &     grid%gfcn(psi)%tlev(1)%data,
     &     grid%gfcn(alp)%tlev(1)%data,
     &     grid%gfcn(cux)%tlev(1)%data,
     &     grid%gfcn(cuy)%tlev(1)%data,
     &     grid%gfcn(cuz)%tlev(1)%data,
     &     grid%gfcn(rg)%tlev(1)%data,
     &     grid%gfcn(uxx)%tlev(1)%data,
     &     grid%gfcn(uxy)%tlev(1)%data,
     &     grid%gfcn(uxz)%tlev(1)%data,
     &     grid%gfcn(uyy)%tlev(1)%data,
     &     grid%gfcn(uyz)%tlev(1)%data,
     &     grid%gfcn(uzz)%tlev(1)%data,
     &     grid%gfcn(gxx)%tlev(1)%data,
     &     grid%gfcn(gxy)%tlev(1)%data,
     &     grid%gfcn(gxz)%tlev(1)%data,
     &     grid%gfcn(gyy)%tlev(1)%data,
     &     grid%gfcn(gyz)%tlev(1)%data,
     &     grid%gfcn(gzz)%tlev(1)%data,
     &     grid%gfcn(qxx)%tlev(1)%data,
     &     grid%gfcn(qxy)%tlev(1)%data,
     &     grid%gfcn(qxz)%tlev(1)%data,
     &     grid%gfcn(qyy)%tlev(1)%data,
     &     grid%gfcn(qyz)%tlev(1)%data,
     &     grid%gfcn(qzz)%tlev(1)%data,
     &     grid%gfcn(dxuxx)%tlev(1)%data,
     &     grid%gfcn(dxuxy)%tlev(1)%data,
     &     grid%gfcn(dxuxz)%tlev(1)%data,
     &     grid%gfcn(dxuyy)%tlev(1)%data,
     &     grid%gfcn(dxuyz)%tlev(1)%data,
     &     grid%gfcn(dxuzz)%tlev(1)%data,
     &     grid%gfcn(dyuxx)%tlev(1)%data,
     &     grid%gfcn(dyuxy)%tlev(1)%data,
     &     grid%gfcn(dyuxz)%tlev(1)%data,
     &     grid%gfcn(dyuyy)%tlev(1)%data,
     &     grid%gfcn(dyuyz)%tlev(1)%data,
     &     grid%gfcn(dyuzz)%tlev(1)%data,
     &     grid%gfcn(dzuxx)%tlev(1)%data,
     &     grid%gfcn(dzuxy)%tlev(1)%data,
     &     grid%gfcn(dzuxz)%tlev(1)%data,
     &     grid%gfcn(dzuyy)%tlev(1)%data,
     &     grid%gfcn(dzuyz)%tlev(1)%data,
     &     grid%gfcn(dzuzz)%tlev(1)%data)

           
      write(6,*) 'DONE*DONE*DONE*DONE'

C     call user_function

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

c      DO igf=1,grid%n_gfcn
c        grid%gfcn(igf)%t(1)=newt
c      ENDDO
       write(6,*) maxval( grid%gfcn(gxy)%tlev(1)%data),
     &           minval( grid%gfcn(gxy)%tlev(1)%data)

       return
       end




































