anchor.diff:              do m = 1,meqn
                  qvec(m) = q(i,1-k,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
                  flux2(i,m,k) = flux2(i,m,k) + dt*el2(i)*flux(m)
  
positive.code:             do m = 1,meqn
                 qvec(m) = q(1-k,j,m)
              enddo
              call rpn2_cons(meqn,maux,idir,qvec,
       &               auxvec_center,auxvec_edge,flux)         
              do m = 1,meqn
                 flux0(j,m,k) = flux0(j,m,k) + dt*el0(j)*flux(m)
  
negative.code:          enddo
  
  c        # right side 0 = in; 1 = out
           iface = 1
           do k = 0,1
              do m = 1,maux
  c              # Cell centered values                
  