anchor.diff:              do m = 1,meqn
                  qvec(m) = q(mx+k,j,m)
               enddo
  -            call rpn2_cons(meqn,maux,idir,qvec,
  +            call rpn2_cons(meqn,maux,idir,iface,qvec,
        &              auxvec_center,auxvec_edge,flux)         
               do m = 1,meqn
                  flux1(j,m,k) = flux1(j,m,k) + dt*el1(j)*flux(m)
  
positive.code:             do m = 1,meqn
                 qvec(m) = q(i,1-k,m)
              enddo
              call rpn2_cons(meqn,maux,idir,qvec,
       &               auxvec_center, auxvec_edge,flux)         
              do m = 1,meqn
                 flux2(i,m,k) = flux2(i,m,k) + dt*el2(i)*flux(m)
  
negative.code:       double precision qvec(meqn),flux(meqn)
        double precision auxvec_center(maux), auxvec_edge(maux)
  
        integer i,j,k,m, idir, iface
  
        do j = 1,my
  c        # left side k=0 = in; k=1 = out
          idir = 0
          iface = 0
           do k = 0,1
              do m = 1,maux
  c              # Cell centered values;  Each cell-centered value
  