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