anchor.diff:                  np2 = elem_nodes(2,ne2)
                   new_xyz(:) = node_xyz(:,np0)*0.5_dp + node_xyz(:,np1)*0.25_dp + node_xyz(:,np2)*0.25_dp
                   node_xyz(:,np) = new_xyz(:)
  -             endif
  -          endif
  -       enddo
  -    enddo
  -    
  -    call enter_exit(sub_name,2)
  -
  -  end subroutine smooth_1d_tree
  +             ENDIF
  +          ENDIF
  +       ENDDO
  +    ENDDO
  +
  +    CALL enter_exit(sub_name,2)
  +
  +  END SUBROUTINE smooth_1d_tree
   
   
     !###############################################################
     !
  -  !*split_seed_points:* divides a set of seed points into two subsets 
  +  !*split_seed_points:* divides a set of seed points into two subsets
     ! using the plane that contains the parent branch and the seed point
     ! centre of mass. Decides which side of a plane a seed point is on by calculating
     ! the distance between two parallel planes: one which is defined by
     ! the parent and grandparent branch, and the other which contains a
     ! seed point.
     !
  -  subroutine split_seed_points(map_seed_to_elem,ne1,ne,np1,np2,np3,COFM,enough_points)
  +  SUBROUTINE split_seed_points(map_seed_to_elem,ne1,ne,np1,np2,np3,COFM,enough_points)
       !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_SPLIT_SEED_POINTS" :: SPLIT_SEED_POINTS
  -    
  -    use diagnostics, only: enter_exit
  -    implicit none
  -    
  -    integer :: map_seed_to_elem(*),ne,ne1,np1,np2,np3
  -    real(dp) :: COFM(3)
  -    logical :: enough_points
  +
  +    USE diagnostics, ONLY: enter_exit
  +    IMPLICIT NONE
  +
  +    INTEGER :: map_seed_to_elem(*),ne,ne1,np1,np2,np3
  +    REAL(dp) :: COFM(3)
  +    LOGICAL :: enough_points
   
       !Local variables
  -    integer :: DAT1,DAT2,nd,ND1_1ST,ND2_1ST,nsp,NPOINTS
  -    real(dp) :: DIST,NORML(4),P(3),Q(3),R(3)
  -    logical :: COLINEAR
  -    
  -    character(len=60) :: sub_name
  +    INTEGER :: DAT1,DAT2,nd,ND1_1ST,ND2_1ST,nsp,NPOINTS
  +    REAL(dp) :: DIST,NORML(4),P(3),Q(3),R(3)
  +    LOGICAL :: COLINEAR
  +
  +    CHARACTER(len=60) :: sub_name
   
       sub_name = 'split_seed_points'
  -    call enter_exit(sub_name,1)
  -    
  +    CALL enter_exit(sub_name,1)
  +
       R = COFM ! split based on cofm and branch
       P(1:3) = node_xyz(1:3,np2) ! point at start of parent branch
       Q(1:3) = node_xyz(1:3,np1) ! point at end of parent branch
   
  -!!! check whether the centre of mass and the parent start & end branches 
  -!!! are co-linear. if so, will need to use 'aunt' branch for split    
  +!!! check whether the centre of mass and the parent start & end branches
  +!!! are co-linear. if so, will need to use 'aunt' branch for split
       colinear = check_colinear_points(P,Q,R)
  -    if(colinear) R(1:3) = node_xyz(1:3,np3) !split based on parent and aunt
  -    call make_plane_from_3points(NORML,1,P,Q,R) !calculate plane
  -    
  +    IF(colinear) R(1:3) = node_xyz(1:3,np3) !split based on parent and aunt
  +    CALL make_plane_from_3points(NORML,1,P,Q,R) !calculate plane
  +
       NPOINTS=0
       DAT1=0
       DAT2=0
       ND1_1ST=0
       ND2_1ST=0
  -    do nd=1,num_data
  +    DO nd=1,num_data
          nsp=map_seed_to_elem(nd) !space # that random point belongs to
  -       if(nsp.eq.ne1)then !random point belongs to this element space
  +       IF(nsp.EQ.ne1)THEN !random point belongs to this element space
             NPOINTS=NPOINTS+1
             dist = -scalar_product_3(norml,data_xyz(1,nd)) - norml(4) ! distance between two planes
  -          if(dist.ge.0.0_dp)then
  -             if(dat1.eq.0) nd1_1st = nd
  +          IF(dist.GE.0.0_dp)THEN
  +             IF(dat1.EQ.0) nd1_1st = nd
                DAT1=DAT1+1
                map_seed_to_elem(nd)=ne+1
  -          else if(DIST.lt.0.0_dp)then
  -             if(DAT2.eq.0) ND2_1ST=nd
  +          ELSE IF(DIST.LT.0.0_dp)THEN
  +             IF(DAT2.EQ.0) ND2_1ST=nd
                DAT2=DAT2+1
                map_seed_to_elem(nd)=ne+2
  -          endif
  -       endif
  -    enddo !nd
  -    
  -    if(dat1.eq.0.and.dat2.eq.0)then
  -       enough_points = .false.
  -       write(*,'('' Zero seed points associated with parent'',I6)') ne1
  -       read(*,*)
  -    else
  -       if(dat1.eq.0)then
  +          ENDIF
  +       ENDIF
  +    ENDDO !nd
  +
  +    IF(dat1.EQ.0.AND.dat2.EQ.0)THEN
  +       enough_points = .FALSE.
  +       WRITE(*,'('' Zero seed points associated with parent'',I6)') ne1
  +       READ(*,*)
  +    ELSE
  +       IF(dat1.EQ.0)THEN
             map_seed_to_elem(nd2_1st) = ne+1
             dat1 = dat1+1
             dat2 = dat2-1
  -       elseif(dat2.eq.0)then
  +       ELSEIF(dat2.EQ.0)THEN
             map_seed_to_elem(nd1_1st) = ne+2
             dat2 = dat2+1
             dat1 = dat1-1
  -       endif
  -       
  -       enough_points = .true.    
  -    endif
  -
  -    call enter_exit(sub_name,2)
  -
  -  end subroutine split_seed_points
  -  
  -  
  +       ENDIF
  +
  +       enough_points = .TRUE.
  +    ENDIF
  +
  +    CALL enter_exit(sub_name,2)
  +
  +  END SUBROUTINE split_seed_points
  +
  +
     !###############################################################
     !
  -  !*split_seed_points_initial:* divides a set of seed points into N subsets 
  +  !*split_seed_points_initial:* divides a set of seed points into N subsets
     ! to match N terminal branches, using the plane that is orthogonal to the branching plane
     ! of child branches, and that passes mid-way between child branches.
     !
  -  subroutine split_seed_points_initial(map_array,ne_stem)
  +  SUBROUTINE split_seed_points_initial(map_array,ne_stem)
       !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_SPLIT_SEED_POINTS_INITIAL" :: SPLIT_SEED_POINTS_INITIAL
   
  -    use diagnostics, only: enter_exit
  -    implicit none
  -    
  -    integer :: map_array(:),ne_stem
  +    USE diagnostics, ONLY: enter_exit
  +    IMPLICIT NONE
  +
  +    INTEGER :: map_array(:),ne_stem
   
       !Local variables
  -    integer :: local_parent(20),local_parent_temp(20),M,nd,ne_parent, &
  +    INTEGER :: local_parent(20),local_parent_temp(20),M,nd,ne_parent, &
            ne1,ne2,np0,np1,np2,nsp,num_points,num_next_parents,num_parents
  -    real(dp) :: DIST,dist_p1,dist_p2,NORML(4),P(3),Q(3),R(3)
  -    
  -    character(len=60) :: sub_name
  +    REAL(dp) :: DIST,dist_p1,dist_p2,NORML(4),P(3),Q(3),R(3)
  +
  +    CHARACTER(len=60) :: sub_name
   
       sub_name = 'split_seed_points_initial'
  -    call enter_exit(sub_name,1)
  +    CALL enter_exit(sub_name,1)
   
       ne_parent = ne_stem
  -    do while(elem_cnct(1,0,ne_parent).eq.1)
  +    DO WHILE(elem_cnct(1,0,ne_parent).EQ.1)
          ne_parent = elem_cnct(1,1,ne_parent) ! get the next element in a refined branch
  -    enddo
  +    ENDDO
       map_array(:) = ne_parent ! initialise that all seed points map to the stem branch
   
       num_next_parents = 1
       local_parent(1) = ne_parent
  -    
  -    do while(num_next_parents.ne.0) !while still some parent branches with seed points
  +
  +    DO WHILE(num_next_parents.NE.0) !while still some parent branches with seed points
          num_parents = num_next_parents ! update the number of current local parent branches
          num_next_parents = 0 ! reset the number of local parent branches in next generation
  -       do M = 1,num_parents ! for each of the current local parent branches
  +       DO M = 1,num_parents ! for each of the current local parent branches
             ne_parent = local_parent(M) !parent element #
  -          do while(elem_cnct(1,0,ne_parent).eq.1)
  +          DO WHILE(elem_cnct(1,0,ne_parent).EQ.1)
                ne_parent = elem_cnct(1,1,ne_parent) ! get the next element in a refined branch
  -          enddo
  +          ENDDO
             np0 = elem_nodes(2,ne_parent)
             ne1 = elem_cnct(1,1,ne_parent)
  -          do while(elem_cnct(1,0,ne1).eq.1)
  +          DO WHILE(elem_cnct(1,0,ne1).EQ.1)
                ne1 = elem_cnct(1,1,ne1) ! get the next element in a refined branch
  -          enddo
  +          ENDDO
             ne2 = elem_cnct(1,2,ne_parent)
  -          do while(elem_cnct(1,0,ne2).eq.1)
  +          DO WHILE(elem_cnct(1,0,ne2).EQ.1)
                ne2 = elem_cnct(1,1,ne2) ! get the next element in a refined branch
  -          enddo
  +          ENDDO
             np1 = elem_nodes(2,ne1)
             np2 = elem_nodes(2,ne2)
             P(:) = node_xyz(:,np0) ! point at end of parent branch
             Q(:) = node_xyz(:,np1) ! point at end of child1 branch
             R(:) = node_xyz(:,np2) ! point at end of child2 branch
  -          call make_plane_from_3points(NORML,1,P,Q,R) !calculate plane
  -          
  +          CALL make_plane_from_3points(NORML,1,P,Q,R) !calculate plane
  +
             P(:) = node_xyz(:,np0) ! point at end of parent branch
             Q(:) = 0.5_dp*(node_xyz(:,np1)+node_xyz(:,np2))
             R(1:3) = Q(1:3) + NORML(1:3)
  -          call make_plane_from_3points(NORML,1,P,Q,R) !calculate plane
  -!!! NORML is now the plane between the child branches    
  -          
  +          CALL make_plane_from_3points(NORML,1,P,Q,R) !calculate plane
  +!!! NORML is now the plane between the child branches
  +
             dist_p1 = -scalar_product_3(norml,node_xyz(:,np1)) - norml(4) ! distance between two planes
             dist_p2 = -scalar_product_3(norml,node_xyz(:,np2)) - norml(4) ! distance between two planes
   
  -          do nd = 1,num_data
  +          DO nd = 1,num_data
                nsp = map_array(nd) !space # that random point belongs to
  -             if(nsp.eq.ne_parent)then !random point belongs to this element space
  +             IF(nsp.EQ.ne_parent)THEN !random point belongs to this element space
                   dist = -scalar_product_3(norml,data_xyz(1,nd)) - norml(4) ! distance between two planes
  -                if(dist.ge.0.0_dp.and.dist_p1.ge.0.0_dp)then
  +                IF(dist.GE.0.0_dp.AND.dist_p1.GE.0.0_dp)THEN
                      map_array(nd) = ne1
  -                else if(dist.ge.0.0_dp.and.dist_p1.lt.0.0_dp)then
  +                ELSE IF(dist.GE.0.0_dp.AND.dist_p1.LT.0.0_dp)THEN
                      map_array(nd) = ne2
  -                else if(dist.le.0.0_dp.and.dist_p2.le.0.0_dp)then
  +                ELSE IF(dist.LE.0.0_dp.AND.dist_p2.LE.0.0_dp)THEN
                      map_array(nd) = ne2
  -                else if(dist.le.0.0_dp.and.dist_p2.gt.0.0_dp)then
  +                ELSE IF(dist.LE.0.0_dp.AND.dist_p2.GT.0.0_dp)THEN
                      map_array(nd) = ne1
  -                endif
  -             endif
  -          enddo !nd
  -
  -          num_points = count(map_array.eq.ne1)
  -          if(num_points.eq.0)then
  -             write(*,'('' Warning: number of points for element'',i6,'' is zero'')') ne1
  -             write(*,'('' Press enter to continue; however the code is likely to fail'')')
  -          endif
  -          num_points = count(map_array.eq.ne2)
  -          if(num_points.eq.0)then
  -             write(*,'('' Warning: number of points for element'',i6,'' is zero'')') ne2
  -             write(*,'('' Press enter to continue; however the code is likely to fail'')')
  -          endif
  -       
  -          if(elem_cnct(1,0,ne1).ne.0)then
  +                ENDIF
  +             ENDIF
  +          ENDDO !nd
  +
  +          num_points = COUNT(map_array.EQ.ne1)
  +          IF(num_points.EQ.0)THEN
  +             WRITE(*,'('' Warning: number of points for element'',i6,'' is zero'')') ne1
  +             WRITE(*,'('' Press enter to continue; however the code is likely to fail'')')
  +          ENDIF
  +          num_points = COUNT(map_array.EQ.ne2)
  +          IF(num_points.EQ.0)THEN
  +             WRITE(*,'('' Warning: number of points for element'',i6,'' is zero'')') ne2
  +             WRITE(*,'('' Press enter to continue; however the code is likely to fail'')')
  +          ENDIF
  +
  +          IF(elem_cnct(1,0,ne1).NE.0)THEN
                num_next_parents = num_next_parents+1
                local_parent_temp(num_next_parents) = ne1
  -          endif
  -          if(elem_cnct(1,0,ne2).ne.0)then
  +          ENDIF
  +          IF(elem_cnct(1,0,ne2).NE.0)THEN
                num_next_parents = num_next_parents+1
                local_parent_temp(num_next_parents) = ne2
  -          endif
  -       enddo ! num_parents
  +          ENDIF
  +       ENDDO ! num_parents
          local_parent(1:num_next_parents) = local_parent_temp(1:num_next_parents)
  -    enddo 
  -
  -    call enter_exit(sub_name,2)
  -
  -  end subroutine split_seed_points_initial
  -  
  -  
  +    ENDDO
  +
  +    CALL enter_exit(sub_name,2)
  +
  +  END SUBROUTINE split_seed_points_initial
  +
  +
     !###############################################################
     !
     !*closest_seed_to_node:* returns the closest seed point to a given branch node
     !
  -  function closest_seed_to_node(map_seed_to_elem,np)
  +  FUNCTION closest_seed_to_node(map_seed_to_elem,np)
       !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_CLOSEST_SEED_TO_NODE" :: CLOSEST_SEED_TO_NODE
   
  -    use diagnostics, only: enter_exit
  -    implicit none
  -    
  -    integer,intent(in) :: map_seed_to_elem(*),np
  -    
  +    USE diagnostics, ONLY: enter_exit
  +    IMPLICIT NONE
  +
  +    INTEGER,INTENT(in) :: map_seed_to_elem(*),np
  +
       !Local variables
  -    integer :: nd
  -    integer :: closest_seed_to_node
  -    real(dp) :: distance,min_distance
  -    
  +    INTEGER :: nd
  +    INTEGER :: closest_seed_to_node
  +    REAL(dp) :: distance,min_distance
  +
       min_distance = 1.0e+10_dp
  -    do nd=1,num_data
  -       if(map_seed_to_elem(nd).ne.0)then
  +    DO nd=1,num_data
  +       IF(map_seed_to_elem(nd).NE.0)THEN
             distance = distance_between_points(data_xyz(1,nd),node_xyz(1,np))
  -          if(distance.lt.min_distance)then
  +          IF(distance.LT.min_distance)THEN
                closest_seed_to_node = nd
                min_distance = distance
  -          endif !DIST
  -       endif !map_seed_to_elem
  -    enddo !nd
  -    
  -  end function closest_seed_to_node
  -  
  +          ENDIF !DIST
  +       ENDIF !map_seed_to_elem
  +    ENDDO !nd
  +
  +  END FUNCTION closest_seed_to_node
  +
   
     !###############################################################
     !
     !*closest_seed_to_node_in_group:* finds the closest seed point to a node
     ! that is in a group of seed points currently associated with a specific element
     !
  -  function closest_seed_to_node_in_group(map_seed_to_elem,ne,np)
  +  FUNCTION closest_seed_to_node_in_group(map_seed_to_elem,ne,np)
       !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_CLOSEST_SEED_TO_NODE_IN_GROUP" :: CLOSEST_SEED_TO_NODE_IN_GROUP
  -    
  -    integer,intent(in) :: map_seed_to_elem(*),ne,np
  -    
  +
  +    INTEGER,INTENT(in) :: map_seed_to_elem(*),ne,np
  +
       !Local variables
  -    integer :: nd
  -    integer :: closest_seed_to_node_in_group
  -    real(dp) :: distance,min_distance
  -    
  +    INTEGER :: nd
  +    INTEGER :: closest_seed_to_node_in_group
  +    REAL(dp) :: distance,min_distance
  +
       min_distance = 1.0e+10_dp
  -    do nd=1,num_data
  -       if(map_seed_to_elem(nd).eq.ne)then
  +    DO nd=1,num_data
  +       IF(map_seed_to_elem(nd).EQ.ne)THEN
             distance = distance_between_points(data_xyz(1,nd),node_xyz(1,np))
  -          if(distance.lt.min_distance)then
  +          IF(distance.LT.min_distance)THEN
                closest_seed_to_node_in_group = nd
                min_distance = distance
  -          endif !DIST
  -       endif !map_seed_to_elem
  -    enddo !nd
  -    
  -  end function closest_seed_to_node_in_group
  -  
  +          ENDIF !DIST
  +       ENDIF !map_seed_to_elem
  +    ENDDO !nd
  +
  +  END FUNCTION closest_seed_to_node_in_group
  +
   
     !###############################################################
     !
     !*rotation_angle:* calculates angle between two branching planes
     !
  -  function rotation_angle(np1,np2,np3,np4,np5)
  +  FUNCTION rotation_angle(np1,np2,np3,np4,np5)
       !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_ROTATION_ANGLE" :: ROTATION_ANGLE
  -    
  -    integer,intent(in) :: np1,np2,np3,np4,np5
  +
  +    INTEGER,INTENT(in) :: np1,np2,np3,np4,np5
   
       !Local variables
  -    real(dp) :: norm_1(4),norm_2(4)
  -    real(dp) :: rotation_angle
  -    
  -    call make_plane_from_3points(norm_1,2,node_xyz(1,np1),node_xyz(1,np2),node_xyz(1,np3))
  -    call make_plane_from_3points(norm_2,2,node_xyz(1,np2),node_xyz(1,np4),node_xyz(1,np5))
  +    REAL(dp) :: norm_1(4),norm_2(4)
  +    REAL(dp) :: rotation_angle
  +
  +    CALL make_plane_from_3points(norm_1,2,node_xyz(1,np1),node_xyz(1,np2),node_xyz(1,np3))
  +    CALL make_plane_from_3points(norm_2,2,node_xyz(1,np2),node_xyz(1,np4),node_xyz(1,np5))
       rotation_angle = angle_btwn_vectors(norm_1,norm_2)
  -    
  -  end function rotation_angle
  -  
  +
  +  END FUNCTION rotation_angle
  +
   
     !###############################################################
     !
     !*vector_for_angle_limit:* Calculates the new direction
     ! of an element, given current and target angles. Specifically, solves small
     ! system of equations to get new direction of branch (w), such that the branch
  -  ! remains in-plane (n.w = 0), the angle of w with u is defined 
  -  ! (as u.w = cos(angle_with_u)), and the angle of w with original direction (v) 
  +  ! remains in-plane (n.w = 0), the angle of w with u is defined
  +  ! (as u.w = cos(angle_with_u)), and the angle of w with original direction (v)
     ! is defined (as v.w = cos(angle_with_v)).
     !
  -  function vector_for_angle_limit(U,V,angle_with_u,angle_with_v)
  +  FUNCTION vector_for_angle_limit(U,V,angle_with_u,angle_with_v)
       !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_VECTOR_FOR_ANGLE_LIMIT" :: VECTOR_FOR_ANGLE_LIMIT
   
  -    real(dp),intent(in) :: U(*),V(*),angle_with_u,angle_with_v
  +    REAL(dp),INTENT(in) :: U(*),V(*),angle_with_u,angle_with_v
   
       !Local variables
  -    real(dp) :: A(3,3),N_UV(3),VECTOR(3),W(3)
  -    real(dp) :: vector_for_angle_limit(3)
  +    REAL(dp) :: A(3,3),N_UV(3),VECTOR(3),W(3)
  +    REAL(dp) :: vector_for_angle_limit(3)
   
       N_UV = cross_product(U,V)     ! calculate the normal to the vectors U and V
       N_UV = unit_vector(N_UV) ! unit vector for normal
       A(1,1:3) = N_UV(1:3)
       A(2,1:3) = U(1:3)
       A(3,1:3) = V(1:3)
  -       
  +
       VECTOR(1) = 0.0_dp
       VECTOR(2) = DCOS(angle_with_u)
       VECTOR(3) = DCOS(angle_with_v)
  -       
  +
       w = mesh_a_x_eq_b(A,VECTOR)
       vector_for_angle_limit = unit_vector(W)
   
  -  end function vector_for_angle_limit
  -  
  +  END FUNCTION vector_for_angle_limit
  +
   
     !###############################################################
  -  
  -end module growtree
  +
  +END MODULE growtree
   
   !###############################################################
  -  
  
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: !!! the closest parent end-point to each seed point. This has been replaced by splitting
  !!! seed points using the orthogonal to branching planes of the upper tree.
      map_seed_to_space(1:num_data) = parentlist(1) !#! this is done for the new-style growing (full grow per terminal)
      IF(num_parents.GT.1)THEN
         first_group = .TRUE.
         !       call group_seeds_with_branch(map_seed_to_space,num_next_parents,num_seeds_from_elem,&
         !            num_terminal,local_parent,500.0_dp,first_group)
         CALL split_seed_points_initial(map_seed_to_space,parent_ne)
      ENDIF !parentlist.gt.1
      first_group = .FALSE.
  
      WRITE(*,'(''  parent  #seeds  #terminal'')')
  
      ! Set initial values for local and global nodes and elements
      ne = num_elems !initialise mesh global element #
      np = num_nodes !initialise mesh global node #
  
  !!! loop over the initial parent list (the initial conditions/terminal elements for growing)
  !!! growing is done into 'spaces', where each 'space' is the initial grouping of seed
  !!! points with the closest terminal branch. This grouping can be manipulated to take into account
  !!! the size of the initial terminal branches. e.g. smaller diameter --> smaller set of seeds
  
  
      DO noelem_parent = 1,num_parents
         ne_stem = parentlist(noelem_parent) ! the 'stem' parent element for the 'space'
         map_seed_to_elem = 0 ! initialise the seed mapping array
         num_seeds_in_space = 0 !initialise the number of seed points in the 'space'
         DO nd = 1,num_data ! for all of the seed points (stored in data_xyz array)
            IF(map_seed_to_space(nd).EQ.ne_stem)THEN ! for the points in this space
               map_seed_to_elem(nd) = ne_stem ! record the current element associated with seed point nd
               num_seeds_in_space = num_seeds_in_space+1 ! count number of seed points in the space
            ENDIF
         ENDDO
  
         num_next_parents = 1 ! initialise the number of current local parent branches
         local_parent(1) = ne_stem ! first local parent branch is the 'stem' branch
         num_terminal = 0 ! initialise the number of definite terminal branches
  
  !!! bifurcating distributive algorithm
         DO WHILE(num_next_parents.NE.0) !while still some parent branches with seed points
            num_parents = num_next_parents ! update the number of current local parent branches
            num_next_parents = 0 ! reset the number of local parent branches in next generation
  
            DO M = 1,num_parents ! for each of the current local parent branches
               ne_parent = local_parent(M) !parent element #
               ! Calculate centre of mass of current seed point set
               CALL calculate_seed_cofm(map_seed_to_elem,ne_parent,COFM)
  
               ne_grnd_parent = elem_cnct(-1,1,ne_parent) !grandparent global element #
               np_start = elem_nodes(2,ne_parent) !parent global end node #
               np_prnt_start = elem_nodes(1,ne_parent) !parent global start node #
               np_grnd_start = elem_nodes(1,ne_grnd_parent) !grandparent global start node #
  
               length_parent = elem_field(ne_length,ne_parent)
  
  !!! Split each set of seed points using the plane defined by the
  !!! parent branch and the centre of mass. Seed points get associated with NEW elements (ne+1,ne+2)
               CALL split_seed_points(map_seed_to_elem,ne_parent,ne,np_start,&
                    np_prnt_start,np_grnd_start,COFM,enough_points)
  
  !!! check whether enough seed points remaining in BOTH seed groups for branching to be done
  !!! (note: this could be improved to continue branching in one set of seeds)
               IF(enough_points)THEN
                  DO N = 1,2 !for each of the two new branches
                     ! Set up arrays for new element and node
                     ! after create_new_node the current element == ne and current node == np
                     CALL create_new_node(ne,ne_parent,np,np_start,.TRUE.)
                     ! find the centre of mass of seed points
                     CALL calculate_seed_cofm(map_seed_to_elem,ne,COFM)
                     ! Generate a branch directed towards the centre of mass. Returns location
                     ! of end node in candidate_xyz (adjusted below based on length and shape criteria)
                     CALL branch_to_cofm(map_seed_to_elem,ne,np_start,&
                          COFM,branch_fraction,length_limit,length_parent,shortest_length,&
                          candidate_xyz,make_branch)
                     node_xyz(1:3,np) = candidate_xyz(1:3) ! the new node location is as returned by 'branch_to_cofm'
                     CALL calc_branch_direction(ne) ! calculate direction of the new branch
                     elem_field(ne_length,ne) = distance_between_points(node_xyz(1,np_start),node_xyz(1,np))
                     ! Check whether this is a new parent branch or a terminal branch
                     IF(make_branch.AND.enough_points)THEN ! meets all criteria for continuing branching
                        num_next_parents = num_next_parents+1 ! increment the number of next parents
                        local_parent_temp(num_next_parents) = ne !records the elements that are parents
                     ELSE ! this is a terminal branch
                        num_terminal = num_terminal+1 ! increment the number of terminal branches
                     ENDIF
                  ENDDO !N (for both new branches)
  
                  IF(limit_branching_angle)THEN
                     ! Check that the branch angles are not too large or too small
                     ! Correct such that branches stay in the original branching plane
                     CALL limit_branch_angles(ne,ne_parent,np,&
                          np_prnt_start,np_start,angle_max,angle_min)
  
                     internal = point_internal_to_surface(num_vertices,triangle,node_xyz(1:3,np),&
                          vertex_xyz)
                     IF(.NOT.internal)THEN ! halve the length, and make terminal
                        elem_field(ne_length,ne) = 0.5_dp*distance_between_points&
                             (node_xyz(1,np_start),node_xyz(1,np))
                        node_xyz(1:3,np) = node_xyz(1:3,np_start) + 0.5_dp*&
  