Here's a simple set of examples that are special cases of the family of solutions that I was referring to. (Describing all of the examples would take longer than I have at the moment.)
Let $I\subset\mathbb{R}$ be an open interval and consider metrics on $I\times\mathbb{R}^n$ of the form $$ g = u(t)^2\bigl((\mathrm{d}t)^2 + (\mathrm{d}x^1)^2 + \cdots +(\mathrm{d}x^n)^2\bigr), $$ where $u:I\to\mathbb{R}^+$ is a positive function on $I$. The condition that $g$ have constant scalar curvature is an autonomous nonlinear second order ODE for $u$ that is solvable in terms of Jacobian elliptic functions (i.e., inversion of an elliptic integral). (I won't go into details about this until I have more time. When $n=2$, for example, the equation is $$2\,u(t)\,u''(t)-u'(t)^2+3c\,u(t)^4=0,\tag1$$ where $6c$ is the scalar curvature of $g$.)
For the generic solution $u$, the metric $g$ is not Einstein. In fact, $\mathrm{Ric}(g) = a(t)\,g + b(t)\,\mathrm{d}t^2$ for functions $a$ and $b$ determined in terms of $u$. Generically, $b$ will not be zero. However, since $g$ is conformally flat and has constant scalar curvature, $g$ has harmonic curvature. (This is a well-known special case of harmonic curvature. In dimension $3$, harmonic curvature is equivalent to being conformally flat with constant scalar curvature.)
Now consider a function $f:I\to\mathbb{R}$ regarded as a function on $I\times\mathbb{R}^n$ and ask that $\mathrm{Ric}(g)+\mathrm{Hess}_g(f) = \lambda(t)\,g$ for some $\lambda(t)$. Since $\mathrm{Hess}_g(f) = p(t)\,g + q(t)\,\mathrm{d}t^2$ for some functions $p$ and $q$ determined in terms of $f$ and $u$, this is a single linear second order equation for $f$ involving $u(t)$ that always has solutions on the interval $I$ of definition of $u$. For example, when $n=2$, the equation takes the form
$$
u(t)\,f''(t)-2\,u(t)u'(t)\,f'(t)+2\,u'(t)^2-u(t)u''(t)=0.
\tag2
$$
Thus, $f$ can be found by quadrature once the function $u$ has been specified. For the generic pair $(u,f)$ satisfying (1) and (2), the corresponding $\lambda(t) = \bigl(u'(t)f'(t)-u''(t)\bigr)u(t)^{-3}$ will be nonconstant, thus producing examples of the desired type.
This can be generalized slightly by replacing $\mathbb{R}^n$ and its flat metric by $(N^n,h)$ where $h$ is a metric on $N$ of constant sectional curvature $k$. The equations for $u(t)$ and $f(t)$ change slightly, but they are still solvable using elliptic functions. (In fact, when $k>0$, there are periodic solutions $u(t)$ defined on the entire line $I=\mathbb{R}$ to the equation corresponding to (1), and one can use this to construct examplees defined on compact manifolds.)
In dimension $3$, it turns out that every conformally flat metric $(M^3,g)$ with constant scalar 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)$. It was after I worked that out using an exterior differential system analysis that I realized that the above construction would work in any dimension.
Addendum: It turns out that one can avoid having to solve the equation for $u$ by making a change of variables to eliminate the (non-geometric) parameter $t$. The result is the following: 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. If $I\subset\mathbb{R}^+$ is an interval on the $u$-line where $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$, Moreover, there will be functions $f(u)$ and $\lambda(u)$ such that $\mathrm{Ric}(g) + \mathrm{Hess}_g(f) = \lambda\,g$. For most values of $a$ and $b$, $\lambda(u)$ 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.