c------------------------------------------------------------
c        Update code that is interfaced to the RNPL-
c        generated update routine, via the UPDATES statement 
c        in burgers_rnpl (note that RNPL supplies the routine 
c        header, declaration of routine arguments, as well 
c        as the RETURN and END statements).
c
c        Implements second-order (in time) Runge Kutta time
c        stepping algorithm.
c------------------------------------------------------------
         character*4   cdnm
         parameter   ( cdnm = 'step' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       i,  Nx

         real*8        delta_t,   delta_tb2

         logical       independent_residual
         parameter   ( independent_residual = .true. )

c============================================================
c============================================================
c-------------------------------------------------------------
c Use a more convenient variable for the number of grid points.
c Notice the RNPL variable is g1_nx.
c-------------------------------------------------------------

         nx = g1_nx
         delta_t   = lambda * (x(2)-x(1))
         delta_tb2 = 0.5d0 * delta_t
         if ( ltrace ) then
            write(0,*) cdnm,':  Nx=',Nx
            write(0,*) cdnm,':  Ng=',Ng
            write(0,*) cdnm,':  dt=',delta_t
         endif

c-------------------------------------------------------------
c First step. Computes quantities at half-advanced time step.
c On return from this section, q_npl contains an estimate of 
c the dynamical variables at the half step.
c-------------------------------------------------------------

         call calc_flux(q_n,   x, FF_n, Nx, Ng)
         call update_q (q_np1, q_n, x, FF_n,  delta_tb2, Nx, Ng)
         call update_boundary(q_np1, Nx, Ng)

c-------------------------------------------------------------
c Second step. Computes quantities at full-advanced time step.
c Numerical fluxes are computed using half-step values.
c-------------------------------------------------------------

         call calc_flux(q_np1, x, FF_n, Nx, Ng)
         call update_q(q_np1, q_n, x, FF_n,  delta_t,   Nx, Ng)
         call update_boundary(q_np1, Nx, Ng)

c-------------------------------------------------------------
c If "independent residual" option is enabled, compute said
c residuals.
c-------------------------------------------------------------
         if ( independent_residual ) then
             call indep_res( q_np1, q_n, x, res_n, Nx, Ng, lambda)
         endif
         return