anchor.diff:         Q_l%data(i+1,iqdens) = 0.0_dp_t
          Q_r%data(i,iqdens) = 0.0_dp_t
   
  -       do n = 1, nprim
  +       do n = 1, nwaves
             if (eval(n) >= 0.0_dp_t) then
                Q_l%data(i+1,iqdens) = Q_l%data(i+1,iqdens) + &
  -                  beta_xp(n)*rvec(n,1)
  +                  beta_xp(n)*rvec(n,iqdens)
             endif
   
             if (eval(n) <= 0.0_dp_t) then
                Q_r%data(i,iqdens) = Q_r%data(i,iqdens) + &
  -                  beta_xm(n)*rvec(n,1)
  +                  beta_xm(n)*rvec(n,iqdens)
             endif
          enddo
   
  
positive.code:        ! can reach.  Do the same from the right interface in the cell,
         ! defining Iplus.  See Almgren et al. 2010 (Eq. 30) or Colella
         ! & Sekora (2008), or Miller & Colella (2002), Eq. 90.
         do m = 1, nprim
            sigma = abs(eval(m))*dtdx
            do n = 1, nprim
  
  
negative.code:     type(gridvar_t) :: Q
      type(gridedgevar_t) :: Q_l, Q_r
  
      real (kind=dp_t) :: rvec(nwaves,nprim), lvec(nwaves,nprim), eval(nwaves)
  
      real (kind=dp_t) :: r, ux, p, cs
      real (kind=dp_t) :: beta_xm(nwaves), beta_xp(nwaves)
  
      real (kind=dp_t) :: dq0, dqp
      real (kind=dp_t) :: Iminus(nwaves,nprim), Iplus(nwaves,nprim)
      real (kind=dp_t) :: Qref_xm(nprim), Qref_xp(nprim)
  
      type(gridedgevar_t) :: Qminus, Qplus
  