anchor.diff:  c
           end subroutine flux_jump
   c
  -        subroutine stability_term(ri,Kij)
  +        subroutine kinematic_condition(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 :: this_kcy(npro)
              real*8, dimension(npro,nflow,nflow,nsd,nsd) :: CKij
  -c
  -           call calc_CKij(CKij,Kij)
  +           real*8,dimension(npro,nflow,nsd) :: cy_jump, kcy
  +c
  +           do iflow = 1,nflow
  +c
  +             cy_jump(:,iflow,:) = zero
  +c
  +             do jflow = 1,nflow
  +               cy_jump(:,iflow,1) = cy_jump(:,iflow,1) + cmtrx(:,iflow,jflow)*y_jump(:,jflow,1)
  +               cy_jump(:,iflow,2) = cy_jump(:,iflow,2) + cmtrx(:,iflow,jflow)*y_jump(:,jflow,2)
  +               cy_jump(:,iflow,3) = cy_jump(:,iflow,3) + cmtrx(:,iflow,jflow)*y_jump(:,jflow,3)
  +             enddo
  +c
  +           enddo
   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
  +             this_kcy = zero
  +c
  +             do jflow = 1,nflow
  +               do jsd = 1,nsd
  +                 this_kcy = this_kcy + Kij(:,iflow,jflow,isd,jsd)*cy_jump(:,jflow,jsd)
                  enddo
  -c
  -               ri(:,nflow*(isd-1)+iflow) = ri(:,nflow*(isd-1)+iflow) 
  -     &          + pt50 * s * this_sum
  +             enddo
  +c
  +             ri(:,nflow*(isd-1)+iflow) = ri(:,nflow*(isd-1)+iflow) + pt50 * s * this_kcy
  +c      write(*,*) 'KINEMATIC: ',ri(1,nflow*(isd-1)+iflow),this_kcy(1)
   c
                enddo
              enddo
   c
  -        end subroutine stability_term
  -c
  -        subroutine kinematic_condition(ri,y0,y1)
  +        end subroutine kinematic_condition
  +c
  +        subroutine dg_penalty(ri,y0,y1)
   c
              real*8, dimension(:,:), intent(inout) :: ri
              real*8, dimension(:,:), intent(in) :: y0,y1
  
positive.code: c
                call set_lhs_matrices
  c
               call calc_egmass(egmass00,egmass01,
       &                         A0_0, A0_1, Ai0, Ai1,
       &                         Kij0, Kij1,
       &                         AiNa0,AiNa1,KijNaj0,KijNaj1,KijNajC0,KijNajC1,
       &                         shp0,nv0,nv1,WdetJif0,prop0,nshl0,nshl1)
  
               call calc_egmass(egmass11,egmass10,
       &                         A0_1, A0_0, Ai1, Ai0,
       &                         Kij1, Kij0,
       &                         AiNa1,AiNa0,KijNaj1,KijNaj0,KijNajC1,KijNajC0,
       &                         shp1,nv1,nv0,WdetJif1,prop1,nshl1,nshl0)
  
  c      call calc_egmass_(egmass00,AiNa0,KijNaj0,KijNaj0,KijNajC0,shp0,shp1,nv0,nv1,WdetJif0,nshl0,nshl1)
  c      call calc_egmass_(egmass01,AiNa1,KijNaj0,KijNaj1,KijNajC1,shp0,shp1,nv0,nv1,WdetJif0,nshl0,nshl1)
  c      call calc_egmass_(egmass10,AiNa0,KijNaj1,KijNaj0,KijNajC0,shp1,shp0,nv1,nv0,WdetJif1,nshl1,nshl0)
  c      call calc_egmass_(egmass11,AiNa1,KijNaj1,KijNaj1,KijNajC1,shp1,shp0,nv1,nv0,WdetJif1,nshl1,nshl0)
  c
              endif
  c
  
negative.code:             call calc_cmtrx
              call calc_y_jump
  c
              call kinematic_condition(ri0,Kij0)
              call kinematic_condition(ri1,Kij1)
  c
  c ... kinematic condition term:
  c     set the mu coeff to the max of the 0,1 materials 
  