What you did wrong: the function $f$ should be treated as $f = f(\gamma(t)) = f(u(t),v(t))$. So $\frac{d}{dt}f = \partial_u f \frac{d}{dt}u + \partial_vf\frac{d}{dt}v$ using the chain rule. By the parametrization, you $\partial_uf = 0$, while $\partial_v f$ is what you wrote $f'$ originally. (Whereas in the above you implicitly wrote $\frac{d}{dt}f^2 = 2f f'$, which is wrong.)
This is an instance where you let notation get in your way. If would be clearer if you reserve the $\prime$ for $f'$, the derivative of the function $f$ relative to the parameter $v$ and $g'$ for the derivative of the function $f$ relative to parameter $v$, and use explicitly either $\frac{d}{dt}$ or $\cdot$ to denote time derivatives along the geodesic (as the geodesic equations that you originally copied down does).
For deriving the "conservation of energy", it may help if instead of the second equation in the two geodesic equations (also, I think there may be a sign error in it, but I am not 100% certain), you look at that equation multiplied by $(f'^2 + g'^2)\frac{dv}{dt}$, that is
$$ (f'^2 + g'^2) \dot{v} \ddot{v} + ff'\dot{v}(\dot{u})^2 + (f'f'' + g'g'')(\dot{v})^3 = 0 $$
Incidentally, the second thing about angle against the parallels is called "Clairaut's relation".
In my text, it writes a surface of revolution as ${\bf x} = (r(t)cos(\theta),r(t)sin(\theta),z(t))$. Then, the meridians are the $t-curves$, and the circles of latitude are the $\theta - curves$. We get the $t-curves$ by fixing $\theta$, and vice versa.
We know that a curve ${\bf \gamma}$ is a geodesic if its second derivative is everywhere normal to the surface along ${\bf \gamma}$. Thus, what you should do is calculate ${\bf x_1}$ and ${\bf x_2}$, the partial derivatives of the surface with respect to $t$ and $\theta$, and then show that $<{\bf \gamma} '',{\bf x_i}>=0$. This implies that $\gamma ''$ is normal to the surface, since the normal vector ${\bf n} = {\bf x_1}\times{\bf x_2}$ , and thus ${\bf \gamma}$ is a geodesic.
For the $t-curves$, you just have to prove that $<{\bf \gamma} '',{\bf x_i}>=0$.
For the $\theta-curves$, you should calculate $<{\bf \gamma} '',{\bf x_i}>$, and determine what condition would make that innerproduct $0$.
I believe the condition should be that $r'=0$, meaning that the function $r(t)$ is in a local maximum or minimum. This makes sense, as we know that the equator of a sphere is a geodesic, and the equator can be viewed as a $\theta-curve$ for a surface of revolution of a half-circle.
Best Answer
Note $f' = \frac{df}{dv}$, $g' = \frac{dg}{dv}$ and similar for $f'', g''$. Hence (note that $v = v(t)$)
\begin{split} f' f'' + g'g'' &= \frac{1}{2} \frac{d}{dv}\big( (f')^2 + (g')^2\big) \\ \end{split} Hence $$\frac{d}{dv} \sqrt{(f')^2 + (g')^2} = \frac{f'f'' + g'g''}{\sqrt{(f')^2 + (g')^2}}$$ Hence the geodesic equation becomes (write $\ell(t) = \sqrt{(f')^2 + (g')^2}$)
$$ \ell (t) \ddot v(t) + \dot \ell(t) \dot v(t) = 0,$$ or $$ \ell (t) \dot v(t) = C.$$
Write $\ell (t) = \ell(v(t))$ (that is, $\ell(v) = \sqrt{(f'(v))^2 + (g'(v))^2}$), then
$$ t'(v) = \frac{C}{v'(t)} = \frac{C}{\ell(v)}$$ and one can integrate this to obtain $t(v)$, and hence $v(t)$.
Remark of course the above looks familiar: this is the ODE that you solve when you want to find an arc-length parametrization of a curve $\gamma(v) = (f(v), g(v))$. Indeed, since the surface of revolution (as a set) is the same when we reparametrize $(f(v), g(v))$, we may WLOG assume that $\gamma(v)$ is of arc-lenght parametrization to start with, so
$$ 1=|\gamma'(v)| = \sqrt{(f')^2 + (g')^2}.$$ So the geodesic equations becomes
\begin{split} \ddot u + \frac{2f'}{f} \dot u\dot v &=0, \\ \ddot v - ff' (\dot u)^2 &= 0 \end{split}
and when it is clear that $(u, v) = (u_0, at+b)$ are solutions to the geodesic equations.