def METHOD_NAME(Nphi, Ntheta, Nr, k, dealias, dtype, basis, output_scales, parallel):
    # Basis
    c, d, b, phi, theta, r, x, y, z = basis(Nphi, Ntheta, Nr, k, dealias, dtype)
    # Fields
    u = field.Field(name='u', dist=d, bases=(b,), dtype=dtype)
    u.preset_scales(b.domain.dealias)
    u['g'] = np.sin(x) * np.sin(y) * np.sin(z)
    # Problem
    dt = operators.TimeDerivative
    problem = problems.IVP([u])
    problem.add_equation((dt(u), 0))
    # Solver
    solver = solvers.InitialValueSolver(problem, timesteppers.RK222, matrix_coupling=[False, False, True])
    # Output
    tasks = [u(phi=np.pi), u(theta=np.pi/2), u(r=1.0), u]
    with tempfile.TemporaryDirectory(dir='.') as tempdir:
        tempdir = pathlib.Path(tempdir).stem
        output = solver.evaluator.add_file_handler(tempdir, iter=1, parallel=parallel)
        for task in tasks:
            output.add_task(task, layout='g', name=str(task), scales=output_scales)
        solver.evaluate_handlers([output])
        # Check solution
        errors = []
        with h5py.File(f'{tempdir}/{tempdir}_s1.h5', mode='r') as file:
            for task in tasks:
                task_saved = file['tasks'][str(task)][-1]
                task = task.evaluate()
                if not isinstance(task, field.LockedField):
                    task.change_scales(output_scales)
                print(task['g'].shape, task_saved.shape, dealias, output_scales)
                errors.append(np.max(np.abs(task['g'] - task_saved)))
    assert np.allclose(errors, 0)