anchor.diff:            real*8, dimension(npro,nflow) :: temp0, temp1 ! local temporary array
             integer :: iel                                                  
   c
  -          factor = zero
  -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
  -          u_ref_0(:,3) = 1.0d0 * 1.0d2 ! hacking, rho_ref*v_ref
  -          u_ref_0(:,4) = 1.0d0 * 1.0d2 ! hacking, rho_ref*v_ref
  -          u_ref_0(:,5) = 1.0d0 * (1.0d2)**two ! hacking, rho_ref*v_ref^2
  -c
  -          u_ref_1(:,1) = 1.0d2 ! hacking, rho_ref, gas phase
  -          u_ref_1(:,2) = 1.0d2 * 1.0d0 ! hacking, rho_ref*v_ref
  -          u_ref_1(:,3) = 1.0d2 * 1.0d0 ! hacking, rho_ref*v_ref
  -          u_ref_1(:,4) = 1.0d2 * 1.0d0 ! hacking, rho_ref*v_ref
  -          u_ref_1(:,5) = 1.0d2 * (1.0d0)**two ! hacking, rho_ref*v_ref^2
  +          factor = one
  +c
  +          u_ref_0(:,1) = 1.000000000000000d0 ! hacking, rho_ref, gas phase
  +          u_ref_0(:,2) = 1.000000000000000d2 ! hacking, rho_ref*v_ref
  +          u_ref_0(:,3) = 1.000000000000000d2 ! hacking, rho_ref*v_ref
  +          u_ref_0(:,4) = 1.000000000000000d2 ! hacking, rho_ref*v_ref
  +          u_ref_0(:,5) = 2.114165517241379d5 ! hacking, rho_ref*v_ref^2
  +c
  +          u_ref_1(:,1) = 1.000000000000000d2 ! hacking, rho_ref, gas phase
  +          u_ref_1(:,2) = 1.000000000000000d2  ! hacking, rho_ref*v_ref
  +          u_ref_1(:,3) = 1.000000000000000d2 ! hacking, rho_ref*v_ref
  +          u_ref_1(:,4) = 1.000000000000000d2 ! hacking, rho_ref*v_ref
  +          u_ref_1(:,5) = 4.204805000000000d7! hacking, rho_ref*v_ref^2
   c... diag(U1_ref ... U5_ref) * flux_jump          
             temp0(:,1) = ( one/u_ref_0 (:,1) )*f_jump(:,1)
             temp0(:,2) = ( one/u_ref_0 (:,2) )*f_jump(:,2)
  
positive.code: 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
  
negative.code:           do iel = 1,npro
              pt_g0(iel,:,:) = matmul(proj(iel,:,:), giju_f0(iel,:,:))
              pt_g1(iel,:,:) = matmul(proj(iel,:,:), giju_f1(iel,:,:))
  c... calculate proj_{ik} g^ij proj_{jm}
              pt_g_p0(iel,:,:) = matmul(pt_g0(iel,:,:),proj(iel,:,:))
              pt_g_p1(iel,:,:) = matmul(pt_g1(iel,:,:),proj(iel,:,:))             
            enddo         
  