Find all entire functions $f$ such that $|f|$ is harmonic.

You have demonstrated that constant functions are the only entire functions for which $|f|^2$ is harmonic.

Now let us investigate the case where $f$ is entire and $ |f|$ is harmonic. We can assume that $f$ is not identically zero.

It is surely possible to calculate $\Delta |f|$ directly, but one can simplify the work by showing first that $f$ has no zeros, so that $f=e^g$ for some entire function $g$: Assume that $f(z_0) = 0$, then $$ 0 = |f(z_0)| = \frac{1}{2\pi} \int_0^{2\pi} |f(z_0+re^{it}| \, dt $$ for all $r > 0$, which implies that $f$ is identically zero, in contrast to the assumption.

So we have $|f| =|e^g| = e^{\operatorname{Re} g} = e^v$ where $v = \operatorname{Re} g$ is harmonic. Then $$ 0 = \Delta |f| = \Delta(e^v) = e^v \left(v_{xx} + (v_x)^2 + v_{yy} + (v_y)^2\right) = e^v \left( (v_x)^2 + (v_y)^2\right) $$ so that $v_x$ and $v_y$ are identically zero. It follows that $v$ is constant. Consequently, $|f|$ is constant, which implies that $f$ is constant.

$f(z)=f(x+iy)=u+iv$. Then from analyticity of $f$: $$ u_x=-v_y\textrm{ and }u_y=v_x.\tag 1 $$ Hence if $|f|=\sqrt{u^2+v^2}$ (away from the zeros of $f$), then $$ |f|_x=\frac{u_x+v_x}{\sqrt{u^2+v^2}}\Rightarrow |f|_{xx}=\frac{(u_{xx}+v_{xx})|f|-\frac{(u_x+v_x)^2}{|f|}}{|f|^2}\Rightarrow $$ $$ |f|_{xx}=\frac{u_{xx}+v_{xx}}{|f|}-\frac{(u_x+v_x)^2}{|f|^3}.\tag 2 $$ In the same way $$ |f|_{yy}=\frac{u_{yy}+v_{yy}}{|f|}-\frac{(u_y+v_y)^2}{|f|^3}.\tag 3 $$ Hence using $u_{xx}+v_{yy}=0$, $v_{xx}+v_{yy}=0$ and equations (1), we get $$ \Delta|f|=|f|_{xx}+|f|_{yy}=-\frac{1}{|f|^3}\left((u_x+u_y)^2+(u_y-u_x)^2\right)=-2|f|^{-3}\left((u_x)^2+(u_y)^2\right).\tag 4 $$ Hence $$ \Delta|f|=0\textrm{ iff }(u_x=u_y=v_x=v_y=0)\textrm{ iff }f=const.\tag 5 $$ In case that $f(z)$ is zero in a set $A$, then $A$ will be discrete and from (5) constant in $\textbf{C}-A$. Hence $f(z)$ zero everywhere.