subroutine test_real_3D_threads(input, output, array_ref)
    implicit none
    double precision ::  input(N, N, N)
    double complex :: output(N, N, N)
    double complex :: array_ref(N, N, N)
    integer*8 :: plan
    integer :: i
    integer :: id, nb_threads
    integer, external :: omp_get_num_threads, omp_get_max_threads

!$OMP PARALLEL SHARED(nb_threads)
    nb_threads = omp_get_num_threads()
    print *, nb_threads
!$OMP END PARALLEL

    call dfftw_plan_with_nthreads(omp_get_max_threads())
    !call dfftw_plan_with_nthreads(3)
    call dfftw_plan_dft_r2c_3d(plan, N,N,N, input, output, FFTW_ESTIMATE)
    id = omp_get_thread_num()
    call init(input)
    ! Don't forget to initialize
    output(:,:,:) = 0

    call dfftw_execute_dft_r2c(plan, input, output)

    call compare(output, array_ref)
    call dfftw_destroy_plan(plan)

  end subroutine