Since $x^2+y^2+z^2$ is conserved, we can switch to the stereographic projection, i. e. substitute $x=\frac{2X}{X^2+Y^2+1}R$, $y=\frac{2Y}{X^2+Y^2+1}R$, $z=\frac{X^2+Y^2-1}{X^2+Y^2+1}R$, with constant $R$. This gives
\begin{align*}
\frac{dX}{dt}&=\frac12\cos(\omega t)(X^2-Y^2+1)-\sin(\omega t)XY,\\
\\
\frac{dY}{dt}&=\frac12\sin(\omega t)(X^2-Y^2-1)+\cos(\omega t)XY.
\end{align*}
Denoting $Z=X+iY$, we reduce to a single equation
$$
\frac{dZ}{dt}=\frac12\left(e^{i\omega t}Z^2+e^{-i\omega t}\right),
$$
solved as
$$
Z(t)=e^{-i\omega t}\left(\sqrt{1+\omega^2}\tan\left(\frac12\sqrt{1+\omega^2}t+C_0+iC_1\right)-i\omega\right)
$$
(with arbitrary real constants $C_0$, $C_1$).

The latter can be found by e. g. further switching to $F(t)=e^{i\omega t}Z(t)$ which satisfies
$$
\frac{dF}{dt}=\frac12\left(1+2i\omega F+F^2\right),
$$
giving
$$
t=C+\int\frac{2dF}{1+2i\omega F+F^2}=C_0+iC_1+\frac2{\sqrt{1+\omega^2}}\arctan\left(\frac{F+i\omega}{\sqrt{1+\omega^2}}\right)
$$

Returning to $(x,y,z)$ it is easy to see that this describes circular motion on a sphere around some axis like this

with the sphere itself at the same time spinning around the $z$ axis (because of the $e^{i\omega t}$ multiplier in $Z=e^{-i\omega t}F$), which means that there must be an easier method to solve the original system. In fact I suspect these are a case of spinning top equations...