anchor.diff:                                                            ! proj is symetric
             real*8, dimension(npro,nsd,nsd) :: pt_g_p0, pt_g_p1 ! proj_{ik} g^ij proj_{jm}
             integer :: iel
  +c... get the tangential projector          
  +          call calc_projector(proj, nv0) ! get the tangential projector
  +c... get the c^h
  +          call calc_ch(ch0, ch1, f_jump)
   c          
  -          call calc_projector(proj, nv0) ! get the tangential projector
             call e3giju (giju0, dxidx0, npro, nsd, lcsyst) ! get g^{ij}, 0 side
             call e3giju (giju1, dxidx1, npro, nsd, lcsyst) ! get g^{ij}, 0 side
   c.... get the full g^ij matrix
  
positive.code:           real*8, dimension(npro,nflow) :: temp0, temp1 ! local temporary array
            integer :: iel                                                  
  c
            factor = one
  c
            u_ref_0(:,1) = 1.0d0 ! hacking, rho_ref, gas phase
            u_ref_0(:,2) = 1.0d0 * 1.0d2 ! hacking, rho_ref*v_ref
  
negative.code:           do iel = 1,npro
              pt_g_p0(iel,:,:) = matmul(pt_g0(iel,:,:),proj(iel,:,:))
              pt_g_p1(iel,:,:) = matmul(pt_g1(iel,:,:),proj(iel,:,:))             
            enddo         
  c... calculate the local residual
            call e3if_dc_res(ri0, var0, A0_0, ch0,pt_g_p0)
            call e3if_dc_res(ri1, var1, A0_1, ch1,pt_g_p1)
  c... calculate the local stiffness matrix
            call e3if_dc_egmass()
            call e3if_dc_egmass()          
  