Theorem: Let $\Omega \subset \mathbb{R}^N$, $u \in C(\Omega)$ be such that $$\frac{1}{|B(x_0,R)|}\int_{B(x_0,R)}u(y)\ dy = u(x_0) = \frac{1}{|\partial B(x_0,R)|}\int_{\partial B(x_0,R)}u\ dS$$ for every ball $\overline{B(x_0,R)} \subset \Omega$. Then $u \in C^{\infty}(\Omega)$ and it is harmonic
Proof: Consider the standard mollifier: $$\rho(x) := \begin{cases}Ce^{-\frac{1}{1 - \|x\|^2}} & \text{if $\|x\|$ < 1} \\0 & \text{otherwise.}
\end{cases}$$ Here $C$ is a constant such that $\|\rho\|_{L^1} = 1.$ Let $\epsilon > 0$ and consider $$\rho_{\epsilon}(x) = \epsilon^{-N}\rho(x\epsilon^{-N}).$$
Set $\Omega_{\epsilon} = \{x \in \Omega : \text{dist}(x,\partial \Omega) > \epsilon\}$ and define for $x \in \Omega_{\epsilon}$$$u_{\epsilon}(x) = \rho_{\epsilon} * u(x) = \int_{\Omega}\rho_{\epsilon}(x - y)u(y)\ dy.$$ The following is a well know theorem in analysis, if it is new to you you can look for a proof Analysis by Lieb and Loss or anywhere else.
**Theorem:***If $u \in C(\Omega)$, then $u_{\epsilon} \to u$ uniformly on compact subsets of $\Omega$, $u_{\epsilon} \in C^{\infty}(\Omega_{\epsilon})$ and for any multindex $\alpha$ we have $$\frac{\partial^{\alpha}u_{\epsilon}}{\partial x^{\alpha}}(x) = \int_{\Omega}\frac{\partial^{\alpha}\rho_{\epsilon}}{\partial x^{\alpha}}(x - y)u(y)\ dy.$$*
Finally we can proceed with the proof!
Fix $x_0 \in \Omega_{\epsilon}$.
$$u_{\epsilon}(x_0) = \int_{B(x_0,\epsilon)}\rho_{\epsilon}(x - y)u(y)\ dy = \int_{B(0,\epsilon)}\rho_{\epsilon}(z)u(x_0 - z)\ dz = $$ $$ = \int_0^{\epsilon}r^{N - 1}\int_{\partial B(0,1)}\rho_{\epsilon}(rw)u(x_0 - rw)\ dS(w)dr = $$ $$ \int_0^{\epsilon}r^{N - 1}\rho(r)\int_{\partial B(0,1)}u(x_0 - rw)\ dS(w)dr = \int_0^{\epsilon}r^{N-1}\rho_{\epsilon}(r)\frac{\alpha_N N}{|\partial B(x_0,r)|}\int_{\partial B(x_0,r)}u(y)\ dS(y)dr $$ $$ = u(x_0)\|\rho\|_{L^1} = u(x_0).$$
This proves that $u = u_{\epsilon}$ and hence $u \in C^{\infty}(\Omega_{\epsilon})$, for every $\epsilon$. We are left to prove that $u$ is harmonic. To this end consider $$f(r) = \frac{1}{|\partial B(x_0,r)|}\int_{\partial B(x_0,r)}u\ dS.$$
By assumption, $f$ is constant, hence $f' \equiv 0$. By the divergence theorem it is easy to show (and I am sure that you have already seen this result since you have proved that if $u$ is harmonic then it satisfies the mean value property) $$0 = f'(r) = \text{constant}\int_{B(x_0,r)}\Delta u(y)\ dy \to \Delta u(x_0),\ \text{if we let}\ r \to 0^+.$$ Thus $u$ is harmonic. QED
Best Answer
Here's a hint to get you started:
Define $A(r) = \int_{\partial B(y,r)} u(x) \, dx$, and show that $A(r)$ is constant by differentiating with respect to $r$, and showing that the derivative is zero by using the divergence theorem to replace the integrand with a $\Delta u$.
For the subharmonic and superharmonic cases, the same technique leads to $A(r)$ increasing and decreasing, respectively.