Empirical Variogram from multiple realizations

I am wondering if there is a methodology within R to calculate a single
variogram based on multiple modeled realizations of the same process.  I am
performing regressions on residual fields from multiple climate models.  I
would like to estimate a covariance function common to all models.  I've
tried using the variog function in geoR, however it does not work with
multiple values at single coordinates.  Any advice would be appreciated

By "common to all models" I assume you refer to "pooled", or "averaged
over all models". In that case, you can:

library(gstat)
loadMeuse()
v = fit.variogram(variogram(log(zinc)~1,meuse),vgm(1, "Sph", 900, 1))
sim = krige(log(zinc)~1, meuse, meuse.grid, v, nsim=20, nmax=30)
sim.stacked = stack(sim)
coordinates(sim.stacked) = ~x+y
v.pooled = variogram(values~ind, sim.stacked, dX = 0) # wait a while...
plot(v.pooled, v, ylim = c(0, .7))

please look into the dX argument of ?variogram (to get pooling over
different data sets), and look closely what stack() does. Note that in
~ind, in combination with dX ind is not used as predictor, but as
criterium to exclude point pairs coming from different models (those
with different values for ind).

Best regards,

Edzer Pebesma
Institute for Geoinformatics (ifgi), University of M?nster
Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de