anchor.diff:         ! compute the sound speed
          cs = sqrt(gamma*p/r)
   
  -       ! compute the eigenvalues
  -       eval(1) = ux - cs
  -       eval(2) = ux
  -       eval(3) = ux + cs
  -
  -
  -       ! compute the left eigenvectors       
  -       lvec(1,:) = [ 0.0_dp_t, -0.5_dp_t*r/cs, 0.5_dp_t/(cs*cs)  ]   ! u - c
  -       lvec(2,:) = [ 1.0_dp_t, 0.0_dp_t,       -1.0_dp_t/(cs*cs) ]   ! u
  -       lvec(3,:) = [ 0.0_dp_t, 0.5_dp_t*r/cs,  0.5_dp_t/(cs*cs)  ]   ! u + c
  -
  -       ! compute the right eigenvectors
  -       rvec(1,:) = [ 1.0_dp_t, -cs/r,    cs*cs    ]   ! u - c
  -       rvec(2,:) = [ 1.0_dp_t, 0.0_dp_t, 0.0_dp_t ]   ! u
  -       rvec(3,:) = [ 1.0_dp_t, cs/r,     cs*cs    ]   ! u + c
  +
  +       ! get the eigenvalues and eigenvectors
  +       call eigen(r, ux, p, cs, lvec, rvec, eval)
   
   
          ! Define the reference states (here xp is the right interface
  
positive.code:        ! first compute beta_xm and beta_xp -- these are the
         ! coefficients to the right eigenvectors in the eigenvector
         ! expansion (see Colella 1990, page 191)
         do m = 1, nprim
            sum = 0.0_dp_t
  
            ! dot product of the current left eigenvector with the
            ! primitive variable jump
            do n = 1, nprim
               sum = sum + lvec(m,n)*dQ(n)
            enddo
  
            ! here the sign() function makes sure we only add the right-moving
            ! waves
  
negative.code: 
      dtdx = dt/U%grid%dx
  
      ! The basic idea here is that we do a characteristic
      ! decomposition.  The jump in primitive variables (Q) can be
      ! transformed to a jump in characteristic variables using the left
  