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}$. 

**Addendum:** Given GB's comment, I should add some explanation of why the principal curvatures are $\pm 1$. If we take a region $\Omega \subset \Sigma_c$, and let $A_t$ denote the area of the image of $\Omega$ at time $t$ under the flow, then $\partial A_t/\partial t = \int_\Omega H dA$, where $H$ is the mean curvature. Since $\nabla f$ preserves the area of $\Omega$, we see that $H=0$, so $\Sigma_c$ is a minimal surface. 

The second variation of area is given by $\partial^2 A_t/\partial t^2 = \int_\Omega (2 Det(B)-Ric(\nabla f)) dA =0$, where $B$ is the second fundamental form, and $Ric$ is the Ricci curvature. Since $\Sigma_c$ is minimal, $Det(B)=-p^2$, where $p$ is a principal curvature. Also, $Ric(\nabla f)=-2$, since $|\nabla f|=1$. So we get that $p=1$, and the principal curvatures are $\pm 1$. 

One can also see more directly that the principal curvatures of $\Sigma_t$ evolve like the curvature of plane hyperbolic curves (this can be seen by considering osculating spheres). The only way that the principal curvatures can remain opposite is if they are $\pm 1$. 


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