character*4   cdnm
         parameter   ( cdnm = 'step' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       i,  Nx

         real*8        set_floor

         real*8        delta_t,   delta_tb2

         logical       independent_residual
         parameter   ( independent_residual = .true. )

         integer       ret, calc_prim  

c============================================================
c============================================================
c-------------------------------------------------------------
c Use a better 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, tau_np1, etc. contain an 
c estimate of the dynamical variables at the half step.
c-------------------------------------------------------------

c-------------------------------------------------------------
c Tracing output to unit 80, which, by default, will go to 
c an ASCII file called 'fort.80'
c-------------------------------------------------------------

         write(80,*) '***************************************'
         write(80,*) '*** 1. SPATIAL COORDINATES      *******'
         write(80,*) '***************************************'
         call dvdump(x,Nx,'r',80)
         write(80,*)

         write(80,*) '***************************************'
         write(80,*) '*** 1. BEFORE CALL TO CALC_FLUX *******'
         write(80,*) '***************************************'
         call dvdump(rho_n,Nx,'rho_n',80)
         call dvdump(v_n,Nx,'v_n',80)
         call calc_flux(rho_n, v_n,
     &                  x,  FF1_n, FF2_n, Nx, Ng, Gamma, floor)
         write(80,*)

         write(80,*) '***************************************'
         write(80,*) '*** 1. BEFORE CALL TO UPDATE_Q  *******'
         write(80,*) '***************************************'
         call dvdump(FF1_n,Nx,'FF1_n',80)
         call dvdump(FF2_n,Nx,'FF2_n',80)
         call dvdump(tau_np1,Nx,'tau_np1',80)
         call dvdump(tau_n,Nx,'tau_n',80)
         call dvdump(S_np1,Nx,'S_np1',80)
         call dvdump(S_n,Nx,'S_n',80)
         call update_q (tau_np1, S_np1, tau_n, S_n, x, 
     &                  FF1_n, FF2_n, delta_tb2, Nx, Ng)
         write(80,*)

         write(80,*) '***************************************'
         write(80,*) '*** 1. BEFORE CALL TO UPDATE_BOUNDARY *'
         write(80,*) '***************************************'
         call dvdump(tau_np1,Nx,'tau_np1',80)
         call dvdump(tau_n,Nx,'tau_n',80)
         call dvdump(S_np1,Nx,'S_np1',80)
         call dvdump(S_n,Nx,'S_n',80)
         call update_boundary(tau_np1, S_np1, Nx, Ng)
         write(80,*)

         write(80,*) '***************************************'
         write(80,*) '*** 1. AFTER  CALL TO UPDATE_BOUNDARY *'
         write(80,*) '***************************************'
         call dvdump(tau_np1,Nx,'tau_np1',80)
         call dvdump(S_np1,Nx,'S_np1',80)
         write(80,*)

         do i = 1, Nx
            tau_np1(i) = max(tau_np1(i),floor+abs(S_np1(i)))  
            ret= calc_prim(tau_np1(i), S_np1(i), 
     &                     rho_np1(i), v_np1(i), P_np1(i), Gamma)
            if (ret .lt. 0) then
               write(0,*) cdnm, ': Stopped at the 1st R-K iteration'
               write(0,*) cdnm, ': position x(',i,')=', x(i)
               stop
            endif
         enddo

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

         write(80,*) '+++++++++++++++++++++++++++++++++++++++'
         write(80,*) '+++ 2. BEFORE CALL TO CALC_FLUX +++++++'
         write(80,*) '+++++++++++++++++++++++++++++++++++++++'
         call dvdump(tau_np1,Nx,'tau_np1',80)
         call dvdump(S_np1,Nx,'S_np1',80)
         call dvdump(rho_np1,Nx,'rho_np1',80)
         call dvdump(v_np1,Nx,'v_np1',80)
         call dvdump(P_np1,Nx,'P_np1',80)
         write(80,*)
         call calc_flux(rho_np1, v_np1,
     &                  x,  FF1_n, FF2_n, Nx, Ng, Gamma, floor)

         write(80,*) '+++++++++++++++++++++++++++++++++++++++'
         write(80,*) '+++ 2. BEFORE CALL TO UPDATE_Q  +++++++'
         write(80,*) '+++++++++++++++++++++++++++++++++++++++'
         call dvdump(FF1_n,Nx,'FF1_n',80)
         call dvdump(FF2_n,Nx,'FF2_n',80)
         call dvdump(tau_np1,Nx,'tau_np1',80)
         call dvdump(tau_n,Nx,'tau_n',80)
         call dvdump(S_np1,Nx,'S_np1',80)
         call dvdump(S_n,Nx,'S_n',80)
         write(80,*)
         call update_q (tau_np1, S_np1, tau_n, S_n, x, 
     &                  FF1_n, FF2_n,  delta_t, Nx, Ng)

         write(80,*) '+++++++++++++++++++++++++++++++++++++++'
         write(80,*) '+++ 2. BEFORE CALL TO UPDATE_BOUNDARY +'
         write(80,*) '+++++++++++++++++++++++++++++++++++++++'
         write(80,*)
         call dvdump(tau_np1,Nx,'tau_np1',80)
         call dvdump(S_np1,Nx,'S_np1',80)
         write(80,*)
         call update_boundary(tau_np1, S_np1, Nx, Ng)

         write(80,*) '***************************************'
         write(80,*) '*** 2. AFTER  CALL TO UPDATE_BOUNDARY *'
         write(80,*) '***************************************'
         write(80,*)
         call dvdump(tau_np1,Nx,'tau_np1',80)
         call dvdump(S_np1,Nx,'S_np1',80)
         write(80,*) 

         do i = 1, Nx
            tau_np1(i) = max(tau_np1(i),floor+abs(S_np1(i)))  
            ret= calc_prim(tau_np1(i), S_np1(i), 
     &                     rho_np1(i), v_np1(i), P_np1(i), Gamma)
            if (ret .lt. 0) then
               write(0,*) cdnm, ': Stopped at the 2nd R-K iteration'
               write(0,*) cdnm, ': position x(',i,')=', x(i)
               stop
            endif
         enddo   
         write(80,*) '+++++++++++++++++++++++++++++++++++++++'
         write(80,*) '+++ 2. AFTER PRIMIITVE COMPUTATION    +'
         write(80,*) '+++++++++++++++++++++++++++++++++++++++'
         write(80,*)
         call dvdump(tau_np1,Nx,'tau_np1',80)
         call dvdump(S_np1,Nx,'S_np1',80)
         call dvdump(rho_np1,Nx,'rho_np1',80)
         call dvdump(v_np1,Nx,'v_np1',80)
         call dvdump(P_np1,Nx,'P_np1',80)

c        if ( independent_residual ) then
c            call indep_res( q_np1, q_n, x, res_n, Nx, Ng, lambda)
c        endif
         return