You have not told us what your working definition of Riemann integrability is.
If integrability is defined in terms of lower and upper sums, and $f$ is unbounded above, then every single upper sum is undefined (or $=\infty$, if you prefer). Basta; no further epsilontics needed.
If integrability is defined in terms of quantities like $\|\Delta f\|_{I_j}:=\sup_{x, \ y\in I_j}|f(x)-f(y)|$ then again for every partition at least one of these quantities is $=\infty$, hence you can never make $\sum_j\|\Delta f\|_{I_j}\>|I_j|<\epsilon$. But this is required before you can even think of computing "general Riemann sums" $\sum_{j=1}^N f(\tau_j)\>|I_j|$ as approximations to the intended integral.
Consider a partition $a = x_0 < x_1 < \ldots < x_n = b$ and for a subinterval $I_j = [x_{j-1},x_j]$ define
$$D_{j}(g) = \sup_{x \in I_j}g(x) - \inf_{x \in I_j}g(x) = \sup_{x,y \in I_j}|g(x) - g(y)| \\ D_{j}(f \circ g ) = \sup_{x \in I_j}f(g(x)) - \inf_{x \in I_j}f(g(x)) = \sup_{x,y \in I_j}|f(g(x)) - f(g(y))|$$
Since $f$ is continuous it is uniformly continuous and bounded (by extension if necessary) on a closed interval $[c,d]$ such that $g([a,b]) \subset [c,d].$
Hence, $|f(x)| \leqslant M$ for $x \in [c,d]$ and for every $\epsilon >0$ there exists $\delta > 0$ such that if $|x_1 - x_2| < \delta$ then $|f(x_1) - f(x_2)| < \epsilon/(2(b-a))$ .
Since $g$ is integrable, if the partition norm $\|P\|$ is sufficiently small we have
$$U(P,g) - L(P,g) = \sum_{j=1}^n D_j(g) (x_j - x_{j-1}) < \frac{ \delta \epsilon}{4M}.$$
We can split the upper-lower sum difference $U(P,f \circ g) - L(P, f \circ g)$ into two sums as given by
$$\tag{1}U(P,f \circ g) - L(P, f \circ g) = \sum_{D_j(g) \geqslant \delta} D_j(f \circ g)(x_j - x_{j-1}) + \sum_{D_j(g) < \delta} D_j(f \circ g)(x_j - x_{j-1})$$
In the second sum on the RHS of (1) we have $D_j(f \circ g) < \epsilon/(2(b-a)$ since by uniform continuity $D_j(g) < \delta \implies |g(x) - g(y)| < \delta \implies |f(g(x)) - f(g(y))| < \epsilon/(2(b-a)$ for all $x,y \in I_j$.
Thus,
$$\tag{2}\sum_{D_j(g) < \delta} D_j(f \circ g)(x_j - x_{j-1}) < \frac{\epsilon}{2}$$
Considering the first sum on the RHS of (1), first note that
$$\sum_{D_j(g) \geqslant \delta} (x_j - x_{j-1}) \\ < \delta^{-1}\sum_{D_j(g) \geqslant \delta} D_j(g)(x_j - x_{j-1}) < \delta^{-1} [U(P,g) - L(P,g)] < \delta^{-1} \frac{\delta \epsilon}{4M} = \frac{\epsilon}{4M} .$$
Hence,
$$\tag{3}\sum_{D_j(g) \geqslant \delta} D_j(f \circ g)(x_j - x_{j-1}) < \sum_{D_j(g) \geqslant \delta} 2M(x_j - x_{j-1}) < \frac{\epsilon}{2}.$$
From (1), (2) and (3) we obtain
$$U(P,f \circ g) - L(P, f \circ g) < \epsilon,$$
and conclude that $f \circ g$ is integrable.
Best Answer
The definition of upper and lower integrals is $\mathcal{U}(f) = \inf_{P'} U(f,P')$ and $\mathcal{L}(f) = \sup_{P'} L(f,P')$ where the infimum and supremum are taken over all partitions $P'$ of $[a,b]$.
In other words, the upper integral $\mathcal{U}(f)$ is the greatest lower bound of all upper sums and for any particular partition $P$ we have $U(f,P) \geqslant \mathcal{U}(f)$. Similarly we have $L(f,P) \leqslant \mathcal{L}(f)$.
Thus, $$\tag{*}U(f,P) - L(f,P) \geqslant \mathcal{U}(f) - \mathcal{L}(f)$$
In proving the reverse implication you first assume that $f$ is not Riemann integrable. This implies that $\mathcal{U}(f) > \mathcal{L}(f)$ and so $\alpha := \mathcal{U}(f) - \mathcal{L}(f) > 0$.
Take $\epsilon = \alpha/2$. By (*) it follows that for any partition $P$ we have
$$U(f,P) - L(f,P) \geqslant \mathcal{U}(f) - \mathcal{L}(f) = \alpha > \frac{\alpha}{2} = \epsilon$$
This contradicts the hypothesis that for any $\epsilon$ there is a partition $P$ such that $U(f,P) - L(f,P) < \epsilon$. Therefore, $f$ must be Riemann integrable.