anchor.diff:         sum_xm = ZERO
          sum_xp = ZERO
          do n = 1, nwaves
  -          sum_xm = sum_xm + beta_xm(n)*rvec(n,1)
  -          sum_xp = sum_xp + beta_xp(n)*rvec(n,1)
  +          sum_xm = sum_xm + beta_xm(n)*rvec(n,iqdens)
  +          sum_xp = sum_xp + beta_xp(n)*rvec(n,iqdens)
          enddo
   
          Q_l%data(i+1,iqdens) = r_xp + sum_xp
  
positive.code:        ldu = ldelta%data(i,iqxvel)
         ldp = ldelta%data(i,iqpres)
  
         dQ(1) = ldr
         dQ(2) = ldu
         dQ(3) = ldp
  
  
         ! compute the sound speed
         cs = sqrt(gamma*p/r)
  
negative.code:        sum_xm = ZERO
         sum_xp = ZERO
         do n = 1, nwaves
            sum_xm = sum_xm + beta_xm(n)*rvec(n,iqxvel)
            sum_xp = sum_xp + beta_xp(n)*rvec(n,iqxvel)
         enddo
  
         Q_l%data(i+1,iqxvel) = u_xp + sum_xp
  