anchor.diff:     integer :: ne_min,ne_max,nindex,ne,n_max_ord,n,ne_start,&
         inlet_count
      real(dp) :: radius
  -
  - 
      character(len=60) :: sub_name
   
      sub_name = 'define_rad_from_geom'
  
positive.code:     call enter_exit(sub_name,2)
  
    END subroutine define_node_geometry
  
  !
  !###################################################################################
  !
    subroutine define_1d_elements(ELEMFILE)
  !*define_1d_elements:* reads in an element ipelem file to define a geometry
      use arrays,only: elem_direction,elem_field,elems,elem_cnct,elem_nodes,&
           elem_ordrs,elem_symmetry,elems_at_node,elem_units_below,&
           expansile,node_xyz,num_elems,num_nodes
      use indices
      use diagnostics, only: enter_exit
      implicit none
  
      character(LEN=*),intent(in) :: ELEMFILE
      !     Local Variables
      integer :: ibeg,iend,ierror,i_ss_end,j,ne,ne_global,&
           nn,np,np1,np2,np_global
      character(LEN=132) :: ctemp1
      character(LEN=40) :: sub_string
      character(len=60) :: sub_name
  
      sub_name = 'define_1d_elements'
      call enter_exit(sub_name,1)
  
      open(10, file=ELEMFILE, status='old')
  
      read_number_of_elements : do
         read(unit=10, fmt="(a)", iostat=ierror) ctemp1
         if(index(ctemp1, "elements")> 0) then
            call get_final_integer(ctemp1,num_elems)
            exit read_number_of_elements
         endif
      end do read_number_of_elements
  
  !!! allocate memory for element arrays
      if(.not.allocated(elems)) allocate(elems(num_elems))
      if(.not.allocated(elem_cnct)) allocate(elem_cnct(-1:1,0:2,0:num_elems))
      if(.not.allocated(elem_nodes)) allocate(elem_nodes(2,num_elems))
      if(.not.allocated(elem_ordrs)) allocate(elem_ordrs(num_ord,num_elems))
      if(.not.allocated(elem_symmetry)) allocate(elem_symmetry(num_elems))
      if(.not.allocated(elem_units_below)) allocate(elem_units_below(num_elems))
      if(.not.allocated(elems_at_node)) allocate(elems_at_node(num_nodes,0:3))
      if(.not.allocated(elem_field)) allocate(elem_field(num_ne,num_elems))
      if(.not.allocated(elem_direction)) allocate(elem_direction(3,num_elems))
      if(.not.allocated(expansile)) allocate(expansile(num_elems))
  
  !!! initialise element arrays
      elems=0
      elem_nodes=0
      elem_symmetry = 1
      elem_field = 0.0_dp
      expansile = .false.
  
      ne=0
  
      read_an_element : do
         !.......read element number
         read(unit=10, fmt="(a)", iostat=ierror) ctemp1
         if(index(ctemp1, "Element")> 0) then
            call get_final_integer(ctemp1,ne_global) !get element number
            ne=ne+1
            elems(ne)=ne_global
  
            read_element_nodes : do
               read(unit=10, fmt="(a)", iostat=ierror) ctemp1
               if(index(ctemp1, "global")> 0) then !found the correct line
                  iend=len(ctemp1)
                  ibeg=index(ctemp1,":")+1 !get location of first integer in string
                  sub_string = adjustl(ctemp1(ibeg:iend)) ! get the characters beyond : remove leading blanks
                  i_ss_end=len(sub_string) !get the end location of the sub-string
                  ibeg=1
                  do nn=1,2
                     iend=index(sub_string," ") !get location of first blank in sub-string
                     read (sub_string(ibeg:iend-1), '(i7)' ) np_global
                     call get_local_node(np_global,np) ! get local node np for global node
                     elem_nodes(nn,ne)=np ! the local node number, not global
                     sub_string = adjustl(sub_string(iend:i_ss_end)) ! get chars beyond blank, remove leading blanks
                  end do
                  exit read_element_nodes
               endif !index
            end do read_element_nodes
            if(ne.ge.num_elems) exit read_an_element
         endif
  
      end do read_an_element
  
      close(10)
  
      ! calculate the element lengths and directions
      do ne=1,num_elems
         np1=elem_nodes(1,ne)
         np2=elem_nodes(2,ne)
         elem_field(ne_length,ne) = DSQRT((node_xyz(1,np2) - &
              node_xyz(1,np1))**2 + (node_xyz(2,np2) - &
              node_xyz(2,np1))**2 + (node_xyz(3,np2) - &
              node_xyz(3,np1))**2)
         do j=1,3
            elem_direction(j,ne) = (node_xyz(j,np2) - &
                 node_xyz(j,np1))/elem_field(ne_length,ne)
         enddo !j
      enddo
  
      call element_connectivity_1d
      call evaluate_ordering
  
      call enter_exit(sub_name,2)
  
    END subroutine define_1d_elements
  
  !
  !###################################################################################
  !
    subroutine define_rad_from_file(FIELDFILE,TYPE)
  !* define_rad_from_file:* reads in a radius field associated with an aiway tree and assigns radius information
  ! to each element, also calculates volume of each element
      use arrays,only: dp,elem_field,elem_cnct,elem_nodes,&
           elems_at_node,num_elems,num_nodes,node_field
      use indices,only: ne_a_A,ne_length,ne_radius,ne_vol,&
  
negative.code:     call enter_exit(sub_name,2)
  
    end subroutine evaluate_ordering
  !
  !###################################################################################
  !
  !>*set_initial_volume:* assigns a volume to terminal units appended on a tree structure
  !>based on an assumption of a linear gradient in the gravitational direction with max
  !> min and COV values defined.
    subroutine set_initial_volume(Gdirn,COV,total_volume,Rmax,Rmin)
  