anchor.diff:            real*8, dimension(npro), intent(in) :: ch
             real*8, dimension(npro,nsd,nsd),intent(in) :: gI !g^{ij}_{I} = 
                                                              ! proj^T g^{ij} proj
  +          real*8, dimension(npro),intent(in) :: WdetJ
             integer, intent(in)  :: nshl
   c
             real*8, dimension(npro,nsd) :: shga_g
  
positive.code: c
                  egmass(iel, a_row+1:a_row+nflow, b_col+1:b_col+nflow) = 
       &          egmass(iel, a_row+1:a_row+nflow, b_col+1:b_col+nflow)
       &        + ch(iel) * shga_g_shgb(iel) * A0(iel,:, :)                  
                enddo
              enddo 
            enddo
  
negative.code: c                
          end subroutine e3if_dc_res
  c
          subroutine e3if_dc_egmass(egmass,nshl, shg,A0, ch,gI, WdetJ)  
  c..............................................................................
  c  calculation of the contribution of DC operator to the local stiffness matrix
  c  for both phases
  