Solution (source code)

= Solution

\b[The statement as printed is false without an integrability condition.] On any nonempty open set, the constant functions $f_n\equiv n$ are harmonic and therefore subharmonic, but their supremum is $+\infty$, which is not locally integrable.

A valid version assumes $g=\sup_n f_n\in L^1_{\rm loc}(\Omega)$. First, the maximum of finitely many subharmonic functions is subharmonic. For a distributional proof, mollify two functions on smaller interior domains and replace their maximum by
$$
M_\delta(s,t)=\frac12\left(s+t+\sqrt{(s-t)^2+\delta^2}\right).
$$
This is smooth, convex, and nondecreasing in each argument. The <chain rule> shows that its Laplacian applied to two smooth subharmonic functions is nonnegative: the first-derivative terms multiply their nonnegative Laplacians, and the Hessian term is nonnegative by convexity. <Mollifier> convergence in $L^1_{\rm loc}$ and the uniform bound $|M_\delta-\max|\leq\delta/2$ allow both regularizations to be removed. Continuity of <distributional derivatives> under this convergence proves the finite-maximum assertion.

Now $g_N=\max_{n\leq N}f_n$ increases to $g$ and satisfies $|g_N|\leq|f_1|+|g|$ almost everywhere. The <dominated convergence theorem> gives convergence in $L^1_{\rm loc}$. For every nonnegative test function,
$$
\int g\,\Delta\varphi=\lim_{N\to\infty}\int g_N\,\Delta\varphi\geq0.
$$
Thus $g$ is subharmonic in the distributional sense. This is the <locally integrable supremum theorem for subharmonic functions>. For the pointwise upper semicontinuous convention, take its canonical representative; the raw supremum need not already have that regularity.