anchor.diff:      node_xyz(1:3,np1) = node_xyz(1:3,np0)+elem_field(ne_length,ne)*elem_direction(1:3,ne)
   
       np0 = np1
  -    do i=1,elem_cnct(1,0,ne) ! for each child
  +    DO i=1,elem_cnct(1,0,ne) ! for each child
          ne1 = elem_cnct(1,i,ne)
          np1 = elem_nodes(2,ne1)
          elem_field(ne_length,ne1) = 0.8_dp*elem_field(ne_length,ne1)
          node_xyz(1:3,np1) = node_xyz(1:3,np0)+elem_field(ne_length,ne1)*elem_direction(1:3,ne1)
  -       do j = 1,elem_cnct(1,0,ne1)
  +       DO j = 1,elem_cnct(1,0,ne1)
             ne2 = elem_cnct(1,j,ne1)
             np2 = elem_nodes(2,ne2)
             elem_field(ne_length,ne2) = 0.8_dp*elem_field(ne_length,ne2)
             node_xyz(1:3,np2) = node_xyz(1:3,np1)+elem_field(ne_length,ne2)*elem_direction(1:3,ne2)
  -       enddo
  -    enddo
  -
  -    call enter_exit(sub_name,2)
  -
  -  end subroutine shorten_branch_and_children
  +       ENDDO
  +    ENDDO
  +
  +    CALL enter_exit(sub_name,2)
  +
  +  END SUBROUTINE shorten_branch_and_children
   
   
     !###############################################################
  
positive.code:                       internal = point_internal_to_surface(num_vertices,triangle,node_xyz(1:3,np),&
                             vertex_xyz)
                        kount = 0
                        do while(.not.internal) 
                           kount = kount+1
                           call shorten_branch_and_children(ne_parent)
                           internal = point_internal_to_surface(num_vertices,triangle,node_xyz(1:3,np),&
                                vertex_xyz)
                           if(kount.ge.3)then
                              call shorten_branch_and_children(elem_cnct(-1,1,ne_parent))
                              internal = point_internal_to_surface(num_vertices,triangle,node_xyz(1:3,np),&
                                   vertex_xyz)
                           endif
                           if(kount.gt.5.and.(.not.internal))then
                              write(*,'('' WARNING: element'',i6,'' not internal'')') ne
                              internal = .true.
                           endif
                        enddo
                        if(inlist(ne,local_parent_temp))then ! set to be terminal
                           local_parent_temp(num_next_parents) = 0
                           num_next_parents = num_next_parents-1 ! decrement the number of next parents
                           num_terminal = num_terminal+1
                           nd_min = closest_seed_to_node_in_group(map_seed_to_elem,ne,np) ! closest seed point                         
                           map_seed_to_elem(nd_min) = 0 ! remove seed point from list
                           map_seed_to_space(nd_min) = ne ! recording element number
  
                           if(to_export) then
                           write(40,*) nd_min,ne
                           endif
  
                        endif
                     endif !.not.internal
  
                     internal = point_internal_to_surface(num_vertices,triangle,node_xyz(1:3,np-1),&
                          vertex_xyz)
                     if(.not.internal)then ! halve the length, and halve the angle from parent
                        elem_field(ne_length,ne-1) = 0.5_dp*distance_between_points&
                             (node_xyz(1,np_start),node_xyz(1,np-1))
                        node_xyz(1:3,np-1) = node_xyz(1:3,np_start) + 0.5_dp*&
                             elem_field(ne_length,ne-1)*elem_direction(1:3,ne-1)
                        internal = point_internal_to_surface(num_vertices,triangle,node_xyz(1:3,np-1),&
                             vertex_xyz)
                        do while(.not.internal)
                           kount = kount+1
                           call shorten_branch_and_children(ne_parent)
                           internal = point_internal_to_surface(num_vertices,triangle,node_xyz(1:3,np),&
                                vertex_xyz)
                           if(kount.ge.3)then
                              call shorten_branch_and_children(elem_cnct(-1,1,ne_parent))
                              internal = point_internal_to_surface(num_vertices,triangle,node_xyz(1:3,np),&
                                   vertex_xyz)
                           endif
                           if(kount.gt.4.and.(.not.internal))then
                              write(*,'('' WARNING: element'',i6,'' not internal'')') ne-1
                              internal = .true.
                           endif
                        enddo
                        if(inlist(ne-1,local_parent_temp))then ! set to be terminal
                           do i = 1,num_next_parents
                              if(local_parent_temp(i).eq.ne-1)then
                                 do j = i,num_next_parents-1
                                    local_parent_temp(j) = local_parent_temp(i+1)
                                 enddo
                              endif
                           enddo
                           num_next_parents = num_next_parents-1 ! decrement the number of next parents
                           num_terminal = num_terminal+1
                           nd_min = closest_seed_to_node_in_group(map_seed_to_elem,ne-1,np-1) ! closest seed point
                           map_seed_to_elem(nd_min) = 0 ! remove seed point from list
                           map_seed_to_space(nd_min) = ne ! recording element number
  
                           if(to_export) then
                           write(40,*) nd_min,ne
                           endif
                        endif
                     endif ! .not.internal
  
                  endif
  
                  if(limit_branching_plane)then
                     ! Check the angle of rotation between child and parent branching planes.
                     ! If absolute angle is larger than a user-specified limit, adjust to be
                     ! the limit value. This is used for making sure that the CFD geometry turns out ok.
                     call check_branch_rotation_plane(map_seed_to_elem,ne,ne_grnd_parent,ne_parent, &
                          local_parent_temp,num_next_parents, &
                          np,np_start,np_prnt_start,np_grnd_start,num_terminal,rotation_limit)
                     internal = point_internal_to_surface(num_vertices,triangle,node_xyz(1:3,np),vertex_xyz)
                  endif
  
               else
                  write(*,*) 'terminal, not enough points',ne !!! never happens!!!
                  read(*,*)
                  ! Not enough seed points in the set during the split.
                  ! Find the closest seed point to node np_start, and remove from seeds
                  num_terminal = num_terminal+1 ! increment number of terminal branches
                  nd_min = closest_seed_to_node(map_seed_to_elem,np_start) ! closest seed point
                  map_seed_to_elem(nd_min) = 0 ! remove seed point from list
                  map_seed_to_space(nd_min) = ne ! record the element to data point mapping
               endif
            enddo ! for all current parent branches
            ! Copy the temporary list of branches to local_parent. These become the
            ! parent elements for the next branching   
            local_parent(1:num_next_parents) = local_parent_temp(1:num_next_parents)
            ! Regroup the seed points with the closest current parent 
            call group_seeds_with_branch(map_seed_to_elem,num_next_parents,num_seeds_from_elem,&
                 num_terminal,local_parent,DISTANCE_LIMIT,.FALSE.)
  
         enddo ! while still parent branches
  
         write(*,'(I7,I8,I9)') ne_parent,num_seeds_in_space,num_terminal
  
      enddo ! for each initial parent
      endif
      close(40)
     
      
  !!! set new total numbers of nodes and elements
      num_nodes=np !highest node # in nr
      num_elems=ne !highest element # in nr
      
  !!! update the tree connectivity
      call element_connectivity_1d
  !!! calculate branch generations and orders
      call evaluate_ordering
  !!! deallocate temporary arrays
      deallocate(vertex_xyz)
      deallocate(triangle)
      deallocate(local_parent_temp)
      deallocate(local_parent)
      deallocate(map_seed_to_elem)
      deallocate(map_seed_to_space)
      deallocate(num_seeds_from_elem)
      
      call enter_exit(sub_name,2)
  endif
    end subroutine grow_tree
  
    
    !###############################################################
    !
    !*limit_branch_angles:* checks both new branches at once to make sure the angles
    ! to parent branch not too large or too small. Both branches are checked to make sure
    ! they are not co-linear with parent. One branch is checked for its angle, and then 
    ! the other is made to be in-plane with the parent and sibling branch.
    !
    subroutine limit_branch_angles(ne,ne_parent,np,&
         np_prnt_start,np_start,angle_max,angle_min)
      !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_LIMIT_BRANCH_ANGLES" :: LIMIT_BRANCH_ANGLES
  
      use diagnostics, only: enter_exit
      implicit none
      
      integer,intent(in) :: ne,ne_parent,np,np_prnt_start,np_start
      real(dp) :: angle_max,angle_min
      
      !Local variables
      integer :: elem_check,elem_in_plane,node_check,node_in_plane
      real(dp),dimension(3) :: B,NRML,U,V,W
      real(dp),dimension(3,3) :: A
      real(dp) :: angle_uw,angle_uv,angle_vw,length,length_w
      
      character(len=60) :: sub_name
  
      sub_name = 'limit_branch_angles'
      call enter_exit(sub_name,1)
  
      elem_check = ne-1  ! the default element to check angle for is the '1st' branch
      node_check = np-1  ! node associated with ne-1
  
negative.code:   ! 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)
      !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
  
      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
                  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
                  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
                  ne1 = elem_cnct(1,1,ne)
                  ne2 = elem_cnct(1,2,ne)
                  np0 = elem_nodes(1,ne)
  