anchor.diff:  !
   subroutine calculate_resistance(density,gamma,viscosity)
       use arrays,only: dp,num_elems,num_nodes,elems,elem_nodes,elem_cnct,node_xyz,units,&
  -        num_units,elem_field
  +        num_units,elem_field,elem_ordrs
       use other_consts
       use indices
       use diagnostics, only: enter_exit
       real(dp)::density,gamma,viscosity
   !local variables
  -    integer :: ne,np1,np2
  +    integer :: ne,np1,np2,noelem
       real(dp) :: resistance,reynolds,zeta
       character(len=60) :: sub_name
   
  
positive.code:         ENDIF
       ENDDO
    ELSE !Add terminal pressure BC for all terminal branches
      write(*,*) num_units,units
       DO nonode=1,num_units
          np=elem_nodes(2,units(nonode)) !Second node in element = terminal node
          ny1=variables_to_node(np,1,1) !for fixed pressure BC
          FIX(ny1)=.TRUE. !set fixed
  !! NB// Add gravitational factor in here
          IF(np_in.eq.0) then
  
negative.code:                     IF(nrc.NE.0) THEN
                         IF(ny.GT.variable_totals(nrc,nc)) variable_totals(nrc,nc)=ny
                      ENDIF
                      variables_at_node(np,nrc,nc)=ny
                      IF(nrc.NE.1.OR.nc.EQ.1) THEN
                         ! Don't set mesh_from_variables for nrc=1(rows), nc<>1(LHS) cases
                         mesh_from_variables(0,ny,nrc)=1 !mesh dof is nodal based
                         mesh_from_variables(1,ny,nrc)=np
                         mesh_from_variables(2,ny,nrc)=nc
                      ENDIF !nrc.NE.1
                   ENDIF !variables_at_node.EQ.0
                ENDDO !nn (np)
                ny=ny+1
                IF(ny.GT.mesh_dof) THEN
  