anchor.diff:  c
               climit = zero
   c            climit = one
  -c            climit = 1.e-1
  +CC            climit = 1.e-1
               alpha_LF(iel) = climit * max(abs(dot_product(u0(iel,:)-um0(iel,:),nv0(iel,:))-c0(iel)),
        &                  abs(dot_product(u1(iel,:)-um1(iel,:),nv1(iel,:))-c1(iel)))
               alpha = alpha_LF(iel)
  
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:               enddo 
              enddo
  c
  CC      ri0(iel,16:20) = ri0(iel,16:20) + pt50*alpha*jump_u(1:5)
  CC      ri1(iel,16:20) = ri1(iel,16:20) - pt50*alpha*jump_u(1:5)
  c      write(*,11) 'ri0 a: ',iel,ri0(iel,16:20)
  c      write(*,11) 'ri1 a: ',iel,ri1(iel,16:20)
        ri0(iel,16:20) = ri0(iel,16:20) + alpha*A0_jump_y(1:5)
  