anchor.diff:      real(c_double) :: half_width
   
       ! Allocate the new nodes and call sub-divide.
  -    allocate(next_candidates(num_next_candidates - 3)%nodes(num_nodes, 2))
  -    allocate(next_candidates(num_next_candidates - 2)%nodes(num_nodes, 2))
  -    allocate(next_candidates(num_next_candidates - 1)%nodes(num_nodes, 2))
  -    allocate(next_candidates(num_next_candidates)%nodes(num_nodes, 2))
  +    ! NOTE: This **assumes** but does not check that if the nodes are
  +    !       allocated, they are also the correct shape.
  +    if (.NOT. allocated(next_candidates(num_next_candidates - 3)%nodes)) then
  +       allocate(next_candidates(num_next_candidates - 3)%nodes(num_nodes, 2))
  +    end if
  +    if (.NOT. allocated(next_candidates(num_next_candidates - 2)%nodes)) then
  +       allocate(next_candidates(num_next_candidates - 2)%nodes(num_nodes, 2))
  +    end if
  +    if (.NOT. allocated(next_candidates(num_next_candidates - 1)%nodes)) then
  +       allocate(next_candidates(num_next_candidates - 1)%nodes(num_nodes, 2))
  +    end if
  +    if (.NOT. allocated(next_candidates(num_next_candidates)%nodes)) then
  +       allocate(next_candidates(num_next_candidates)%nodes(num_nodes, 2))
  +    end if
   
       call subdivide_nodes( &
            num_nodes, 2, candidate%nodes, degree, &
  
positive.code:   implicit none
    private &
         LocateCandidate, MAX_LOCATE_SUBDIVISIONS, LOCATE_EPS, &
         newton_refine_solve, split_candidate
    public newton_refine, locate_point
  
    ! For ``locate_point``.
  
negative.code:     real(c_double), intent(out) :: s_val, t_val
      ! Variables outside of signature.
      real(c_double) :: point(1, 2)
      integer(c_int) :: sub_index
      integer(c_int) :: num_candidates, num_next_candidates
      type(LocateCandidate), allocatable :: candidates_odd(:), candidates_even(:)
      real(c_double) :: s_approx, t_approx
      real(c_double) :: actual(1, 2)
      logical(c_bool) :: is_even
  
      point(1, :) = [x_val, y_val]
      ! Start out with the full curve.
      allocate(candidates_odd(1))  ! First iteration is odd.
      candidates_odd(1) = LocateCandidate(1.0_dp, 1.0_dp, 1.0_dp, nodes)
      num_candidates = 1
      s_val = LOCATE_MISS
  
      is_even = .TRUE.  ! At zero.
      do sub_index = 1, MAX_LOCATE_SUBDIVISIONS + 1
         is_even = .NOT. is_even  ! Switch parity.
         num_next_candidates = 0
  
         if (is_even) then
            ! Read from even, write to odd.
            call update_candidates( &
                 num_nodes, degree, point, num_candidates, candidates_even, &
                 num_next_candidates, candidates_odd)
         else
            ! Read from odd, write to even.
            call update_candidates( &
                 num_nodes, degree, point, num_candidates, candidates_odd, &
                 num_next_candidates, candidates_even)
         end if
  
         ! If there are no more candidates, we are done.
         if (num_next_candidates == 0) then
            return
         end if
  
         num_candidates = num_next_candidates
  
      end do
  