In deriving Planck's blackbody formula, the number density of normal modes (per unit frequency$^\dagger$) is found, given by $$ N(\omega)=\frac{V}{\pi^2c^3}\omega^2, $$ where $V$ is the volume of the blackbody. Then the average energy of a mode of frequency $\omega$ is calculated using Planck's quantization that energy of a normal mode with frequency $\omega$ can only be $n\hbar\omega$ for some whole number $n$.
But hold on a minute!
The energy of a linear quantum oscillator with frequency $\omega$ can be $(n+1/2)\hbar\omega$, and since in deriving Planck's formula, it is the analogy of normal modes of the electromagnetic waves inside the blackbody cavity with the normal modes of the linear oscillator is made, I'd expect that we take into account this "zero-point energy", $\epsilon_\omega = \hbar\omega/2$ for each normal mode. Note that this zero-point energy, $\epsilon_\omega$ depends on $\omega$.
But if I do take this into account then I get bizarre results!
The average energy of normal mode of frequency $\omega$ is then given by $$ \langle E_\omega \rangle = \frac{\hbar\omega}{e^{\hbar\omega/k_BT}-1} + \epsilon_\omega, $$ and this leads to the following energy density.
\begin{align} \rho(T, \omega) &:= \frac{1}{V} N(\omega)\langle E_\omega \rangle\\ &\;= \frac{\hbar}{\pi^2c^3}\frac{\omega^3}{e^{\hbar\omega/k_BT} - 1} + \underbrace{\frac{1}{\pi^2c^3}\omega^2\epsilon_\omega}_{\text{additional term}} \end{align}
Now, this result is catastrophic! Unless $\epsilon_\omega = 0$ for all $\omega$'s (in which case this will coincide with the "correct" Planck's formula), $\rho$ diverges as $\omega\to\infty$.
Questions: So what's the way out? Is the often-presented analogy with quantum oscillator plainly wrong? For electromagnetic radiation, is the zero-point energy exactly zero for all $\omega$'s?
$^\dagger$ By frequency, I mean angular frequency.