There are two approaches I can think of.

(a) Analytic. We need only information about $\zeta(s)$, not about other $L$-functions. Hence we can use the fact that RH has been verified up to a very large height. This should imply
a result of the form

$\sum_{n\; \text{squarefree}} f(n/x) = (1+\epsilon) 2 \zeta(2) + O^*(C_{\epsilon'} x^{1/4+\epsilon'})$,

where $O^*(K)$ means " a quantity of absolute value at most $K$", $\epsilon$ is a very tiny constant, $\epsilon'>0$ is arbitrarily small and $C_{\epsilon'}$ depends only on $\epsilon'$.

Difficulty: giving an explicit, and preferably small, value for $C_{\epsilon'}$ does not seem very easy. We cannot shift the line of integration all the way to $\Re(s)=1/2$ without having a bound on the residues of $1/\zeta(s)$; if we shift the line of integration only to $\Re(s)=1/2+\epsilon'$, we still need upper bounds for $1/|\zeta(s)|$ (i.e., lower bounds for $\zeta(s)$). There is some non-explicit work in this direction, but I do not know of anything explicit with reasonable constants. (I would be very glad to be surprised.)

(b) "Elementary". A fully elementary bound would presumably follow Cohen and Dress, MR0952866, and give bounds of about the same quality (error term = $O(x^{1/2})$, where
the implied constant is not large but also not very small). A more mixed approach would proceed as in the paper by Cohen, Dress and El Marraki cited by Barry Cipra above. Half of the work (sieving out factors m^2 with $m\leq c\cdot \sqrt{x}$) would indeed be improved by the smoothing. The other half would rely on estimates on $M(N) = \sum_{n\leq N} \mu(n)$. Here estimates of the form $|M(N)|\leq c N$, $c$ a small constant (for $N$ larger than a constant) are known, and lead to estimates of the form

$\sum_{n\; \text{squarefree}} f(n/x) = 2 \zeta(2) + O^*(c x^{1/2})$,

with $c$ rather small (but not "very tiny"; we are speaking about $10^{-3}$ rather than $10^{-20}$ , say) provided that $x$ is larger than a constant. Getting useful error-term estimates better than $O^*(c x^{1/2})$ is hard, however, in that it involves estimating $M(x)$: this seems hard for the same reason given above, namely, estimates for the residues of $1/\zeta(s)$ or upper bounds for $1/|\zeta(s)|$ become necessary.

Question: is there a third way? Or is there a way to give good upper bounds for $1/|\zeta(s)|$ for $|\Im(s)|\leq H_0$, assuming that the Riemann hypothesis holds for $|\Im(s)|\leq H_0$ (or a little further)?