**I'm revising my answer to shorten it, since there is a much simpler way to describe these solutions more fully.** 

Let $(N^n,h)$ be a metric of constant sectional curvature $k$ and consider the quadratic form
$$
g = \frac{\mathrm{d}u^2}{k-a\,u^2+ b\,u^{1-n}} + u^2\,h
$$
on $M^{n+1} = \mathbb{R}^+\times N$, where $a$ and $b$ are constants and $u>0$ is the coordinate on $\mathbb{R}^+$.  If $I\subset\mathbb{R}^+$ is an interval on the $u$-line on which $k-a\,u^2+ b\,u^{1-n} >0$, then $g$ is a Riemannian metric on $I\times N$ that is conformally flat and has constant scalar curvature $S = n(n{+}1)a$.  Hence it has harmonic curvature.  The Ricci curvature is
$$
\mathrm{Ric}(g) = \bigl(n\,a - \tfrac{1}{2}\,b\,u^{-n-1}\bigr)\,g 
+ \frac{(n^2{-}1)b\,\mathrm{d}u^2}{2\bigl(b\,u^2+k\,u^{n+1}-a\,u^{n+3}\bigr)},
$$
so $g$ is Einstein if and only if $b=0$.


Moreover, it is now easy to construct (by quadrature) a function $f = f(u)$ on $I$  such that $\mathrm{Ric}(g) + \mathrm{Hess}_g(f) = \lambda\,g$ for some function $\lambda$.  When $b\not=0$, $\lambda$ will not be constant.  Thus, this gives a completely explicit $3$-parameter family of non-trivial almost Ricci solitons with harmonic curvature.

If $I = (r_1,r_2)$ where $r_2>r_1>0$ are simple roots of $k-a\,u^2+ b\,u^{1-n}=0$, then the curve $v^2 = k-a\,u^2+ b\,u^{1-n}$ in the $uv$-plane has a smooth circle component $C$ between the lines $u=r_1$ and $u=r_2$.  In this case, the metric $g$ extends to a smooth complete metric on $C\times N$ (assuming that $(N,h)$ is complete).  In this way, one can construct many complete or compact examples of such metrics.  However, when $b\not=0$, the functions $f$ and $\lambda$ will only be locally defined unless one passes to the simply-connected cover of $C$, so that $M = \mathbb{R}\times N$.  On this covering space, $f$ (and $\lambda$) can be globally defined.

*Remark:* In dimension $3$, it turns out that *every* conformally flat metric $(M^3,g)$ with constant scalar curvature (i.e., every metric in dimension $3$ with harmonic curvature) that admits a 'Ricci potential', i.e., a function $f$ such that $\mathrm{Ric}(g)+\mathrm{Hess}_g(f) = \lambda\,g$ for some function $\lambda$, is locally of the above form for some $(N^2,h)$ with constant curvature.  It was after I worked out this fact, via an exterior differential system analysis, that I realized that the above construction could be used to produce examples in any dimension.  

In dimensions above $3$, it is not likely that every metric with harmonic curvature that admits a Ricci potential in the above sense is of the form given above, but I don't know a classification, even when $n=4$.