anchor.diff:  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,:, :)                  
  +     &        + ch(iel) * shga_g_shgb(iel) * A0(iel,:, :) * WdetJ(iel)
                 enddo
               enddo 
             enddo
  
positive.code:           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
  
negative.code:           call e3if_dc_res(ri1, var1, A0_1, ch1,pt_g_p1)
  c... calculate the local stiffness matrix
            if (lhs_dg .eq. 1) then
              call e3if_dc_egmass(egmass00, nshl0, shg0,A0_0, ch0, pt_g_p0,
       &                          WdetJif0)
              call e3if_dc_egmass(egmass11, nshl1, shg1,A0_1, ch1, pt_g_p1,
       &                          WdetJif1)       
            endif          
  c                                                
            
  