  For positive $p$-harmonic functions on Riemannian manifolds, we derive a
gradient estimate and Harnack inequality with constants depending only on the
lower bound of the Ricci curvature, the dimension $n$, $p$ and the radius of
the ball on which the function is defined. Our approach is based on a careful
application of the Moser iteration technique and is different from Cheng-Yau's
method employed by Kostchwar and Ni, in which a gradient estimate for positive
$p$-harmonic functions is derived under the assumption that the sectional
curvature is bounded from below.
