**NB:** *For some time, I have been meaning to revise this answer to make it more complete (and, to be frank, more accurate), but I never found the time to do it. The main point is that my original answer did not take into account the difference between the cut locus and the conjugate locus, and, of course, this affects the formula for the distance between points.*

I'm aware of a few metrics with non-constant curvature for which one can write the distance function explicitly in terms of the coordinates. The simplest such metric I know is the (incomplete) metric $ds^2 = y\ (dx^2+dy^2)$ on the upper half plane $y>0$. The Gauss curvature of this metric is $K = 1/(2y^3)>0$, so it's not constant.

Every geodesic of this metric in the upper half plane can be parametrized in the form
$$
x = a + b\ t\qquad\qquad y = b^2 + \frac{t^2}{4}
$$
for some constants $a$ and $b$, and, for such a geodesic, the arclength function along the curve is
$$
s = c + b^2\ t + \frac{t^3}{12}\ .
$$
for some constant $c$.

Using these formulae, one finds that two points $(x_1,y_1)$ and $(x_2,y_2)$ are joinable by a geodesic segment if and only if $4y_1y_2 \ge (x_1{-}x_2)^2$. In the case of strict inequality, there are two geodesic segments joining the two points, and the length of the shorter segment is
$$
L_1\bigl((x_1,y_1),(x_2,y_2)\bigr)
= {1\over3}\sqrt{3(x_1{-}x_2)^2(y_1{+}y_2)+4(y_1^3{+}y_2^3)
- (4y_1y_2-(x_1{-}x_2)^2)^{3/2}}\ .
$$
Note that, in a sense, this is better than the constant curvature case. Here, the distance function is algebraic in suitable coordinates, whereas, in the constant nonzero curvature cases, the distance function is not.

*However,* the function $L_1$ does not necessarily give the actual distance between the two points (i.e., the infinimum of the lengths of curves joining the two points), and it's not only because not every pair of points can be joined by a geodesic. To see this, one should complete the upper half plane by adding a point that represents the 'boundary' $y=0$. The Riemannian metric does not extend smoothly across this 'point', of course (after all, the Gauss curvature blows up at you approach this point), but it does extend as a metric space. The vertical lines, which are geodesics, can then be used to join $(x_1,y_1)$ to $(x_2,y_2)$ by going through the singular point, and the total length of this geodesic is
$$
L_2\bigl((x_1,y_1),(x_2,y_2)\bigr)
= \frac{2}{3}\bigl({y_1}^{3/2}+{y_2}^{3/2}\bigr).
$$
(Also, note that $L_2$ is defined for *any* pair of points in the upper half-plane.)
If one doesn't like this path that goes through the singular point, one can easily perturb it slightly to avoid the singular point and not increase the length by much, so it's clear that the infimum of lengths of curves lying strictly in the upper half plane and joining the two points is no more than $L_2$.

This suggests that the true distance function $L$ should be the minimum of $L_1$ and $L_2$ where they are both defined, i.e., where $4y_1y_2 \ge (x_1{-}x_2)^2$, and $L_2$ on the set where $4y_1y_2 < (x_1{-}x_2)^2$.

To get a sense of how these two formulae interact, one can use the fact that $x$-translation preserves the metric while the scalings $(x,y)\mapsto (ax,ay)$ for $a>0$ preserve the metric up to a homothety (and hence preserve the geodesics and scale the distances). These two actions generate a transitive group on the upper half plane, so, it suffices to see how these two functions interact when $(x_1,y_1) = (0,1)$, i.e., to see the conjugate locus and cut locus of this point.

The conjugate locus is easy: It's just $y-x^2/4=0$, which is the boundary of the region $y-x^2/4\ge0$ consisting of the points that can be joined to $(0,1)$ by a geodesic segment. Meanwhile, the cut locus is given by points $(x,y)$ that satisfy $y-x^2/4\ge0$ and for which $L_1\bigl((0,1),(x,y)\bigr) = L_2\bigl((0,1),(x,y)\bigr)$. In fact, one has $L_1\bigl((0,1),(x,y)\bigr) < L_2\bigl((0,1),(x,y)\bigr)$ only when $y > f(x)$, where $f$ is a certain even algebraic function of $x$ that satisfies $f(x) \ge x^2/4$ (with equality only when $x=0$). Moreover, for $|x|$ small, one has
$$
f(x) = \left({\frac{{\sqrt{3}}}{4}}x\right)^{4/3} + O(x^2)
$$
while, for $|x|$ large, one has
$$
f(x) = \left({\frac{\sqrt{3}}{4}}x\right)^{4} + o(x^4).
$$

Thus, all of the geodesics leaving $(x,y)=(0,1)$, other than the vertical ones, meet the cut locus *before* they reach the conjugate locus (and they all do meet the conjugate locus).

Thus, the actual distance function for this metric is explicit (it's essentially the minimum of $L_1$ and $L_2$), but it is only semi-algebraic.

*Remark:* The thing that makes this work is that, while the metric has only a 1-parameter family of symmetries, it has a 2-parameter family of homotheties (as described above), and this extra symmetry of the geodesics is critical for making this work. Of course, there are other such metrics, all the ones of the form $ds^2 = y^{a}\ (dx^2+dy^2)$ ($a$ is a constant) have this property and don't have constant curvature unless $a = 0$ or $a = -2$. You don't get algebraic answers for all values of $a$, of course, but there is a way to get $D$ implicitly defined in terms of a special function (depending on the value of $a$).

More generally, the metrics whose geodesics admit more symmetries than the metric itself does tend to have such formulae. I'm not aware of any other cases in which one can get $D$ so explicitly.