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.
In fact the conjecture becomes true if you add embedded to the
hypothesis according to a theorem of Colding and Minicozzi.
(Colding, Tobias H.; Minicozzi, William P., II The Calabi-Yau conjectures for embedded surfaces. Ann. of Math. (2) 167 (2008), no. 1, 211–243.)
Who knows what Hadamard actually had in mind.
(Sorry this should have been a comment, its certainly not an answer to the
question.)
Best Answer
Yes, you can easily do this with surfaces of revolution: If you take a metric of the form $ds^2 = dr^2 + f(r)^2\,d\theta^2$, where $f$ is an odd function of $r$ satisfying $f'(0)=1$, then the Gauss curvature function $K(r)$ satisfies $$ f''(r) + K(r)\,f(r) = 0,$$ where $K(r)$ is an even function of $r$.
Suppose now that $K_1$ and $K_2$ are any two even functions of $r$ that have a nondegenerate minimum of $0$ at $r=0$. Then there will be a unique smooth, increasing, odd function $\phi$ that satisfies $K_1(r) = K_2(\phi(r))$ (and $\phi(0)=0$). Now consider the mapping of the disk $|r|<\epsilon$ defined by $F(r,\theta) = (\phi(r),\theta)$, and let $f_1(r)$ and $f_2(r)$ be the odd functions of $r$ that satisfy $f'_i(0)=1$ and $$ f''_i(r) + K_i(r)\,f_i(r) = 0.$$ Then $F$ aligns the curvatures of the two metrics $ds_i^2 = dr^2 + f_i(r)^2\,d\theta^2$, but generally isn't an isometry if $\phi$ is not the identity map. Moreover, it is not hard to show that there is no isometry between the two metrics.
The two metrics $ds^2_1$ and $F^*(ds^2_2)$ now satisfy your requirements.
N.B.: In addition to Kulkarni's paper, which the OP cites above, it is probably a good idea to look at Yau's 1974 Annals paper, "Curvature preserving diffeomorphisms" (http://www.jstor.org/stable/1970843), where he discusses why the result does not hold in dimension $3$. I might also add that you could find some useful information by looking at the answers to the MO question Does the curvature determine the metric?