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: c          
          end subroutine calc_ch
  c
          subroutine e3if_dc_res(ri, var, A0, ch, pt_g_p)
  c..............................................................................
  c  calculation of the contribution of DC operator to the integrand of residual
  c  for both phases
  c..............................................................................
            use propar_m, only: npro
            use conpar_m, only: nflow
            use global_const_m, only: nsd
            use e3if_param_m, only: var_t 
            implicit none
  c
            real*8, dimension(npro,nflow*(nsd+1)), intent(inout) :: ri
            real*8, dimension(npro,nflow,nflow), intent(in) :: A0
            type(var_t), dimension(npro), intent(in) :: var
            real*8, dimension(npro), intent(in) :: ch
            real*8, dimension(npro,nsd,nsd),intent(in) :: pt_g_p
  c
            real*8, dimension(npro,nflow,nsd) :: A0_grad_y
            integer :: iel, k, k0
  c
            do iel = 1,npro
  c... A0 Y_{,m}          
              A0_grad_y(iel,:,:) = matmul( A0(iel,:,:),
       &                                   transpose(var(iel)%grad_y(:,:)))
  c... c^h g^{km}_{I} A0 Y_{,m}
              do k = 1,nsd
                k0 = (k-1)*nflow !k0+1, starting index in ri for kth direction
                                 !k0+nflow, ending index in ri for kth direction
                ri(iel,k0+1:k0+nflow) = ri(iel,k0+1:k0+nflow)
       &        + ch(iel) 
       &        * matmul( A0_grad_y(iel,:,:), pt_g_p(iel,k,:) )
              enddo
            enddo 
  c                
          end subroutine e3if_dc_res
  c
          subroutine e3if_dc_egmass
  c... To be filled
          end subroutine e3if_dc_egmass
          
          subroutine e3if_dc
  