Let $f(x)$ be a harmonic function in $\mathbb{R}^n$, and let $Z$ be the zero level set $$ Z = \{(x)\in\mathbb{R}^n \mid f(x) = 0\}. $$ Arguing using the Harnack inequality, I believe one can show that $$ \frac{|\nabla f(x)|}{|f(x)|} \lesssim r^{-1} $$ when $d(x,Z) = \inf_{z\in Z} |x-z| \geq r$. I wonder if the opposite inequality holds.
Do there exist constants $c,C>0$ depending only on $n$ such that $$ c\,d(x,Z)^{-1} \leq \frac{|\nabla f(x)|}{|f(x)|} \leq C\, d(x,Z)^{-1}? $$
This sounds like it is too good to be true, so I would be interested to see a counterexample. I am also curious about the special case that $n=2$ and $f(x)$ is the real part of a polynomial.
Consider $f(z) = \operatorname{Re}(z^2+1)$, which in real notation is $f(x, y) = x^2-y^2+1$. At the point $(0, 0)$ the gradient vanishes, while the distance to the zero set is $1$. So, the lower bound does not hold.
The upper bound, as you said, follows from Harnack's inequality. Let $U^+ = \{x:f(x)>0\}$. Consider a Whitney cube $Q\subset U$, i.e., $\operatorname{diam} Q$ is comparable to $\operatorname{dist}(Q, \partial U)$, let's call it $r$. By Harnack, all values of $f$ on $Q$ are of about the same size, say $\sim M$. By the interior gradient bound applied to $Q$, the gradient of $f$ at the center of $Q$ is controlled by $M/r$, as claimed.