subroutine test_real_3D_threads(input, output, array_ref)
    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

    !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)
!$OMP PARALLEL SHARED(input, output) PRIVATE(plan, id)
    id = omp_get_thread_num()
    call init(input)
    ! Don't forget to initialize
!$OMP END PARALLEL
    output(:,:,:) = 0

    call dfftw_execute_dft_r2c(plan, input, output)

    print *, "thread ", id
    call compare(output, array_ref)
    call dfftw_destroy_plan(plan)

  end subroutine