Typically the solution for the inhomogeneous problem is obtained via the homogeneous problem via Duhamel's principle. In short, the solution to your equation (with zero initial conditions) is
$$
u(t, x) = \int_0^t v(t, x; s)\, ds
$$
where $v(\cdot, \cdot; s) : (0, \infty) \times \mathbb{R}^3 \to \mathbb{R}$ is the solution to the equation
$$
\begin{cases}
(\partial_t^2 - \Delta)v(t, x; s) = 0 \\
v(t, x; s) = 0, v_t(t, x;s) = f(x, s). \tag{$1s$}
\end{cases}
$$
Note that in 3D this $v$ is given by
$$
v(t, x; s) = \frac{1}{|\partial B(x, t)|}\int_{\partial B(x,t)} th(y)\, dS(y).
$$
Therefore
$$
u(t, x) = \int_0^t \frac{1}{|\partial B(x, t)|}\int_{\partial B(x,t)} th(y)\, dS(y)\, dt.
$$
See Evans for more details. One can formulate this question with wave propagators as well; in that case the solution is
$$
u(t, x) = \cos(t\sqrt{-\Delta})u_0(0,x) + \frac{\sin(t\sqrt{-\Delta})}{\sqrt{-\Delta}}u_1(0, x) - \int_0^t \frac{\sin((t - s)\sqrt{-\Delta})}{\sqrt{-\Delta}}F(s)\, ds.
$$
(See Exercise 2.22 of Tao's book on nonlinear dispersive equations.) Here $u_0(x) = u(0, x)$ and $u_1(x) = \partial_t u(0, x)$ are general initial conditions.
Summary of comments:
The integral diverges because of the non-integrable poles on the integration contour. Since these poles are of first order, the Cauchy principle value (PV) of the integral exists, and corresponds to avoiding the poles with semi-circular arcs. The PV is unchanged whether we make an arc above or below the pole. This is a type of regularization. Other regularizations may give different results. Even though we will be integrating above/ below the poles, this is different to the popular $\pm i \epsilon$ prescription whereby one considers a related integrand with poles shifted above or below the axis$^\dagger$ (the $\pm i\epsilon $ regularization leads to results that do depend on the contour taken).
The principle value is
$$\tag{1}
\text{PV}\int_{-\infty}^\infty d\omega \frac{ e^{-i\omega t}}{\omega^2-k^2} = -\frac{\pi}{k}\sin(k|t|)
$$
Consider the case $t>0$, so we close the contour in the lower half plane. I've also set $c=1$ for my convenience. Suppose we make arcs above the poles:
Then we have by residue theorem
$$\tag{2}
I+C_-+C_++C_R=-2\pi i (R_- +R_+)
$$
Where $I$ is the integral in eq (1), $C_{\pm}$ are the contributions due to integrating over the poles at $\pm k$, $C_R$ is the (vanishing) contribution due to the large arc, $R_\pm$ are the residues at $\pm k$, and the minus sign on the RHS is because the contour is CW.
Near the pole $\omega=k$ we can write
$$
\frac{ e^{-i\omega t}}{\omega^2-k^2} = \frac{e^{-ikt}}{2k(\omega-k)} + \text{finite}
$$
To integrate around the pole, parameterize: $\omega=k+\epsilon e^{i\phi}$, then $d\omega=i \epsilon e^{i\phi}d\phi$. Integrating around a CCW semicircle yields
$$\tag{3}
\int\limits_0^\pi \frac{i e^{-ikt}}{2k} d\phi + \mathcal{O}(\epsilon) = i\pi \left(\frac{ e^{-ikt}}{2k}\right) + \mathcal{O}(\epsilon)
$$
When we take the radius of the arc, $\epsilon$ to zero, all the other terms vanish. The result is precisely $i\pi$ times the residue at the pole. This will generally be the case for simple first order poles$^\ddagger$. Substituting this and the similar calculation for the other pole into eq (2) we have
$$
I-i\pi R_- - i \pi R_+ = -2 \pi i (R_-+R_+)
$$
The new minus signs on the LHS are because of the orientation of $C_\pm$. We are left with
$$
I=-i\pi(R_-+R_+)=-i\pi \left(-\frac{e^{ikt}}{2k}+\frac{e^{-ikt}}{2k} \right)=-\frac{\pi}{k}\sin(kt)
$$
Suppose we instead made the arcs beneath the poles (still with $t>0$). The analogue of eq (2) is
$$
I+C'_-+C'_+=0
$$
Where $C'_\pm$ are the contributions due to CCW arcs beneath the poles (already calculated in eq (3), they differ only by a sign from the $C_\pm$). The RHS is zero as this contour encloses no poles. Substituting in
$$
I+i\pi R_-+i\pi R_+=0 \\
I=-i \pi (R_- + R_+)
$$
Precisely as before. You can check that other combinations of above/ below the poles yield the same result. When $t>0$, the solution will be (slightly) different, but using the same method you can show that the answer is independent of whether your chosen contour goes above or below the poles.
$\dagger$ A fact that even some well respected textbooks forget.
$\ddagger$ You can see why this only works for first order poles, if there were any terms more singular than $1/\omega$ in eq (3), there would remain contributions proportional to $1/\epsilon$ which prevent the limit of zero radius being taken safely.
Best Answer
You haven't done anything wrong, rather your question is actually ill-posed. Notice that your problem does not have any boundary conditions so it can't have an unambiguous answer (what happens if you add a constant to $\phi$). If your source term cut off at some point in the past, say $\kappa(z, t) = e^{- i \omega t + i \omega z} \Theta_L(z) \Theta(t - t_0)$ then your integral would look like, $$ \int_{t_0}^{t - |z - z'|} \mathrm{d}t' e^{- i \omega t'} = \frac{i}{\omega} \left( e^{- i \omega (t - |z - z'|)} - e^{- i \omega t_0} \right) $$
which is perfectly well behaved. This corresponds to implicitly imposing boundary conditions $\phi(z, t) = 0$ in the past for $t \le t_0$ (I encourage you to think about how boundary contitions are incorporated into the general solution throught the Green's function). However, in your case you are taking $t_0 \to - \infty$ but we can't impose this sort of boundary condition at $- \infty$ since $\phi(z, -\infty) = 0$ isn't meaningful! This correspond to the fact that, for the solutions $\phi_{t_0}$ computed for a cutoff set at $t_0$, the limit $\lim\limits_{t_0 \to - \infty} \phi_{t_0}(z, t)$ does not exist. This is exactly the same ill-defined limit that you noticed when you couldn't compute the limit in the improper integral, $$ \int_{-\infty}^{t - |z - z'|} \mathrm{d}t' e^{- i \omega t'} = \lim_{t_0 \to - \infty} \int_{t_0}^{t - |z - z'|} \mathrm{d}t' e^{- i \omega t'} = \lim_{t_0 \to - \infty} \frac{i}{\omega} \left( e^{- i \omega (t - |z - z'|)} - e^{- i \omega t_0} \right) $$ Now what would make your problem well-posed while retaininig the same source term $\kappa(z, t)$. Well, suppose at some $t_0$ we know the value of $\phi(z, t_0)$ and $\partial_t \phi(z, t) |_{t_0}$. Then we need to modify our Green's function $G(z,z',t,t')$ to take this information into account. Where $G$ satisfies, $$ (\partial_t^2 - \partial_x^2) G(z, z', t, t') = \delta(t - t') \delta(x - x') $$ and the Green's function must be a function of $t$ and $t'$ (not of the form $G(z - z', t - t')$ as you had before since we require that $G(z, z', t_0, t')$ satisfy the boundary contition for all $t'$). Suppose we know $\phi(z, t_0) = 0$ and $\partial_t \phi(z, t)|_{t_0} = 0$. You will find something piecewise like, $$ G(z, z', t, t') = \begin{cases} \Theta(t - t')\Theta(t - t' - |z - z'|) & t' > t_0 \\ \Theta(t' - t)\Theta(t' - t - |z - z'|) & t' < t_0 \end{cases} $$ Notice the Green's function is advanced before $t_0$ and retarded afterwards conforming to the sort of causality we expect for the propagation of information about the solution at $t_0$ to information about the solution at all $t$. Now our solution takes the form, $$ \phi(z, t) = \int_{-\infty}^{\infty} \mathrm{d}{z'} \int_{-\infty}^{\infty} \mathrm{d}{t'} G(z, z', t, t') \kappa(z', t') $$ which becomes $$ \phi(z, t) = \int_{0}^{L} \mathrm{d}{z'} \begin{cases} \int_{t_0}^{t - |z - z'|} \mathrm{d}{t'} e^{- i \omega t' + i \omega z'} & t > t_0 \\ \int_{t + |z - z'|}^{t_0} \mathrm{d}{t'} e^{- i \omega t' + i \omega z'} & t < t_0 \end{cases} \quad = \frac{i}{\omega} \int_{0}^{L} \mathrm{d}{z'} e^{i \omega z'} \begin{cases} \Theta(t - t_0 - |z - z'|)\left( e^{- i \omega (t - |z - z'|)} - e^{- i \omega t_0} \right) & t > t_0 \\ \Theta(t_0 - t - |z - z'|)\left(-e^{- i \omega (t + |z - z'|)} + e^{- i \omega t_0} \right) & t < t_0 \end{cases} $$ and thus, $$ \phi(z, t) = \frac{i}{\omega} \mathrm{sign}(t - t_0) \left( \int_{0}^{L} \mathrm{d}{z'} \Theta(|t - t_0| - |z - z'|) \left( e^{i \omega z'} e^{\mathrm{sign}(t - t_0) i \omega | z - z'|} e^{- i \omega t} - e^{- i \omega t_0} \right) \right) $$ If we replace our source with a delta function at the origin $\kappa(z, t) = e^{- i \omega t} \delta(z)$ it is easier to see what is going on. In this case we get, $$ \phi(z, t) = \frac{i}{\omega} \mathrm{sign}(t - t_0) \Theta(|t - t_0| - |z|) \left( e^{\mathrm{sign}(t - t_0) i \omega | z |} e^{- i \omega t} - e^{- i \omega t_0} \right) $$ which after $t_0$ is a solution with outgoing waves from the origin and before $t_0$ is a solution with incoming waves exactly absorbed at the origin such that at $t = t_0$ all the waves cancel.