anchor.diff:    ! at the average of the parent end and child end coordinates. This is used to
     ! improve the topology of generated trees, minimising the impact of 'odd' branching
     !
  -  subroutine smooth_1d_tree(num_elem_start,length_limit)
  +  SUBROUTINE smooth_1d_tree(num_elem_start,length_limit)
       !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_SMOOTH_1D_TREE" :: SMOOTH_1D_TREE
   
  -    use diagnostics, only: enter_exit
  -    implicit none
  -
  -    integer,intent(in) :: num_elem_start
  -    real(dp),intent(in) :: length_limit
  -
  -    integer :: n,ne,ne1,ne2,np,np0,np1,np2,n_smoothing_steps = 2
  -    real(dp) :: new_xyz(3)
  -    character(len=60) :: sub_name
  +    USE diagnostics, ONLY: enter_exit
  +    IMPLICIT NONE
  +
  +    INTEGER,INTENT(in) :: num_elem_start
  +    REAL(dp),INTENT(in) :: length_limit
  +
  +    INTEGER :: n,ne,ne1,ne2,np,np0,np1,np2,n_smoothing_steps = 2
  +    REAL(dp) :: new_xyz(3)
  +    CHARACTER(len=60) :: sub_name
   
       sub_name = 'smooth_1d_tree'
  -    call enter_exit(sub_name,1)
  -
  -    do n = 1,n_smoothing_steps
  -       do ne = num_elems,num_elem_start,-1
  -          if(elem_cnct(1,0,ne).eq.0)then ! terminal, check branch length
  -             if(elem_field(ne_length,ne).lt.0.75_dp*length_limit)then
  +    CALL enter_exit(sub_name,1)
  +
  +    DO n = 1,n_smoothing_steps
  +       DO ne = num_elems,num_elem_start,-1
  +          IF(elem_cnct(1,0,ne).EQ.0)THEN ! terminal, check branch length
  +             IF(elem_field(ne_length,ne).LT.0.75_dp*length_limit)THEN
                   elem_field(ne_length,ne) = 0.75_dp*length_limit
                   np1 = elem_nodes(1,ne) ! the start node
                   np2 = elem_nodes(2,ne) ! the end node
                   node_xyz(:,np2) = node_xyz(:,np1) + elem_direction(:,ne)*0.75_dp*length_limit
  -             else if(elem_field(ne_length,ne).gt.1.5_dp*length_limit)then
  +             ELSE IF(elem_field(ne_length,ne).GT.1.5_dp*length_limit)THEN
                   elem_field(ne_length,ne) = 1.5_dp*length_limit
                   np1 = elem_nodes(1,ne) ! the start node
                   np2 = elem_nodes(2,ne) ! the end node
                   node_xyz(:,np2) = node_xyz(:,np1) + elem_direction(:,ne)*1.5_dp*length_limit
  -             endif
  -          else
  -             if(elem_cnct(1,0,ne).eq.2)then
  +             ENDIF
  +          ELSE
  +             IF(elem_cnct(1,0,ne).EQ.2)THEN
                   ne1 = elem_cnct(1,1,ne)
                   ne2 = elem_cnct(1,2,ne)
                   np0 = elem_nodes(1,ne)
  
positive.code:     V = unit_vector(V)
      W = unit_vector(W)
  
      if(abs(angle).gt.a_lim)then !reduce angle
         if(abs(dabs(angle)-pi).lt.loose_tol)then !reduce angle from 180 degrees
            if(Nth.eq.2)then
               V(1:3) = elem_direction(1:3,ne-1) !direction of sibling
            else if(Nth.eq.1)then
               V(1:3) = elem_direction(1:3,ne+1) !direction of sibling
            endif
            angle_sibling = angle_btwn_vectors(U,V)
            W = vector_for_angle_limit(U,V,a_lim,angle_sibling+a_lim)
         else
            W = vector_for_angle_limit(U,V,a_lim,angle-a_lim)
         endif
         node_xyz(1:3,np) = node_xyz(1:3,np1)+W(1:3)*0.5_dp*length
         elem_field(ne_length,ne) = 0.5_dp*length
         elem_direction(1:3,ne) = W(1:3)
  
      elseif(abs(angle).lt.a_min)then
         if(abs(angle).lt.loose_tol)then !increase angle from 0 degrees
            if(Nth.eq.2)then
               V(1:3) = elem_direction(1:3,ne-1) !direction of sibling
            else if(Nth.eq.1)then
               V(1:3) = elem_direction(1:3,ne+1) !direction of sibling
            endif
            angle_sibling = angle_btwn_vectors(U,V) !angle between branch and sibling
            W = vector_for_angle_limit(U,V,a_min,angle_sibling+a_min)
         else
            W = vector_for_angle_limit(U,V,a_min,angle-a_min)
         endif
         node_xyz(1:3,np) = node_xyz(1:3,np1)+W(1:3)*0.5_dp*length
         elem_field(ne_length,ne) = 0.5_dp*length
         elem_direction(1:3,ne) = W(1:3)
         
      endif !abs(angle).gt.a_lim
      
      call enter_exit(sub_name,2)
  
    end subroutine adjust_branch_angle
    
    
    !###############################################################
    !
    !*branch_to_cofm*: Creates a new branch towards the cofm of a set of seed points. 
    ! Used in a volume-filling branching method to create a branch that runs from a
    ! defined point to some fraction along a line towards the centre of mass of a 
    ! collection of seed points.
    !
    subroutine branch_to_cofm(map_seed_to_elem,nen,np1,COFM,branch_fraction,length_limit,&
      length_parent,shortest_length,candidate_xyz,make_branch)
      !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_BRANCH_TO_COFM" :: BRANCH_TO_COFM
  
      use diagnostics, only: enter_exit
      implicit none
      
      integer :: map_seed_to_elem(*),nen,np1
      real(dp) :: COFM(3),branch_fraction,length_limit,length_parent,&
           shortest_length,candidate_xyz(3)
      logical :: make_branch
  
      !Local variables
      integer :: N,NCLOSEST(100),nd,nsp,NUM_CLOSEST,NUM_ND,number_of_points
      real(dp) :: CLOSEST(100),DIST,L_COFM,LENGTH,MIN_DIST,VECTOR(3)
      
      character(len=60) :: sub_name
  
      sub_name = 'branch_to_cofm'
      call enter_exit(sub_name,1)
  
      number_of_points = 0
      do nd = 1,num_data
         if(map_seed_to_elem(nd).eq.nen) number_of_points=number_of_points+1
      enddo
      candidate_xyz(1:3) = node_xyz(1:3,np1) + branch_fraction*(COFM(1:3)-node_xyz(1:3,np1))
      VECTOR(1:3) = COFM(1:3)-node_xyz(1:3,np1)
      LENGTH = distance_between_points(candidate_xyz,node_xyz(1,np1))
      L_COFM = vector_length(VECTOR)
      VECTOR = unit_vector(VECTOR)
      
      if(LENGTH.ge.LENGTH_LIMIT)then !the branch will not be terminal
         make_branch = .TRUE.
      else
         if(elem_ordrs(no_gen,nen).lt.12)then
            make_branch = .true.
            candidate_xyz(1:3) = node_xyz(1:3,np1)+VECTOR(1:3)*0.33_dp*length_parent
         else
            make_branch = .FALSE.
         endif
      endif !LENGTH.ge.L_LIM
  
      !*** Position branch
      MIN_DIST = 1.0e+6_dp
  
negative.code:           DIRECTION(1:3)=(node_xyz(1:3,np4)-node_xyz(1:3,np2))
            DIRECTION = unit_vector(DIRECTION)
            elem_direction(1:3,ne-1)=DIRECTION(1:3)
  
            U(1:3)=elem_direction(1:3,ne) !direction of a branch
            V(1:3)=elem_direction(1:3,ne-1) !direction of its sibling
            U = unit_vector(U)
            V = unit_vector(V)
            NRML = cross_product(U,V) !calculate branching plane
            NRML = unit_vector(NRML)
  
            !...angle between branching planes
            ANGLE=angle_btwn_vectors(NRML,NRML_PARENT)
  
            IF(ABS(ROTATION_LIMIT-ABS(PI/2.0_dp-ANGLE)).LE.0.001_dp)THEN
               COMPLETE=.TRUE.
            ENDIF
  
            IF(IT.GT.ITMAX)THEN
               WRITE(*,*) 'WARNING!!!! rotation angle = ',ANGLE*180.0_dp/PI
            ENDIF
  
         ENDDO !do while not found
  
         !.......Alternate calculation for the rotation angle
         angle = rotation_angle(np1,np2,np00,np3,np4)
      ELSE
         !        write(*,*) 'Not',ANGLE_BETWEEN*180.0_dp/PI
      ENDIF
  
      CALL enter_exit(sub_name,2)
  
    END SUBROUTINE check_rotation_angle
  
  
    !###############################################################
    !
    !*create_new_node:* sets up arrays for a new mesh node and element.
    !
    SUBROUTINE create_new_node(ne,ne_start,np,np_start,MAKE)
      !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_CREATE_NEW_NODE" :: CREATE_NEW_NODE
  
      USE diagnostics, ONLY: enter_exit
      IMPLICIT NONE
  
      INTEGER :: ne,ne_start,np,np_start
      LOGICAL :: MAKE
  
      !Local variables
      CHARACTER(len=60) :: sub_name
  
      sub_name = 'create_new_node'
      CALL enter_exit(sub_name,1)
  
      IF(MAKE)THEN
         ne=ne+1
         elems(ne) = ne ! store global element number
         elem_nodes(1,ne) = np_start
         elems_at_node(np_start,0)=elems_at_node(np_start,0)+1
         elems_at_node(np_start,elems_at_node(np_start,0))=ne
  
         np = np+1
         nodes(np) = np
         elems_at_node(np,0) = 0 !initialise
         elem_nodes(2,ne) = np !end node of new element
         elems_at_node(np,0) = elems_at_node(np,0)+1
         elems_at_node(np,elems_at_node(np,0)) = ne
  
         elem_cnct(1,0,ne)=0 !initialise number of proximal branches
         IF(ne_start.NE.0)THEN
            elem_cnct(-1,0,ne)=1
            elem_cnct(-1,elem_cnct(-1,0,ne),ne)=ne_start
            elem_cnct(1,0,ne_start)=elem_cnct(1,0,ne_start)+1
            elem_cnct(1,elem_cnct(1,0,ne_start),ne_start)=ne
         ENDIF
      ENDIF
  
      CALL enter_exit(sub_name,2)
  
    END SUBROUTINE create_new_node
  
  
    !###############################################################
    !
    !*group_seeds_with_branch:* groups a set of seed points with the
    ! closest candidate parent branches. reassigns data (seed) points
    ! to the closest ending of branches in the current generation.
    !
    SUBROUTINE group_seeds_with_branch(map_array,num_next_parents,num_seeds_from_elem, &
         num_terminal,local_parent,DISTANCE_LIMIT,FIRST)
      !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_GROUP_SEEDS_WITH_BRANCH" :: GROUP_SEEDS_WITH_BRANCH
  
      USE diagnostics, ONLY: enter_exit
      IMPLICIT NONE
  
      INTEGER :: num_next_parents,local_parent(:),map_array(:),num_seeds_from_elem(*),&
           num_terminal
      REAL(dp),INTENT(in) :: DISTANCE_LIMIT
      LOGICAL :: FIRST
  
      !Local variables
      INTEGER :: i,n,m,nd,nd_min,ne,n_elm_temp,ne_min,noelem,np,np_temp
      INTEGER :: size_map
      INTEGER,ALLOCATABLE :: map_array_copy(:),my_closest(:)
      REAL(dp) :: dist,min_dist
  
      CHARACTER(len=60) :: sub_name
  
      sub_name = 'group_seeds_with_branch'
      CALL enter_exit(sub_name,1)
  
      size_map = SIZE(map_array)
      ALLOCATE(my_closest(size_map))
  
      IF(first)THEN
  !!! for the first seed point allocation, need to make sure that every terminal branch has
  !!! at least two seed points. find the closest two data points to each terminal node, and
  !!! assign to terminals. the 'my_closest' array is used to indicate that the points have
  !!! already been allocated (when calculating the distribution of points below). Also need
  !!! to ensure that the two closest seed points are unique (assigned to only one terminal).
         !       n_closest = 0
         !       do i=1,2
         !          do N=1,num_next_parents
         !             ne_min=local_parent(N)
         !             np_temp=elem_nodes(2,ne_min)
         !             min_dist = 1.0e+10_dp
         !             nd_min = 0
         !             do nd = 1,num_data
         !                if(map_array(nd).ne.0)then
         !                   if(.not.inlist(nd,my_closest))then
         !                      dist = distance_between_points(data_xyz(1,nd),node_xyz(1,np_temp))
         !                      if(dist.lt.min_dist)then
         !                         nd_min = nd
         !                         min_dist = dist
         !                      endif !dist
         !                   endif
         !                endif
         !             enddo
         !             n_closest = n_closest + 1
         !             my_closest(n_closest) = nd_min
         !             map_array(nd_min) = ne_min
         !          enddo
         !       enddo !i
  
      ELSE !use the data groupings from previous
         ALLOCATE(map_array_copy(size_map))
         map_array_copy(1:size_map) = map_array(1:size_map)
         DO n=1,num_next_parents
            ne_min = local_parent(n)
            np_temp = elem_nodes(2,ne_min)
            MIN_DIST=1.0e+10_dp
            DO nd=1,num_data
               IF(map_array(nd).EQ.ne_min)THEN ! was associated with this element
                  dist = distance_between_points(data_xyz(1,nd),node_xyz(1,np_temp))
                  IF(dist.LT.min_dist)THEN
                     nd_min = nd
                     min_dist = dist
                  ENDIF !DIST
               ENDIF
            ENDDO
            my_closest(N) = nd_min
            map_array(nd_min) = ne_min
         ENDDO
      ENDIF
  
      DO nd = 1,num_data            ! for all seed/data points
         IF(map_array(nd).NE.0)THEN ! the data point is still in use
            IF(.NOT.inlist(nd,my_closest))THEN
               MIN_DIST=1.0e+10_dp     ! initialise the minimum (closest) distance
               DO noelem = 1,num_next_parents ! for each parent in the next branch generation
                  ne = local_parent(noelem)
                  np = elem_nodes(2,ne)
                  dist = distance_between_points(data_xyz(1,nd),node_xyz(1,np))
                  IF(DIST.LT.MIN_DIST)THEN
                     ne_min = ne
                     MIN_DIST=DIST
                  ENDIF
               ENDDO
               IF(first)THEN
                  map_array(nd)=ne_min
               ELSE
                  IF(MIN_DIST.LT.DISTANCE_LIMIT)THEN !keep seed points
                     map_array(nd)=ne_min
                  ELSE
                     map_array(nd)=0 !too far from branch ends, so discard
                  ENDIF
               ENDIF
            ENDIF
         ENDIF
      ENDDO
  
      num_seeds_from_elem(1:num_elems) = 0 !initialise the count of nd
      DO nd=1,num_data
         IF(map_array(nd).NE.0)THEN
            ne_min = map_array(nd)
            num_seeds_from_elem(ne_min) = num_seeds_from_elem(ne_min)+1
         ENDIF !map_array
      ENDDO !nd
  
  !!! If there is only 0 or 1 seed point grouped with an element then set it as a
  !!! terminal and remove a single seed point. Also involves modifying the local list of parents.
  
      IF(.NOT.first)THEN
         N_ELM_TEMP=num_next_parents
         DO N=1,num_next_parents
            ne_min=local_parent(N)
            IF(num_seeds_from_elem(ne_min).EQ.0)THEN !find closest point to end node
               nd_min = my_closest(N)
               map_array(nd_min)=0
               N_ELM_TEMP=N_ELM_TEMP-1
               local_parent(N)=0
               num_terminal=num_terminal+1
  
            ELSE IF(num_seeds_from_elem(ne_min).EQ.1)THEN
               DO nd=1,num_data
                  IF(map_array(nd).EQ.ne_min)THEN
                     map_array(nd)=0
                     local_parent(N)=0
                     N_ELM_TEMP=N_ELM_TEMP-1
                     num_terminal=num_terminal+1
                  ENDIF
               ENDDO !nd
  
            ENDIF !num_seeds_from_elem
         ENDDO !N
  
         DO N=1,num_next_parents
            IF(local_parent(N).EQ.0)THEN
               I=0
               DO WHILE((N+I.LT.num_next_parents).AND.(local_parent(N+I).EQ.0))
                  I=I+1
               ENDDO
               DO M=N,num_next_parents-I
                  local_parent(M)=local_parent(M+I)
               ENDDO !M
            ENDIF !local_parent
         ENDDO !N
         num_next_parents = N_ELM_TEMP
  
         CALL sort_integer_list(num_next_parents,local_parent)
  
      ENDIF
  
      DEALLOCATE(my_closest)
      IF(.NOT.first) DEALLOCATE(map_array_copy)
  
      CALL enter_exit(sub_name,2)
  
    END SUBROUTINE group_seeds_with_branch
  
  
    !###############################################################
    !
    !*grow_tree:* the main growing subroutine (public). Genertes a volume-filling
    ! tree into a closed surface.
    !
    SUBROUTINE grow_tree(parent_ne,surface_elems,angle_max,angle_min,&
         branch_fraction,length_limit,shortest_length,rotation_limit,to_export,filename)
      !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_GROW_TREE" :: GROW_TREE
  
      USE diagnostics, ONLY: enter_exit
      IMPLICIT NONE
  
      INTEGER,INTENT(in)  :: parent_ne                ! list of end branch elements to grow from
      INTEGER,INTENT(in)  :: surface_elems(:)         ! list of surface elements defining the host region
      REAL(dp),INTENT(in) :: angle_max                ! maximum branch angle with parent; in degrees
      REAL(dp),INTENT(in) :: angle_min                ! minimum branch angle with parent; in degrees
      REAL(dp),INTENT(in) :: branch_fraction          ! fraction of distance (to COFM) to branch
      REAL(dp),INTENT(in) :: length_limit             ! minimum length of a generated branch (shorter == terminal)
      REAL(dp),INTENT(in) :: shortest_length          ! length that short branches are reset to (shortest in model)
      REAL(dp),INTENT(in) :: rotation_limit           ! maximum angle of rotation of branching plane
      LOGICAL,INTENT(in) :: to_export                 ! option to export terminal element mapping to datapoints
      CHARACTER(len=*),INTENT(in) :: filename
  
      !Local variables
      INTEGER,ALLOCATABLE :: local_parent(:)          ! stores current generation of local parent elements
      INTEGER,ALLOCATABLE :: local_parent_temp(:)     ! temporary storage of next generation of local parent elems
      INTEGER,ALLOCATABLE :: map_seed_to_elem(:)      ! records current elem associated w. data points
      INTEGER,ALLOCATABLE :: map_seed_to_space(:)     ! records initial elem associated w. data points (the 'space')
      INTEGER,ALLOCATABLE :: num_seeds_from_elem(:)   ! records # of seeds currently grouped with an elem
      INTEGER,ALLOCATABLE :: triangle(:,:)
      CHARACTER(len=100) :: writefile
  
      INTEGER :: i,j,kount,M,N,nd,nd_min,ne,ne_grnd_parent,ne_parent,ne_stem,&
           noelem_parent,np,np_start,np_prnt_start,np_grnd_start,num_seeds_in_space,num_next_parents, &
           num_parents,num_triangles,num_vertices,num_elems_new,num_nodes_new,num_terminal
  
      REAL(dp),ALLOCATABLE :: vertex_xyz(:,:)
      REAL(dp),DIMENSION(3) :: COFM,candidate_xyz
      REAL(dp) :: distance_limit = 300.0_dp,length_parent
  
      LOGICAL :: make_branch,enough_points,first_group,internal, &
           limit_branching_angle = .TRUE., &  ! option to restrict branch angle
           limit_branching_plane = .FALSE.    ! option to restrict angle between branching planes
  
      CHARACTER(len=60) :: sub_name
  
      sub_name = 'grow_tree'
      CALL enter_exit(sub_name,1)
  
  
      IF(to_export)THEN
  !!! export vertices as nodes
         writefile = TRIM(filename)//'.txt'
         OPEN(40, file = writefile, status='replace')
         WRITE(40,'('' Data point number          Terminal element number'')')
      ENDIF
  
  
  
      CALL triangles_from_surface(num_triangles,num_vertices,surface_elems,triangle,vertex_xyz)
  
  !!! We can estimate the number of elements in the generated model based on the
  !!! number of data (seed) points. i.e. N = 2*N_data - 1. So the total number of
  !!! elements following tree generation will be ~ num_elems + 2*num_data. Use this estimate
  !!! to increase the node and element arrays.
      num_elems_new = num_elems + 2*num_data + 100
      num_nodes_new = num_nodes + 2*num_data + 100
      CALL reallocate_node_elem_arrays(num_elems_new,num_nodes_new)
  
  !!! Allocate memory for temporary arrays (need a more intelligent way of estimating size!)
      ALLOCATE(local_parent_temp(num_elems_new))
      ALLOCATE(local_parent(num_elems_new))
      ALLOCATE(num_seeds_from_elem(num_elems_new))
      ALLOCATE(map_seed_to_elem(num_data))
      ALLOCATE(map_seed_to_space(num_data))
  
  !!! Initialise local_parent to the list of parent elements, and num_parents (current
  !!! number of parent branches) to the number of parent branches.
      local_parent(1:SIZE(parentlist)) = parentlist(1:SIZE(parentlist))
      num_parents = COUNT(parentlist.NE.0) !initial number of 'terminal' parent branches
  
      NUM_SEEDS_FROM_ELEM = 0
      num_next_parents = num_parents
  
  
  !!! Calculate the initial grouping of data points with terminal elements
  !!! this defines the 'space' with which each seed is associated
  