Actually there's another conversion when you reuse an xarray dimension coordinate in array-like computations:

```python ds = xr.Dataset(coords={"x": np.array([1.2, 1.3, 1.4], dtype=np.float16)})

coordinate data is a wrapper around a pandas.Index object

(it keeps track of the original array dtype)

ds.variables["x"]._data

PandasIndexingAdapter(array=Float64Index([1.2001953125, 1.2998046875, 1.400390625], dtype='float64', name='x'), dtype=dtype('float16'))

This coerces the pandas.Index back as a numpy array

np.asarray(ds.x)

array([1.2, 1.3, 1.4], dtype=float16)

which is equivalent to

ds.variables["x"]._data.array()

array([1.2, 1.3, 1.4], dtype=float16)

```

The round-trip conversion preserves the original dtype so different execution times may be expected.

I can't tell much why the results are different (how much are they different?), but I wouldn't be surprised if it's caused by rounding errors accumulated through the computation of a complex formula like haversine.