$$
\newcommand{\qcl}{q_{\rm cl}}
\newcommand{\ket}[1]{| #1 \rangle}
\newcommand{\bra}[1]{\langle #1 |}
$$
This is a source of very sloppy work which appears in many textbooks. You are completely correct that it makes no sense to divide by this diverge factor ad hoc. The reason they are doing this is because they weren't careful enough with the measure of the path integral which essentially provides the perfect divergent factor you need.
When calculating this quantity, you decomposed variations around the classical path into Fourier modes. This change of variables in the path integral comes with an associated (divergent) Jacobian factor $J_N$. Before getting into the Harmonic oscillator, let's start with the free Hamiltonian when $\omega = 0$. Because this Jacobian factor doesn't depend on the Hamiltonian, we can use the well known expression for the heat kernel of the free Hamiltonian to solve for it. After we extract this factor we will then move to the Harmonic oscillator and use it to compute the transition amplitude exactly.
(While I will present this work in terms of real time and a transition amplitude instead of Euclidean time and the traced partition function, it is a simple matter to go between the two.)
We will follow Nakahara's "Geometery, Topology, and Physics second edition" section 1.4.
Let us begin with the well known expression for the Heat kernel of the free particle Hamitlonian $H_{\rm free}$.
\begin{align}\tag{1}\label{eq1}
\bra{q_f}e^{- i H_{\rm free} T / \hbar} \ket{q_i} &= \left( \frac{m/\hbar}{2 \pi i T} \right)^{1/2} \exp( \frac{i}{\hbar} \frac{m}{2 T} (q_f - q_i)^2 ) \\
&= \left( \frac{m/\hbar}{2 \pi i T} \right)^{1/2} e^{i S_{\rm free}[q_{\rm cl} ]/\hbar}.
\end{align}
Recall also that the measure of our path integral is defined as
\begin{equation}
\int \mathcal{D} q(t) = \lim_{N \to \infty} \left( \frac{m/ \hbar }{2 \pi i \epsilon}\right)^{(N+1)/2} \int dq_1 \ldots dq_N
\end{equation}
where
\begin{equation}
\epsilon = \frac{T}{N+1}.
\end{equation}
For the free particle, the action is
\begin{equation}
S_{\rm free}[q] = \frac{m}{2} \int_0^T dt \; \dot q^2.
\end{equation}
If we break up our paths as
\begin{equation}
q = \qcl + y
\end{equation}
then, for this simple quadratic action,
\begin{align}
S_{\rm free}[q] = S_{\rm free}[\qcl] + S_{\rm free}[y]
\end{align}
and
\begin{align}
\bra{q_f}e^{- i H_{\rm free} T / \hbar} \ket{q_i} &= \int_{q(0) = q_i}^{q(T) = q_f} \mathcal{D} q(t) e^{i S_{\rm free}[q] / \hbar} \\
&= e^{i S_{\rm free}[\qcl] / \hbar} \int_{y(0) = 0}^{y(T) = 0} \mathcal{D}y(t) e^{i S_{\rm free}[y] / \hbar}.
\end{align}
Combining this with our expression \eqref{eq1}, we get
\begin{equation}\label{heatkernel2}\tag{2}
\left( \frac{m/\hbar}{2 \pi i T} \right)^{1/2} = \int_{y(0) = 0}^{y(T) = 0} \mathcal{D} y(t) e^{i S_{\rm free}[y]/\hbar}
\end{equation}
which is our starting point.
We now consider $y(t)$ in terms of its Fourier coefficients.
\begin{align}
y(t) = \sum_{n = 1}^N a_n \sin( \frac{n \pi t}{T}).
\end{align}
We have only allowed $N$ coefficients for the following reason: our integral is an integral over the $N$ variables $y_k$ for $k = 1 \ldots N$, and When we change the variables of integration to $a_n$, we must preserve the number of variables we are integrating over. The Jacobian of this transformation, which our main object of interest, is
\begin{equation}
J_N \equiv \det( \frac{\partial y_k}{\partial q_n} ) = \det( \sin(\frac{n \pi t_k}{T} ) ).
\end{equation}
The action is now
\begin{align}
S_{\rm free}[y] = \frac{m}{2}\sum_{n=1}^N\frac{T}{2} \left( \frac{a_n n \pi}{T} \right)^2
\end{align}
and Eq. \ref{heatkernel2} becomes
\begin{align}
\left( \frac{m/\hbar}{2 \pi i T} \right)^{1/2} &= J_N \left( \frac{m /\hbar}{2 \pi i \epsilon} \right)^{(N+1)/2} \int d a_1 \ldots d a_N \exp( \frac{i m}{\hbar} \sum_{n=1}^N \frac{a_n^2 \pi^2 n^2}{4 T} ) \\
&= J_N \left( \frac{m /\hbar}{2 \pi i \epsilon} \right)^{(N+1)/2} \prod_{n=1}^N \left( - \pi \frac{ 4 T \hbar}{i m \pi^2 n^2} \right)^{1/2} \\
&= J_N \left( \frac{m /\hbar}{2 \pi i \epsilon} \right)^{(N+1)/2} \left( \frac{i 4 \pi T \hbar}{m \pi^2} \right)^{N/2} \frac{1}{N!}.
\end{align}
After some clean up, this becomes
\begin{equation}
J_N =2^{-N/2} \pi^N (N+1)^{-(N+1)/2} N!
\end{equation}
which happily is dimensionless.
Now let's see how we can use $J_N$ to get a finite result for the transition amplitude of the Harmonic oscillator, which has the action of
\begin{equation}
S[q] = \frac{m}{2} \int_0^T \left( \dot q^2 - \omega^2 q^2 \right) dt.
\end{equation}
The simple and quadratic nature of this action still means that
\begin{equation}
S[\qcl + y] = S[\qcl] + S[y].
\end{equation}
The action of our Fourier expanded paths is
\begin{align}
S[y] &= \frac{m}{2} \sum_{n = 1}^N \frac{T}{2} a_n^2 \left( \left(\frac{n \pi}{T}\right)^2 - \omega^2 \right) \\
&= \frac{T}{4} \sum_{n = 1}^N a_n^2 \lambda_n
\end{align}
where we have defined
\begin{equation}\label{eigenvaluefirst}
\lambda_n = m \left( \left(\frac{n \pi}{T}\right)^2 - \omega^2 \right).
\end{equation}
Therefore,
\begin{align}
\bra{q_f} e^{- i H T/\hbar} \ket{q_i} =& e^{\tfrac{i}{\hbar} S[\qcl] } \int \mathcal{D} y(t) e^{\tfrac{i}{\hbar} S[y] } \\
=& e^{\tfrac{i}{\hbar} S[\qcl] } J_N \left( \frac{m/ \hbar }{2 \pi i \epsilon}\right)^{(N+1)/2} \int d a_1 \ldots d a_N \exp( \frac{i T}{4 \hbar} \sum_{n = 1}^N a_n^2 \lambda_n ) \\
=& e^{\tfrac{i}{\hbar} S[\qcl] } J_N \left( \frac{m /\hbar}{2 \pi i \epsilon} \right)^{(N+1)/2} \left(\frac{i 4 \pi \hbar}{T} \right)^{N/2} \left( \prod_{n = 1}^N \lambda_n \right)^{-1/2} \\
=& e^{\tfrac{i}{\hbar} S[\qcl] } \left( \frac{m/\hbar}{2 \pi i T} \right)^{1/2} \left( \frac{m \pi^2 }{T^2} \right)^{N/2} N! \left( \prod_{n = 1}^N \lambda_n \right)^{-1/2}.
\end{align}
Note the terms divergent in $N$. If we were using the free action, with $\omega = 0$, then
\begin{equation}
\lambda^{\omega = 0}_n = m \left( \frac{n \pi}{T} \right)^2, \hspace{1 cm} \left(\prod_{n=1}^N \lambda^{\omega =0}_n \right)^{1/2} = \frac{m^{N/2} \pi^N}{T^N} N!
\end{equation}
and thus
\begin{equation}\label{ratioeigenvalues}
\bra{q_f} e^{- i H T/\hbar} \ket{q_i} = e^{\tfrac{i}{\hbar} S[\qcl] } \left( \frac{m/\hbar}{2 \pi i T} \right)^{1/2} \frac{\left( \prod_{n = 1}^N \lambda_n \right)^{-1/2}}{\left( \prod_{n = 1}^N \lambda_n^{\omega = 0} \right)^{-1/2}}.
\end{equation}
We now compute
\begin{equation}
\prod_{n=1}^N\frac{\lambda_n}{\lambda_n^{\omega=0}} = \prod_{n=1}^N \left( 1 - \left( \frac{\omega T}{n \pi} \right)^2 \right) = \frac{\sin(\omega T) }{\omega T}
\end{equation}
making
\begin{equation}
\bra{q_f} e^{- i H T/\hbar} \ket{q_i} = \left( \frac{m \omega}{\hbar \pi} \frac{1}{2 i \sin ( \omega T) }\right)^{1/2} e^{\tfrac{i}{\hbar} S[\qcl] }.
\end{equation}
