anchor.diff:         Q_l%data(i+1,iqxvel) = 0.0_dp_t
          Q_r%data(i,iqxvel) = 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,iqxvel) = Q_l%data(i+1,iqxvel) + &
  -                  beta_xp(n)*rvec(n,2)
  +                  beta_xp(n)*rvec(n,iqxvel)
             endif
   
             if (eval(n) <= 0.0_dp_t) then
                Q_r%data(i,iqxvel) = Q_r%data(i,iqxvel) + &
  -                  beta_xm(n)*rvec(n,2)
  +                  beta_xm(n)*rvec(n,iqxvel)
             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:        Q_l%data(i+1,iqdens) = 0.0_dp_t
         Q_r%data(i,iqdens) = 0.0_dp_t
  
         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,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,iqdens)
            endif
         enddo
  
  