Here is a solution for counting ideals of a given norm in $\mathcal O_K$ with $K$ being a general number field. By relating the count of *proper* ideals with a given norm in a non-maximal order to counting ideals with a given norm in the maximal order containing it, I suspect the case of non-maximal orders could be deduced from the case of maximal orders, but I haven't thought that through.

For ${\rm Re}(s) > 1$, the Dirichlet series for $\zeta_K(s)$ can be written as
$\sum_{n \geq 1} a_n/n^s$ where $a_n$ is the number of ideals in $\mathcal O_K$ with norm $n$. We want to get an upper bound on $a_n$. Write the Euler product for $\zeta_K(s)$ as a product where all prime ideals dividing a given rational prime appear together:
$$
\zeta_K(s) = \prod_{p} \prod_{\mathfrak p \mid p} \frac{1}{1 - 1/{\rm N}(\mathfrak p)^s}.
$$
For each $\mathfrak p$ dividing $p$, ${\rm N}(\mathfrak p) \geq p$ and there are at most $[K:\mathbf Q]$ prime ideals in $\mathcal O_K$ dividing $p$, so
for *real* $s > 1$ we have $1/(1 - 1/{\rm N}(\mathfrak p)^s) \leq 1/(1 - 1/p^s)$ and therefore
$$
\zeta_K(s) \leq \prod_{p} \left(\frac{1}{1 - 1/p^s}\right)^m = \zeta(s)^m,
$$
where $m = [K:\mathbf Q]$. Thus the number $a_n = |\{\mathfrak a \subset \mathcal O_K : {\rm N}(\mathfrak a) = n\}|$ is at most the coefficient of $1/n^s$ in the Dirichlet series for $\zeta(s)^m$.

Multiplying Dirichlet series uses "Dirichlet convolution" on the coefficients:
$$
\sum_{n \geq 1} \frac{b_n}{n^s} \sum_{n \geq 1} \frac{c_n}{n^s} =
\sum_{n \geq 1} \left(\sum_{d \mid n} b_dc_{n/d}\right)\frac{1}{n^s}.
$$
In particular, taking $b_n = c_n = 1$, the sum $\sum_{d \mid n} b_dc_{n/d}$ is
the number of (positive) divisors of $n$. Writing that as $\tau(n)$,
we get $\zeta(s)^2 = \sum_{n \geq 1} \tau(n)/n^s$, so $a_n \leq \tau(n)$ when $[K:\mathbf Q] = 2$. That tells you what you want when $\mathcal O_d = \mathcal O_K$ for all quadratic $K$ (not just real quadratic $K$).

The Euler product of $\zeta(s)^m$ is $\prod_{p} 1/(1-1/p^s)^m$, so by the power series $1/(1-x)^m = \sum_{k \geq 0} \binom{k+m-1}{k}x^k$ for $|x| < 1$ with $x = 1/p^s$, the coefficient of $1/n^s$ in $\zeta(s)^m$ can be expressed in terms of the prime factorization $n = p_1^{k_1}\cdots p_r^{k_r}$ as
$$
\binom{k_1 + m-1}{k_1} \cdots \binom{k_r+m-1}{k_r}.
$$
Thus in general,
$$
[K:\mathbf Q] = m \text{ and } n = p_1^{k_1}\cdots p_r^{k_r} \Longrightarrow a_n \leq \binom{k_1 + m-1}{k_1} \cdots \binom{k_r+m-1}{k_r}.
$$
When $m = 2$, the upper bound is $(k_1+1)\cdots (k_r+1)$, which is $\tau(n)$.
For all $m$ we have $\binom{k+m-1}{k} \leq (k+1)^{m-1}$, so
$a_n \leq \tau(n)^{m-1}$ for all $n$. That is a small savings in the exponent compared to the bound $a_n \leq \tau(n)^m$ that I mentioned from Borevich and Shafarevich's *Number Theory* in a comment to the question, although I think the same savings could probably be obtained from B&S's own proof.

Edit: The use of Dirichlet series is not strictly necessary, as the number of ideals in $\mathcal O_K$ with norm $n = p_1^{k_1}\cdots p_r^{k_r}$ is the product of the number of ideals of norm $p_i^{k_i}$ for each $i$, and the count of ideals with prime-power norm $p^k$ can be bounded by the count that could occur if $p$ split completely (to maximize the number of possibilities). Nevertheless, the convenience of using Dirichlet series (really, an Euler product) is the same reason that it's often convenient to use generating functions to answer combinatorial questions.

Number Theoryis a proof that in all number fields $K$ (not just real quadratic fields), the number of ideals in $\mathcal O_K$ with norm at most $a$ -- which they write as $\psi(a)$ -- is $\leq \tau(a)^{[K:\mathbf Q]}$, where $\tau(a)$ is the number of divisors of $a$. For quadratic $K$ this upper bound on $\psi(a)$ is $\tau(a)^2$, not $\tau(a)$. In the proof, they make averyweak estimate near the end (cont). $\endgroup$