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: !
  !*boundary_conditions:* Defines boundary conditions for prq problems
   subroutine boundary_conditions(ADD,FIX,bc_type,g_type,grav_vect,density,inletbc,outletbc,&
         variables_to_node,variables_to_elem,prq_solution,mesh_dof)
   use arrays,only: dp,num_elems,num_nodes,elems,elem_nodes,elem_cnct,node_xyz,units,&
          num_units
      use diagnostics, only: enter_exit
  
      integer :: mesh_dof
      integer :: variables_to_elem(0:2,2,num_elems)
      integer :: variables_to_node(num_nodes,0:2,2)
      real(dp) :: prq_solution(mesh_dof,2),inletbc,outletbc,density,grav_vect(3)
      logical:: ADD
      logical :: FIX(mesh_dof)
  
negative.code:    call calculate_resistance(density,gamma,viscosity)
          
  !! Calculate sparsity structure for solution matrices
  !estimate nonzeros
  NonZeros=250000! need to minus num of knowns
  MatrixSize=0 !=number of unknowns
  DO variable=1,mesh_dof
    IF(.NOT.FIX(variable)) THEN
      MatrixSize=MatrixSize+1
    ENDIF
  ENDDO !no2
  
     write(*,*) 'MatrixSize',MatrixSize,num_elems,num_nodes,mesh_dof
      allocate (SparseCol(NonZeros), STAT = AllocateStatus)
      IF (AllocateStatus /= 0) STOP "*** Not enough memory for SparseCol array ***"
      allocate (SparseRow(MatrixSize+1), STAT = AllocateStatus)
      IF (AllocateStatus /= 0) STOP "*** Not enough memory for SparseRow array ***"
      allocate (SparseVal(NonZeros), STAT = AllocateStatus)
      IF (AllocateStatus /= 0) STOP "*** Not enough memory for SparseVal array ***"
      allocate (RHS(MatrixSize), STAT = AllocateStatus)
      IF (AllocateStatus /= 0) STOP "*** Not enough memory for RHS array ***"
  
     call calc_sparse_1dtree(mesh_dof,variables_at_elem,variables_at_node,FIX,NonZeros,MatrixSize,SparseCol,&
  	SparseRow,SparseVal,RHS,prq_solution)!(SparseCol,SparseRow,FIX,SparseVal,RHS,NDIAG,NPLIST,&
     !     variables_at_node,variables_at_elem,nz_esed,prq_solution)
  
      allocate (solver_solution(MatrixSize), STAT = AllocateStatus)
      IF (AllocateStatus /= 0) STOP "*** Not enough memory for solver_solution array ***"
      solver_solution=0
  !! Initialise solution vector based on bcs and rigid vessel resistance
     if(bc_type.eq.'pressure')then
        call tree_resistance(total_resistance)
        write(*,*) 'Total Resistance:',total_resistance
        call initialise_solution(inletbc,outletbc,(inletbc-outletbc)/total_resistance,mesh_dof,prq_solution,&
             variables_at_node,variables_at_elem,FIX)
          !put into initial guess
  	no=0
       	DO variable=1,mesh_dof
          IF(.NOT.FIX(variable)) THEN
             no=no+1
             solver_solution(no)=prq_solution(variable,1) 
          ENDIF
       ENDDO !no2
     else
     endif
  
    !  print *,prq_solution(:,1)  
      print *,solver_solution
  
      call BICGSTAB_LinSolv(MatrixSize,NonZeros,RHS,solver_solution,SparseCol,SparseRow,SparseVal,1.d-9,500)
    !  print *,solver_solution(:)
  ! Put solution values into YP
       no=0
       DO variable=1,mesh_dof
          IF(.NOT.FIX(variable)) THEN
             no=no+1
             prq_solution(variable,2)=prq_solution(variable,1) !temp storage of previous solution
             prq_solution(variable,1)=solver_solution(no) !new pressure & flow solutions
             IF(DABS(prq_solution(variable,1)).GT.0.d-6)THEN
                ERR=ERR+(prq_solution(variable,2)-prq_solution(variable,1))**2.d0/prq_solution(variable,1)**2
             ENDIF
          ENDIF
       ENDDO !no2
      write(*,*) 'ERROR', ERR
  
      
      deallocate (mesh_from_variables, STAT = AllocateStatus)
      deallocate (variables_at_elem, STAT = AllocateStatus)
      deallocate (variables_at_node, STAT = AllocateStatus)
      deallocate (prq_solution, STAT = AllocateStatus)
      deallocate (FIX, STAT = AllocateStatus)
      deallocate (solver_solution, STAT = AllocateStatus)
      deallocate (SparseCol, STAT = AllocateStatus)
      deallocate (SparseVal, STAT = AllocateStatus)
      deallocate (SparseRow, STAT = AllocateStatus)
      deallocate (RHS, STAT = AllocateStatus)
      call enter_exit(sub_name,2)
    end subroutine evaluate_prq
  !
  