IncrementIndex Subroutine

    subroutine IncrementIndex(vec, rank, limit, tReachedLimit)
        implicit none
        integer, intent(inout) :: vec(:)
        integer, intent(in) :: rank, limit
        logical, intent(inout) :: tReachedLimit
        integer :: i

        tReachedLimit = .false.

        ! start from the left, increment any element which differs in value by at least 2
        ! with its right-hand neighbour

        ! if we are starting from empty indices:
        if (vec(1) == 0) then
            call resetIndices(vec, rank)
            return
        end if

        do i = 1, rank
            if (i == rank) then
                if (vec(rank) < limit) then
                    vec(rank) = vec(rank) + 1
                    call ResetIndices(vec, rank - 1)
                    return
                else
                    tReachedLimit = .true.
                    return
                end if
            else if (vec(i + 1) - vec(i) > 1) then
                vec(i) = vec(i) + 1
                call ResetIndices(vec, i - 1)
                return
            end if
        end do
    end subroutine IncrementIndex