Such a harmonic function does not exist. Here's a sketch of a proof. Consider a level set $\Sigma_c$ of $f=c$. The gradient $\nabla f$ is a constant length vector perpendicular to $\Sigma_c$ at each point of $\Sigma_c$ (let's assume $|\nabla f|=1$).
Also, the vector field $\nabla f$ is divergence-free, since $f$ is harmonic, and thus the flow by $\nabla f$ is volume preserving. Now, flow by the vector field $\nabla f$ for time $t$ takes $\Sigma_c$ into $\Sigma_{c+t}$. Since the vectors are constant length, this gives a local orthogonal coordinate system about $\Sigma_c$, which therefore is Fermi coordinates. So the flow lines of $\nabla f$ are geodesics by the Gauss lemma. 
Since the flow is volume preserving and $|\nabla f|=1$, it also preserves the area of $\Sigma_{c+t}$. This implies that $\Sigma_c$ is a minimal surface, since the derivative of the variation of area under orthogonal deformation is the trace of the second fundamental form. One then computes that the principal curvatures of $\Sigma_c$ at each point are $\pm 1$. This is because under the orthogonal flow for time $t$, the second fundamental form at time $t$ is determined uniquely by the second fundamental form at time $0$. The only way that it can remain trace $0$ is if the principal curvatures are $\pm 1$. Thus, by Gauss' equation, $\Sigma_c$ is isometric to $\mathbb{H}^2_{-2}$. However, [by a result of Doug Moore][1] (generalizing a classic result of Hilbert), there is no isometric immersion of $\mathbb{H}^2_{-2}$ into $\mathbb{H}^3_{-1}$. 


  [1]: http://www.ams.org/mathscinet-getitem?mr=305312