anchor.diff:        endif
   #endif
   c
  -            cycle element_loop
   c...UPWIND????
   c
   C      ri0(iel,16:20) = ri0(iel,16:20) + f1n0(1:5)
  
positive.code: c
          end subroutine flux_jump
  c
          subroutine stability_term(ri,Kij)
  c
             real*8, dimension(:,:), intent(inout) :: ri
             real*8, dimension(:,:,:,:,:), intent(in) :: Kij
  c
             integer :: iflow,jflow,kflow,isd,jsd
             real*8 :: this_sum(npro)
             real*8, dimension(npro,nflow,nflow,nsd,nsd) :: CKij
  c
             call calc_CKij(CKij,Kij)
  c
             do iflow = 1,nflow
               do isd = 1,nsd
  c
                 this_sum = zero
  c
                 do jflow = 1,nflow
                   do jsd = 1,nsd
                     this_sum = this_sum + CKij(:,iflow,jflow,isd,jsd)*y_jump(:,jflow,jsd)
                   enddo
                 enddo
  c
                 ri(:,nflow*(isd-1)+iflow) = ri(:,nflow*(isd-1)+iflow) 
       &          + pt50 * s * this_sum
  c
               enddo
             enddo
  c
          end subroutine stability_term
  c
          subroutine kinematic_condition(ri,y0,y1)
  c
             real*8, dimension(:,:), intent(inout) :: ri
             real*8, dimension(:,:), intent(in) :: y0,y1
  
negative.code: c
        deallocate(tmpmu0,tmpmu1)
  c
  CC            call dg_penalty(ri0,y0,y1)
  CC            call dg_penalty(ri1,y1,y0)
  c
  c...LHS calculations...
  c
  