!***********************************************************************
      subroutine timestep
!
!  fnew =
!
!   r = dt/dtold
!
!  The third order Hyman leap frog scheme is stable.  Courant time
!  step less than or about 0.50 - 0.60.
!
! 29.06.92/kg     : included different step lenght in timesteps.
! 07.04.93/aake   : 2nd order first time step (1st order for debug)
! 08.04.93/aake   : simplified Courant
! 13.02.94/aake   : no dt change for nt <= 3 (for debugging)
! 15.02.96/aake   : fixed bug in periodic case in Courant
!
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cwork.inc'

      real, dimension(mx,my,mz)::
     & rold,pxold,pyold,pzold,eold,
     & rsav,pxsav,pysav,pzsav,esav
      common /chyman/
     & rold,pxold,pyold,pzold,eold,
     & rsav,pxsav,pysav,pzsav,esav

      data ifirst/2/
      save ifirst

!-----------------------------------------------------------------------
!  For debugging, force first order when nt  <= 3.
!-----------------------------------------------------------------------
      if (nstep .le. 3) ifirst=1
      call pdez
      call Courant

!-----------------------------------------------------------------------
!  First order time step
!-----------------------------------------------------------------------
      if (ifirst .eq. 1) then

        ifirst = 0
        rold  = r
        r     = r  + dt*drdt
        pxold = px
        px    = px + dt*dpxdt
        pyold = py
        py    = py + dt*dpydt
        pzold = pz
        pz    = pz + dt*dpzdt
        eold  = e
        e     = e  + dt*dedt
        t = t + dt
        return

!-----------------------------------------------------------------------
!  2nd order first step
!-----------------------------------------------------------------------
      elseif (ifirst .eq. 2) then

        ifirst = 0
        rold  = r
        r     = r  + (0.5*dt)*drdt
        pxold = px
        px    = px + (0.5*dt)*dpxdt
        pyold = py
        py    = py + (0.5*dt)*dpydt
        pzold = pz
        pz    = pz + (0.5*dt)*dpzdt
        eold  = e
        e     = e  + (0.5*dt)*dedt
        t = t + 0.5*dt

        call pdez
        r    = rold  + dt*drdt
        px   = pxold + dt*dpxdt
        py   = pyold + dt*dpydt
        pz   = pzold + dt*dpzdt
        e    = eold  + dt*dedt
        t = t + 0.5*dt
        return
      end if

!-----------------------------------------------------------------------
!  Predictor step
!-----------------------------------------------------------------------
      ratio = dt/dtold
      a1 = ratio**2
      b1 = dt*(1.+ratio   )
      a2 = 2.*(1.+ratio   )/(2.+3.*ratio)
      b2 = dt*(1.+ratio**2)/(2.+3.*ratio)
      c2 = dt*(1.+ratio   )/(2.+3.*ratio)

      call pred (r , rold , rsav , drdt , a1, b1, a2, b2)
      call pred (px, pxold, pxsav, dpxdt, a1, b1, a2, b2)
      call pred (py, pyold, pysav, dpydt, a1, b1, a2, b2)
      call pred (pz, pzold, pzsav, dpzdt, a1, b1, a2, b2)
      call pred (e , eold , esav , dedt , a1, b1, a2, b2)
      t = t + dt

!-----------------------------------------------------------------------
!  Corrector step
!-----------------------------------------------------------------------
      call pdez
      r  = rsav  + c2*drdt 
      px = pxsav + c2*dpxdt
      py = pysav + c2*dpydt
      pz = pzsav + c2*dpzdt
      e  = esav  + c2*dedt  
      nflop = nflop+mvar*2
!
      end
!***********************************************************************
      subroutine pdez
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cwork.inc'

      drdt  = 0.
      dpxdt = 0.
      dpydt = 0.
      dpzdt = 0.
      dedt  = 0.
      if (mvar .gt. 5) then
        dBxdt = 0.
        dBydt = 0.
        dBzdt = 0.
      end if

      call pde

      end
!***********************************************************************
      subroutine pred (f, fold, fsav, dfdt, a1, b1, a2, b2)
      include 'cparam.inc'
      real f(mx,my,mz), fold(mx,my,mz),  fsav(mx,my,mz),  dfdt(mx,my,mz)
      
      fsav = f
!      f    = a1*fold + (1.-a1)*f + b1*dfdt
      f    = f + a1*(fold-f) + b1*dfdt
      fold = fsav
      fsav = f + a2*(fsav-f) + b2*dfdt
      nflop = nflop+8
      end