Note that $\mathrm{grad} f$ acts on the tangent space to a point on $S^2$, so in local coordinates it would be a map from a two-dimensional vector space to itself, whereas what you have calculated is the action of $\mathrm{grad} f$ on the tangent space of a point within $\mathbb{R}^3$. In particular, if you had written $\mathbb{R}^3$ instead of $S^2$ everywhere, then your answer would be correct.
Instead, what's happening is that the integral curves (by definition) are confined to be on $S^2$. Nonetheless, your observation that the curves are moving in the vertical direction at constant speed helps us to cook up a pretty good guess at what the integral curve through a given point looks like.
In particular, we claim that the integral curves are of the form
$$\gamma(t) = (\alpha(t) \cos(u), \alpha(t) \sin(u), z(t))$$
for $u \in [0, 2\pi)$ and $t \in [-1, 1]$. Here, $\alpha(t)$ is whatever function ensures that $\gamma(t)$ is actually a point on the unit sphere. That is, for each $t$, $1 = \lVert \gamma(t) \rVert^2 = \alpha(t)^2 + z(t)^2$, so concretely we have $\alpha(t) = \sqrt{1 - z(t)^2}$. Let's get back to what $\alpha(t)$ and $z(t)$ actually look like in a second.
As you've noted, what we actually want to solve is $\gamma'(t) = (\mathrm{grad} f)_{\gamma(t)}$.
Now, if you define $\hat{f} : \mathbb{R}^3 \to \mathbb{R}$ by $\hat{f}(x, y, z) = z$, then we can use the gradient of $\hat{f}$ to find the gradient of $f$. In particular (see e.g. this answer), $(\mathrm{grad} f)_{\gamma(t)}$ is the orthogonal projection of $(\mathrm{grad} \hat{f})_{\gamma(t)}$ onto $T_{\gamma(t)} S^2$. You've already noticed that $(\mathrm{grad} \hat{f})_{\gamma(t)} = (0, 0, 1)$, so what remains is to find the tangent space. The good thing about having $S^2$ inside $\mathbb{R}^3$, is that for any point $p \in S^2$, $T_p S^2$ is the orthogonal complement to $p$ in $\mathbb{R}^3$, so what we are looking for is the projection of $(0, 0, 1)$ onto the plane orthogonal to $\gamma(t)$. That is, if we write $\gamma(t) = (x(t), y(t), z(t))$, then
$$(\mathrm{grad} f)_{\gamma(t)} = (0, 0, 1) - \frac{(0, 0, 1) \cdot \gamma(t)}{\lVert \gamma(t) \rVert^2} \gamma(t) = (0, 0, 1) - z(t) \gamma(t) = (-x(t)z(t), -y(t)z(t), 1-z(t)^2).$$
Now we can finally get to work. Integrality tells us that
$$(-x(t)z(t), -y(t)z(t), 1-z(t)^2) = (x'(t), y'(t), z'(t)),$$
a collection of a couple of differential equations. Here we have a few special cases to get rid of first as we find a couple of trivial solutions, $\gamma(t) = (0, 0, -1)$ and $\gamma(t) = (0, 0, 1)$. Intuitively, these tell us that if we start off at any of the poles, the integral curves will not take us anywhere.
More interestingly, we find that $z(t) = \frac{e^{2t} -1}{e^{2t} + 1}$ solves the equation $1-z(t)^2 = z'(t)$. Using this particular solution, we find that $- \alpha(t) z(t) = \alpha'(t), \alpha(0) = 1,$ is solved by
$$\alpha(t) = \frac{2e^t}{e^{2t} +1}.$$
At this point, every point on $S^2$ is in the image of one of our solutions, so summing everything up, we have that the integral curves are of the form
$$t \mapsto (0, 0, 1), \\ t \mapsto (0, 0, -1), \\ t \mapsto (\alpha(t) \cos(u), \alpha(t) \sin(u), z(t)), \quad u \in [0, 2\pi),$$
where
$$\alpha(t) = \frac{2e^t}{e^{2t} +1}, \quad z(t) = \frac{e^{2t} -1}{e^{2t} + 1}.$$


matplotlib. – fuglede Dec 05 '16 at 08:47Axes3D.quiver. – fuglede Dec 05 '16 at 10:47