Call a function $F$ nice, if $$\DeclareMathOperator{\Dm}{d\!} \begin{align} &\zeta(2)\sum_{\frac{a}{b}\in\mathbb{Q}_n} \frac{F(\log \frac{a}{b})}{\sqrt{abn}}\xrightarrow[n\to \infty]{}\int F(x)\Dm x \qquad \label{1}\tag{1}\\ \textrm{where } \quad &\mathbb{Q}_n = \{\tfrac{a}{b}\in \mathbb{Q}^+,\textrm{ gcd}(a,b)=1, ab\leq n\} \quad . \end{align} $$ Obviously, nice functions form an $\mathbb{R}$-linear space. Also, assume that for some function $F$ for every $\varepsilon>0$ there exist nice functions $F_1,F_2$ such that $F_1\leqslant F\leqslant F_2$ and $\int (F_2-F_1)<\varepsilon$. Then $F$ is nice itself. Indeed, $$ \begin{split} \limsup_{n\to \infty} \zeta(2)\sum_{\frac{a}{b}\in\mathbb{Q}_n} \frac{F(\log \frac{a}{b})}{\sqrt{abn}} & \leqslant \limsup_{n\to \infty} \zeta(2)\sum_{\frac{a}{b}\in\mathbb{Q}_n} \frac{F_2(\log \frac{a}{b})}{\sqrt{abn}}\\ &=\int F_2\Dm x \leqslant \varepsilon+\int F_1 \Dm x\\ &\leqslant \varepsilon+\int F\Dm x, \end{split} $$ and, since $\varepsilon>0$ is arbitrary, the limsup of LHS of \eqref{1} does not exceed $\int F$. Analogously, liminf is not less than $\int F$. Therefore the liminf and limsup are both equal to $\int F$, thus $F$ is nice.
If we prove that for every real $p<q$ the characteristic function of the segment $[p,q]$ is nice, by above observation automatically we get that every properly Riemann integrable function is nice.
So, fix $p<q$ and consider $F=\chi_{[p,q]}$. Then \eqref{1} reads as $$\zeta(2)\sum_{e^p\leqslant a/b\leqslant e^q,\frac{a}{b}\in\mathbb{Q}_n} \frac{1}{\sqrt{ab}}\sim (q-p)\sqrt{n}.\label{2}\tag{2}$$ Fix positive but small constant $c$. Partition the sum in \eqref{2} by two parts: the sum $S_c$ corresponds to those points $(a,b)$ for which $\min(a,b)\geqslant c\sqrt{n}$, the sum $R_c$ ($R$ for remainder) to those points $(a,b)$ for which $\min(a,b)< c\sqrt{n}$. Both sums are over regions $\sqrt{n}\cdot \Omega_s,\sqrt{n}\cdot \Omega_r$ respectively which are homothetic images of fixed regions $\Omega_s$, $\Omega_r$.
By the uniform distribution of coprime points, we get $$S_c\sim \sqrt{n}\int_{\Omega_c}\frac{\Dm x\Dm y}{\sqrt{xy}}$$ ($\Omega_c$ is away of all singularities of the integrand, so there are no problems). For $R_c$, we may bound the sum over coprime points by the sum of all integer points with positive coordinates, which may be bounded from above by the corresponding integral. Since $\int_{x,y>0}\frac{\Dm x\Dm y}{\sqrt{xy}}$ converges, the bound we get for $R_c$ is something small (when $c$ is small) times $\sqrt{n}$, and the integral over $\Omega_c$ is close (again if $c$ is small) to the integral $$\int_{x,y>0,e^p<x/y<e^q,xy<1}\frac{\Dm x\Dm y}{\sqrt{xy}}=q-p.$$
