PROBLEM: Suppose $f$ is a smooth non-constant function on a connected Riemann manifold $M$ of dimension 2 such that $f$ is bounded and $\Delta_M f \ge0$. Show $f$ cannot attain its supremum.
I try to reduce the problem to the special case $M=\mathbb R^2$. In this case we can use mean value property.
What's more, can this problem be generalized to any connected Riemann manifold $M$ of an arbitrary dimension n?
Suppose the subharmonic function $u:M\to\mathbb R$ attained its supremum at $x\in M$. Take a small neighborhood $U$ of $x$ homeomorphic to a disc. Then by the uniformization theorem $U$ is conformally equivalent to the unit disc. Let $\phi:D\to U$ be the conformal map. The Laplacian is conformally invariant in two dimensions, so $\Delta u\geq0$ on $U$ implies that $\Delta (u\circ\phi)\geq0$ on the disc $D$. But in the Euclidean plane we have the maximum principle for subharmonic functions, so $u$ has to be constant in $U$.
This shows that the set where $u$ reaches its supremum is open, and by continuity of $u$ it is also closed. Therefore if $u$ reaches its supremum, it has to be constant on the connected component.
For higher dimensions the uniformization theorem cannot be used, and a similar reduction to the Euclidean case is not possible. Something more clever will be needed, but I suspect that the result is still true.
Added later: Let me elaborate a bit on conformal transforms and the Laplacian. We have a conformal map $\phi:D\to U$ and a function $u:U\to\mathbb R$. We would like to be able to compare $\Delta_D(u\circ\phi)$ and $(\Delta_U u)\circ\phi$, both functions $D:\to\mathbb R$. We know that $(\Delta_U u)\circ\phi\geq0$ and we would like to conclude that $\Delta_D(u\circ\phi)\geq0$.
Conformality means that there is a function $\eta:D\to\mathbb R$ so that $g_D=e^{2\eta}\phi^*g_U$. Here $g_D$ and $g_U$ are the Riemannian metrics on $D$ and $U$ ($g_D$ is the Euclidean metric) and $\phi^*$ is the pullback over $\phi$. Many objects behave relatively nicely under conformal changes, as documented in Wikipedia (and other places, of course). In particular, in dimension two the Laplace-Beltrami operators related to these two metrics are related so that $\Delta_D(u\circ\phi)=e^{-2\eta}(\Delta_U u)\circ\phi$. (In the Wikipedia article, find the formula for $\tilde\Delta f$ and notice that it contains an unwanted term with a coefficient $(n-2)$.) This gives us what we wanted.