anchor.diff:    ! 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
  +          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
   
   
     !###############################################################
  
positive.code:     v(1:3) = elem_direction(1:3,ne-1)       ! new branch direction
      length = elem_field(ne_length,ne-1)
      nrml = cross_product(u,v) ! calculate normal
      IF(ABS(nrml(1))+ABS(nrml(2))+ABS(nrml(3)).LT.1.0e-4)THEN !co-linear branches
         ! check whether the other branch is co-linear with parent
         v(1:3) = elem_direction(1:3,ne)
         nrml = cross_product(u,v) ! calculate normal
         IF(ABS(nrml(1))+ABS(nrml(2))+ABS(nrml(3)).LT.1.0e-4)THEN !co-linear points
            WRITE(*,*) 'warning: both branches co-linear with parent'
            READ(*,*)
         ELSE
            elem_check = ne
            node_check = np
            elem_in_plane = ne-1
            node_in_plane = np-1
         ENDIF
      ENDIF
      CALL adjust_branch_angle(1,elem_check,np_start,np_prnt_start,node_check, &
           angle_max,angle_min)
  
  
  
negative.code:        Q(3,2)=2*(Q3*Q2+Q0*Q1)
         Q(3,3)=Q0**2-Q1**2-Q2**2+Q3**2
  
  !       X(1:3) = elem_direction(1:3,np3) ! unit vector
         X(1:3) = elem_direction(1:3,ne) ! unit vector
         length = distance_between_points(node_xyz(1,np3),node_xyz(1,np2))
  
  