There are plenty of solutions that are not flat. In fact, the general solution (up to diffeomorphism) depends on two functions of three variables. Here is how one can write them down:
First, since $p$ and $h$ are nonzero, we can write $f$ in the form $f = -(p/h)x$ where $x$ is now a function on $M$ that satisfies the system
$$
\Delta x + p\,x^2 = |\nabla x|^2 - x^3 = 0.\tag1
$$
Conversely, if $x$ satisfies this system, then $f = -(p/h)x$ satisfies the original system.
Now, applying a moving frame analysis to the system $(1)$, it is easy to show that any point in the open set $U\subset M$ on which $x$ is nonzero (and therefore positive) has a neighborhood $V$ on which there exist coordinates $(x,y,z)$ and functions $E(x,y,z)$, $F(x,y,z)$, and $G(x,y,z)$ that satisfy $EG-F^2 = 1$ such that $g$ on $V$ can be written in the form
$$
g = \frac{{\mathrm{d}x}^2}{x^3} + x^{p-3/2}\bigl(E\,\mathrm{d}y^2 + 2F\,\mathrm{d}y\,\mathrm{d}z + G\,\mathrm{d}z^2\bigr).
$$
(The coordinates $(x,y,z)$ are not unique but $x$ is (of course) and $\mathrm{d}y\wedge\mathrm{d}z$ is uniquely defined up to sign, since it is the area form associated to the (degenerate) quadratic form $x^{-(p+3/2)}\bigl(x^3\,g - \mathrm{d}x^2\bigr)$. Thus, the ambiguity in the choice of coordinates is essentially one function of two variables.)
Conversely, for any $E$, $F$, and $G$ on a domain in $\mathbb{R}^+\times \mathbb{R}^2$ satisfying $EG-F^2 = 1$, the above formula gives a metric that will have the property that the pair $(g,x)$ satisfy $(1)$.
Addendum: In response to a question, I realized that there is a simpler way to derive the above local normal form than using moving frames, so here it is. Assume that $(g,x)$ are a metric and positive function on $M^3$ satisfying $(1)$. The second equation of $(1)$ implies that the $1$-form $\omega_1 = x^{-3/2}\,\mathrm{d}x$ has $g$-norm equal to $1$, so, locally, we can write $g = {\omega_1}^2+{\omega_2}^2+{\omega_3}^2$, where the $\omega_i$ are an orthonormal coframing. For simplicity, fix the orientation of $M$ so that the $\omega_i$ are an oriented orthonormal coframing, and let $\ast_g$ denote the Hodge star operator of $g$ in this orientation. Thus, $\ast_g\omega_1 = \omega_2\wedge\omega_3$ and $\ast_g1 = \omega_1\wedge\omega_2\wedge\omega_3$ is the volume form. Then, by the definition of the Laplacian with respect to $g$, $\mathrm{d}(\ast_g\mathrm{d}x) = -\Delta x\,\,{\ast_g}1$.
Using the equations $(1)$, we have $\mathrm{d} x = x^{3/2}\,\omega_1$, so
$$
\mathrm{d}(x^{3/2}\,\omega_2{\wedge}\omega_3) = \mathrm{d}\bigl(x^{3/2}\,\,{\ast_g}\omega_1\bigr) = \mathrm{d}(\ast_g\mathrm{d}x) = -\Delta x\,\,{\ast_g}1 \\ = px^2\,\omega_1{\wedge}\omega_2{\wedge}\omega_3 = p \frac{\mathrm{d}x}{x} \wedge \bigl(x^{3/2}\,\omega_2{\wedge}\omega_3\bigr).
$$
Consequently, $\mathrm{d}(x^{3/2-p}\,\omega_2{\wedge}\omega_3)=0$. Thus, locally, one can write
$$
x^{3/2-p}\,\omega_2{\wedge}\omega_3 = \mathrm{d}y\wedge\mathrm{d}z\tag2
$$
for some functions $y$ and $z$. Since $\omega_1\wedge\omega_2\wedge\omega_3= x^{p-3}\, \mathrm{d}x\wedge\mathrm{d}y\wedge\mathrm{d}z$ does not vanish, it follows that $x$, $y$, and $z$ are local coordinates on $M$. Finally, $(2)$ implies that
$$
{\omega_2}^2+{\omega_3}^2 = x^{p-3/2}\bigl(E\,\mathrm{d}y^2 + 2F\,\mathrm{d}y\,\mathrm{d}z + G\,\mathrm{d}z^2\bigr)
$$
for some functions $E$, $F$, and $G$ on the domain of the coordinate
chart $(x,y,z)$ that satisfy $EG-F^2 = 1$.
Note that the same technique allows one to write down the solution pairs $(g,x)$ to $(1)$ in any dimension at least two (not just three).
