sorting_doubles_t Subroutine

    subroutine sorting_doubles_t(doubles)
        type(doubles_t), intent(inout) :: doubles(:)
    

        ! merge_sort
        block
            type(doubles_t), dimension(:), allocatable :: tmp
            integer :: current_size, left, mid, right
            integer, parameter :: along = 1, &
                bubblesort_size = 20

            associate(X => doubles)
                ! Determine good number of starting splits
                block
                    integer :: n_splits
                    n_splits = 1
                    do while (size(X, along) / n_splits + merge(0, 1, mod(size(X, along), n_splits) == 0) > bubblesort_size)
                        n_splits = n_splits + 1
                    end do
                    current_size = size(X, along) / n_splits + merge(0, 1, mod(size(X, along), n_splits) == 0)
                end block

                ! Reuse this variable as tmp for swap in bubble_sort
                ! and for merge in merge_sort.
                block
                    integer :: max_buffer_size, n_merges
                    n_merges = ceiling(log(real(size(X, along)) / real(current_size)) / log(2.0))
                    max_buffer_size = current_size * merge(2**(n_merges - 1), 1, n_merges >= 1)
                    allocate(tmp(max_buffer_size))
                end block

                ! Bubble sort bins of size `bubblesort_size`.
                do left = lbound(X, along), ubound(X, along), current_size
                    right = min(left + bubblesort_size - 1, ubound(X, along))
                        ! bubblesort
    block
      integer :: n, i
      associate(V => X(left : right))
        do n = ubound(V, 1), lbound(V, 1) + 1, -1
          do i = lbound(V, 1), ubound(V, 1) - 1
            if (.not.   doubles_comp(V(i), V(i + 1))) then
                  ! swap
    block
        associate(tmp => tmp(1))
            tmp = V(i)
            V(i) = V(i + 1)
            V(i + 1) = tmp
        end associate
    end block
            end if
          end do
        end do
      end associate
    end block
                end do

                do while (current_size < size(X, along))
                    do left = lbound(X, along), ubound(X, along), 2 * current_size
                        right = min(left + 2 * current_size - 1, ubound(X, along))
                        mid = min(left + current_size, right) - 1
                        tmp(: mid - left + 1) = X(left : mid)
                            ! merge
    block
        integer :: i, j, k

        associate(A => tmp(: mid - left + 1), B => X(mid + 1 : right), C => X(left : right))

            if (size(A) + size(B) > size(C)) then
                error stop
            end if

            i = lbound(A, 1)
            j = lbound(B, 1)
            do k = lbound(C, 1), ubound(C, 1)
                if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then
                    if (  doubles_comp(A(i), B(j))) then
                        C(k) = A(i)
                        i = i + 1
                    else
                        C(k) = B(j)
                        j = j + 1
                    end if
                else if (i <= ubound(A, 1)) then
                    C(k) = A(i)
                    i = i + 1
                else if (j <= ubound(B, 1)) then
                    C(k) = B(j)
                    j = j + 1
                end if
            end do
        end associate
    end block
                    end do
                    current_size = 2 * current_size
                end do
            end associate
        end block
    contains
        logical elemental function doubles_comp(p1, p2)
            type(doubles_t), intent(in) :: p1, p2
            associate(idx_1 => [p1%j, p1%i, p1%b, p1%a], &
                      idx_2 => [p2%j, p2%i, p2%b, p2%a])
                doubles_comp = lex_leq(idx_1, idx_2)
            end associate
        end function
    end subroutine sorting_doubles_t