I am interested in a reference for (or counterexample to) a particular $C^0$ estimate for solutions of the Laplace equation with Neumann boundary conditions. More precisely, let $(M,g)$ be a $C^\infty$ compact Riemannian manifold with non-empty smooth boundary $\partial M$. Let $\Delta_g=-div_g\nabla$ denote the Laplace operator and $\nu$ the outward pointing unit normal on $\partial M$.

Consider (strong) solutions $u$ to the following system \begin{cases} \Delta_gu=f&\text{ in }M\\ \partial_\nu u=h&\text{ on }\partial M \end{cases} where $f\in C^\infty(M)$ and $h\in C^\infty(\partial M)$ satisfy $\int_M fdvol_g=\int_{\partial M}hd\sigma_g$ (the smoothness condition is for simplicity but the later condition is necessary for the problem to be well-posed). Here $dvol_g$ denotes the Riemannian volume form and $d\sigma_g$ denotes its pullback under the inclusion $\partial M\to M$. To specify $u$ uniquely, let us also suppose that it has zero average value i.e. $\int_Mudvol_g=0$. My question is this: does there exist $C>0$, independent of $f$ and $h$, such that $\sup_M|u|\leq C(\sup_M|f|+\sup_{\partial M}|h|)$?

Most classical sources (Evans, Gilbarg+Trudinger, and Aubin) either only consider Dirichlet conditions (in which case the analogous estimate is essentially just the maximum principle) or only provide Schauder estimates for the Neumann problem -- this is not quite what I want.

Does anyone know of a resource for this kind of estimate or of a counterexample?