Harmonic-curvature warped products over surfaces

For warped products with harmonic curvature, nonconstant warping functions $\phi$, and compact two-dimensional bases $(M,h)$, we establish a dichotomy: either the Gaussian curvature $K$ of the metric $g=\phi^{-2}h$ is constant and negative, or $\phi$ equals a specific elementary function of $K$, also depending on the dimension $p$ and Einstein constant $\varepsilon$ of the fibre. In both cases the fibre must be an Einstein manifold with $p>1$ and $\varepsilon>0$, while the function $f=\phi^{p/2}$ satisfies a Yamabe-type second-order differential equation on $(M,g)$. We prove that both possibilities are realized on every closed orientable surface of genus greater than $1$, and in the latter case -- which also occurs on the $2$-sphere and real projective plane -- the metrics in question constitute uncountably many distinct homothety types.