gstat::variogram - distance calculation is incorrect

Hello,

I am trying to use 'variogram' from the gstat package to compute an empirical semivariogram for a dataset that spans the North American continent.  I have been unable to get gstat to calculate the great circle distances correctly.  For example, in the example below, the first two lines of vg.foo suggest to me that the calculated distance between foo[1,] and foo[2,] is 4804, and between foo[1,] and foo[3,] is 3047.  I assume those units are km.  

spDistsN1, however tells me that those distances should be 3115 km and 1907 km, respectively.

Could someone suggest to me my error, or how I might use calculate an emprical variogram with correct distances?  I do not think UTM coordinates are an option for me, as my data span approximately 5000 km.

Many thanks,
Tim

=====

code to illustrate the problem:

library(sp)
library(rgdal)
library(gstat)

foo <-
structure(list(z = c(-1.95824831109744, -1.9158901643563, 4.22211761150161, 
3.23356929459598, 1.12038389231868, 0.34613850821113, 1.12589932643631, 
23.517912251617, 3.0519158690268, 3.20261431141517, -2.10947106854739
), lon = c(-125.29228, -82.1556, -98.524722, -99.948333, -104.691741, 
-79.420833, -105.100533, -88.291867, -72.171478, -121.556944, 
-89.34765), lat = c(49.87217, 48.2167, 55.905833, 56.635833, 
53.916264, 39.063333, 48.307883, 40.0061, 42.537756, 44.448889, 
46.242017)), .Names = c("z", "lon", "lat"), row.names = c(NA, 
-11L), class = "data.frame")

coordinates(foo) <- ~lon+lat

proj4string(foo) <- CRS('+proj=longlat')

vg.foo <- variogram(z~1, foo, cloud=TRUE, cutoff=1e10)

cat('==========\nvariogram:\n')
print(head(vg.foo))

cat('==========\nspDistsN1 Distances:\n')
print(spDistsN1(coordinates(foo), coordinates(foo)[1,], longlat=T))