anchor.diff:          ENDIF
        ENDDO
     ELSE !Add terminal pressure BC for all terminal branches
  -    write(*,*) num_units,units
  +    !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
  +        ny1=variables_at_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
  
positive.code:     use diagnostics, only: enter_exit
      !local variables
      integer :: mesh_dof,variable_types
      integer, allocatable :: node_to_variables(:,:,:)
      integer, allocatable :: variables_to_node(:,:,:)
      integer, allocatable :: variables_to_elem(:,:,:)
      integer, dimension(0:2,2) :: variable_totals
      integer :: num_vars
      integer :: AllocateStatus
  
      real(dp), allocatable :: prq_solution(:,:)
      real(dp) :: viscosity,density,inletbc,outletbc,grav_vect(3),gamma
      logical, allocatable :: FIX(:)
      logical :: ADD=.FALSE.
      character(len=60) :: sub_name,mesh_type,vessel_type,mechanics_type,bc_type,g_type
      integer :: g_dirn
  
      sub_name = 'evaluate_prq'
      call enter_exit(sub_name,1)
  
negative.code:             (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),elem_ordrs(2,ne)
      enddo 
  
      call enter_exit(sub_name,2)
  