anchor.diff:          endif
        enddo
     else !Add terminal pressure BC for all terminal branches
  -     do nonode=1,num_units
  -        np=elem_nodes(2,units(nonode)) !Second node in element = terminal node
  -        ny1=depvar_at_node(np,1,1) !for fixed pressure BC
  -        FIX(ny1)=.TRUE. !set fixed
  +    if(mesh_type.eq.'simple_tree')then !NEED TO SET GRAVITY IN THIS CASE
  +      do nonode=1,num_units
  +         np=elem_nodes(2,units(nonode)) !Second node in element = terminal node
  +         ny1=depvar_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
  -           print *,"Warning --> np_in is not set yet, setting to first node as default"
  -           np_in=1 !Setting to first node as default
  -        endif
  +         if(np_in.eq.0) then
  +            print *,"Warning --> np_in is not set yet, setting to first node as default"
  +            np_in=1 !Setting to first node as default
  +         endif
            grav=0.d0
            do nj=1,3
               grav=grav+density*grav_vect(nj)*9810.d0*(node_xyz(nj,np_in)-node_xyz(nj,np))
            enddo
  -        prq_solution(ny1,1)=outletbc-grav !Putting BC value into solution array       
  -        !print *,"BC----",ny1,outletbc,grav,prq_solution(ny1,1)
  -     enddo
  +         prq_solution(ny1,1)=outletbc-grav !Putting BC value into solution array
  +         !print *,"BC----",ny1,outletbc,grav,prq_solution(ny1,1)
  +       enddo
  +     else
  +       do ne=1,num_elems
  +        !ne=elems(noelem)
  +        if (elem_cnct(1,0,ne) == 0) THEN !EXIT ELEMENT
  +              np=elem_nodes(2,ne)
  +              ny1=depvar_at_node(np,1,1) !for fixed pressure BC
  +              FIX(ny1)=.TRUE. !set fixed
  +              prq_solution(ny1,1)=outletbc !Putting BC value into solution array
  +              np_in=np !inlet node set here, gravity reference to this point
  +          endif
  +        enddo
  +
  +     endif
     endif
       call enter_exit(sub_name,2)
     end subroutine boundary_conditions
  
positive.code:     !pressure (at inlet and outlets)
      !flow (flow at inlet pressure at outlet).
  
  mesh_type='simple_tree'
  vessel_type='elastic_g0_beta'
  mechanics_type='linear'
  bc_type='pressure'
  
negative.code:     real(dp) :: prq_solution(mesh_dof,2),inletbc,outletbc,density,grav_vect(3)
      logical:: ADD
      logical :: FIX(mesh_dof)
      character(len=60) ::bc_type,mesh_type
   
    ! local variables
      integer :: nonode,np,ne,ny1,nj,np_in
  