anchor.diff:    !*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,&
  +  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
  +    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,&
  +    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
  +    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)
  +    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
  
positive.code:     !...angle between branching planes
      ANGLE=angle_btwn_vectors(NRML,NRML_PARENT)
      ANGLE_BETWEEN=ANGLE
      IF(ANGLE_BETWEEN.GT.PI/2.0_dp)THEN
         ANGLE_BETWEEN=ANGLE_BETWEEN-PI
         ANGLE=ANGLE_BETWEEN
      ENDIF
      !      ANGLE=PI/2.0_dp-ANGLE
      ANGLE=PI/2.0_dp+ANGLE
  
      IF(ABS(ANGLE_BETWEEN).GT.ROTATION_LIMIT.AND.ABS(ANGLE_BETWEEN) &
           .LT.PI/2.0_dp-ROTATION_LIMIT)THEN
         IF(ANGLE.LT.0.0_dp)THEN
            ROT_ANGLE=-(ANGLE+ROTATION_LIMIT)
         ELSE
            ROT_ANGLE=-(ANGLE-ROTATION_LIMIT)
         ENDIF
  
         ANGLE0=ANGLE
  
  
negative.code: !!! 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
  