anchor.diff:              (PI*elem_field(ne_radius,ne)**4) !laminar resistance
   
          ! element turbulent resistance (flow in bifurcating tubes)
  -      
          reynolds=DABS(elem_field(ne_flow,ne)*2.d0*density/ &
               (PI*elem_field(ne_radius,ne)*viscosity))
          zeta = MAX(1.d0,dsqrt(2.d0*elem_field(ne_radius,ne)* &
               reynolds/elem_field(ne_length,ne))*gamma)
          elem_field(ne_resist,ne) = resistance * zeta
  -       write(*,*),"TESTING RESISTANCE",ne,elem_field(ne_resist,ne),elem_field(ne_radius,ne),zeta
  +       !write(*,*),"TESTING RESISTANCE",ne,elem_field(ne_resist,ne),elem_field(ne_radius,ne),elem_ordrs(2,ne)
       enddo 
   
       call enter_exit(sub_name,2)
  
positive.code:                     IF(nrc.NE.0) THEN
                         IF(ny.GT.variable_totals(nrc,nc)) variable_totals(nrc,nc)=ny
                      ENDIF
                      variables_to_node(np,nrc,nc)=ny
                      IF(nrc.NE.1.OR.nc.EQ.1) THEN
                         ! Don't set node_to_variables for nrc=1(rows), nc<>1(LHS) cases
                         node_to_variables(0,ny,nrc)=1 !mesh dof is nodal based
                         node_to_variables(1,ny,nrc)=np
                         node_to_variables(2,ny,nrc)=nc
                      ENDIF !nrc.NE.1
                   ENDIF !variables_to_node.EQ.0
                ENDDO !nn (np)
                ny=ny+1
                IF(ny.GT.mesh_dof) THEN
  
negative.code:     call enter_exit(sub_name,1)
       
       variable_totals = 0
       mesh_from_variables = 0
       variables_at_elem = 0
       variables_at_node = 0
       
  !nrc = loops from 0,1,2
  
  !  Set up mapping arrays for current region:
  !  includes nested loop creating variables_at_node & variables_at_elem consecutively
       DO nrc=0,2 !row or column no
          ny=0 !variable tag
          DO nc=1,2 !no of dependent variable types
  