I'm looking for explicit examples of Riemannian surfaces (two-dimensional Riemannian manifolds $(M,g)$) for which the distance function d(x,y) can be given explicitly in terms of local coordinates of x,y, assuming that x and y are sufficiently close.  By "explicit", I mean things like a closed form description in terms of special functions, by implicitly solving a transcendental equation or (at worst) by solving an ODE, as opposed to having to solve a variational problem or a PDE such as the eikonal equation, or an inverse problem for an ODE, or to sum an asymptotic series.

The only examples of this that I know of are the constant curvature surfaces, which can be locally modeled either by the Euclidean plane ${\bf R}^2$, the sphere ${\bf S}^2$, or the hyperbolic plane ${\bf H}^2$, for which we have classical formulae for the distance function.

But I don't know of any other examples.  For instance, the distance functions on the surface of the solid ellipsoid or solid torus in ${\bf R}^3$ look quite unpleasant already to write down explicitly.  Presumably Zoll surfaces would be the next thing to try, but I don't know of any tractable explicit examples of Zoll surfaces that are not already constant curvature.

I'll briefly spell out what others have pointed to concerning geodesics on surfaces of revolution (or more generally, surfaces with a 1-parameter group of symmetries), because it's nice and not as widely understood as it should be.

Geodesics on surfaces of revolution conserve angular momentum about the central axis, so the geodesic flow splits into 2-dimensional surfaces having constant energy (~length) and angular momentum  (The more general principle is that the inner product of the tangent to a geodesic with any infinitesimal isometry of a Riemannian manifold is constant). The surfaces are generically toruses. The shadow of these toruses on the surface of revolution is an annulus, a component of a set of $r \ge r_0$, where on each point with $r > r_0$ there are two vectors having the given angular momentum, but they merge at the boundary, both becoming tangent to the boundary of the annulus.  If you sketch the picture, you will see the torus.  The geodesics correspond to the physical phenomenon of the pattern of string or thread mechnically but passively wound around a cylinder.  As string builds up in the middle, geodesics start to oscillate back and forth in a sinusoidal pattern, further amplifying the bulge in the middle.

To find the geodesic from point x to point y, you need to know which angular momentum will take you from x to y.  For any two meridian circles and any choice of angular momentum, the geodesics of given angular momentum map one circle to the other by a rotation.  Both the angle of rotation of the map and the length of the particular family of geodesics traversing the annulus is given by an integral over an interval cutting across the annulus, since the the slope of the vector field at all intervening points is known. I have an aversion to actual symbolic computation so I won't give you example formulas, but I believe this should meet your criterion for explicitness.

But to take a step back:  this question, asking for an explicit formula, has an unstated (and probably unintended) connotation that is worth examining: this use of language implicitly suggests that non-symbolic forms are less worthy.  I don't know the background motivation for the question, but an alternative question for some purposes would be to give example of surfaces where you can exhibit the distance function.  Communication of mathematics is biassed toward symbolic forms.   However, for many people and many purposes, some kind of graphical representation of the distance function, and/or diagrams or explanations of why it is what it is as well as a striaghtforward method for computing it, would often be better than a symbolic answer.

The geodesic flow of course is an ordinary differential equation. It is a vector field on the 3-manifold of unit-length tangent vectors to the surface, defined by very easy equations: the vectors are tangent to the surface, and their derivative (= the 2nd derivative of a geodesic arc) is normal to the surface. The solutions may not always have a nice symbolic form, but they always have a nice and easy-to-compute geometric form.  Finding the distance involves the implicit function theorem, but this is easy and intuitive.  One could, for instance, easily draw a parametric surface that is the graph of distance as a function of position directly from solutions to the ODE (which no doubt sometimes even have reasonable symbolic representations). Both the ODE for the geodesic flow and the inverse function to give distance as a function of position are easy to compute numerically, and easy to understand qualitatively.
