anchor.diff:      logical, intent(in) :: FIX(mesh_dof)
       integer :: NonZeros,MatrixSize
   !local variables
  -    integer :: ne,np1,depvar1,NumDiag,nhs,np2,depvar2,ost2,nzz,nzz_val,&
  -      nzz_row,ost1,ny2c,nn,np,ny2,ne2,noelem2,nonode
  -    integer :: NPLIST(0:num_nodes)
  -    logical :: FLOW_BALANCED
  +    integer :: i,ne,np,fixed_variables, fixed_flows, fixed_pressures
       character(len=60) :: sub_name
       sub_name = 'calc_sparse_size'
       call enter_exit(sub_name,1)
   
    
  -  NPLIST=0
  -  NumDiag=0!number of diagonal entries
  -  !NonZeros=0
  -  nzz=1 !position in SparseCol
  -  nzz_val=1 !position in SparseVal
  -  nzz_row=1 !position in SparseRow
  -  ost1=0!offset
  -  ost2=0!offset
  -
  -  do ne=1,num_elems
  -     !ne=elems(noelem)!elem no
  -     np1=elem_nodes(2,ne) !second node in that element
  -     depvar1=depvar_at_node(np1,1,1) !which entry
  -     if((.NOT.FIX(depvar1)).OR.(elems_at_node(np1,0).LE.1))then !if its not fixed, or it only has one element connected
  -       NumDiag=NumDiag+1!one more diagonal entry
  -        do nhs=1,2 !for each row entry
  -              np2=elem_nodes(nhs,ne) 
  -              depvar2=depvar_at_node(np2,1,1) 
  -              if(.NOT.FIX(depvar2))then !if depvar for node2 is not fixed
  -                 nzz=nzz+1
  -                 if (nhs.EQ.1) THEN !row entry one
  -                    nzz_val=nzz_val+1
  -                 elseif (nhs.EQ.2) THEN !row entry two
  -                    nzz_val=nzz_val+1
  -                 endif
  -             else!notfix depvar2
  -             endif
  -        enddo!for each row entry,nhs
  -           !now looking at flow depvar
  -           ny2c=depvar_at_elem(0,1,ne)
  -           if(.NOT.FIX(ny2c)) THEN !if not fixed
  -              nzz=nzz+1
  -              nzz_val=nzz_val+1
  -           else
  -           endif !FIX
  -        nzz_row=nzz_row+1
  -        !Now flow balance equation
  -        do nn=1,2 !balances at each node of ne
  -           FLOW_BALANCED=.FALSE. !initialise
  -           np=elem_nodes(nn,ne)
  -           do nonode=1,NPLIST(0) !junctions balanced at already
  -              if(np.EQ.NPLIST(nonode)) THEN
  -                 FLOW_BALANCED=.TRUE. !already balanced at np
  -              endif
  -           enddo !nonode
  -           if((elems_at_node(np,0).GT.1).AND.(.NOT.FLOW_BALANCED))THEN !if there is more than one node at a junction an we arent flow balanced
  -              !at an unbalanced junction
  -              do noelem2=1,elems_at_node(np,0) !for elems with
  -                 ne2=elems_at_node(np,noelem2) !subtended branches only
  -                    ny2=depvar_at_elem(1,1,ne2)
  -                    if(.NOT.FIX(ny2))THEN
  -                       nzz=nzz+1
  -                       if(np.EQ.elem_nodes(2,ne2)) THEN !end node
  -                          nzz_val=nzz_val+1
  -                       elseif(np.EQ.elem_nodes(1,ne2)) THEN !start node
  -                          nzz_val=nzz_val+1
  -                       endif
  -                    endif
  -              enddo !noelem2
  -              NPLIST(0)=NPLIST(0)+1 !stores nodes where flow
  -              NPLIST(NPLIST(0))=np !balance been done
  -           endif !elem_from_node(np,0).GT.1
  -           if((.NOT.FLOW_BALANCED).AND.(elems_at_node(np,0).GT.1))THEN
  -              nzz_row=nzz_row+1
  -              NumDiag=NumDiag+1 !# entries on diagonal
  -           endif !.NOT.FLOW_BALANCED
  -        enddo !nn
  -     else
  -     endif
  -
  -  enddo 
  -    NonZeros=nzz-1
  -    MatrixSize=nzz_row-1
  +    fixed_variables = 0
  +    !count fixed variables
  +    do i=1,mesh_dof
  +        if(FIX(i))then
  +            fixed_variables = fixed_variables + 1
  +        endif
  +    enddo
  +    MatrixSize = mesh_dof - fixed_variables
  +
  +    !get count of fixed flows
  +    fixed_flows = 0
  +    do ne=1,num_elems
  +        if(FIX(depvar_at_elem(1,1,ne)))then
  +            fixed_flows = fixed_flows + 1
  +        endif
  +    enddo
  +
  +    fixed_pressures = fixed_variables - fixed_flows
  +
  +    !count of pressure equations = (number of elements * 3 variables in each equation) - fixed pressures - fixed flows
  +    NonZeros = num_elems*3 - fixed_pressures - fixed_flows
  +    !count of conservation of flow equations = sum of elements connected to nodes which have at least 2 connected elements - fixed flows
  +    do np=1, num_nodes
  +        if(elems_at_node(np,0).GT.1)then
  +            NonZeros = NonZeros + elems_at_node(np,0)
  +        endif
  +    enddo
  +    NonZeros = NonZeros - fixed_flows
       call enter_exit(sub_name,2)
     end subroutine calc_sparse_size
   
  
positive.code: !
  !*calc_sparse_1d_tree:* Calculates sparsity structure for 1d tree problems
  
  subroutine calc_sparse_1dtree(density,FIX,grav_vect,mesh_dof,depvar_at_elem,&
          depvar_at_node,NonZeros,MatrixSize,SparseCol,SparseRow,SparseVal,RHS,&
          prq_solution,update_resistance_entries)
  
  
      use indices
  
negative.code:         endif
        else!Need to update just the resistance values in the solution matrix
          do ne=1,num_elems !update for all ne
            if(update_resistance_entries(ne).gt.0)then
              nz=update_resistance_entries(ne)
              SparseVal(nz)=-elem_field(ne_resist,ne) !Just updating resistance
            endif
          enddo
          if(bc_type.eq.'flow')then !update RHS to account for element resistance
            write(*,*) update_flow_nzz_row
            do ne=1,num_elems
              depvar = depvar_at_elem(1,1,ne)
              if(FIX(depvar))then
                RHS(update_flow_nzz_row) = prq_solution(depvar,1)*elem_field(ne_resist,ne)
              endif
            enddo
            !SparseVal(nz)=-elem_field(ne_resist,ne) !Just updating resistance
          endif
        endif!first or subsequent iteration
  !! ----CALL SOLVER----
        call pmgmres_ilu_cr(MatrixSize, NonZeros, SparseRow, SparseCol, SparseVal, &
  