2/3d interpolation from a regular grid to another regular grid

I just read the description in ?Krig in the package fields which  
says:
" Fits a surface to irregularly spaced data. "

Yes, that is the most general case. Regular data location is a subset
of irregular. Anyway, kriging, just one g, after the name of Danie
Krige, the south african statistician who first applied such method
for minig survey.

My problem is simpler

...

So it is really purely numerical.

...

I just hoped that R had that already coded ...

Of course R has ... ;) If your grids are really as simple as the
example you posted above, and you have a really little variability,
all you need is a "moving average", the arithmetic mean of the two
nearest points belonging to grid1 and grid2 respectively. I assume
that your regularly shaped grids are values stored in matrix objects.

The functions comes from the "diff.default" code (downloading the R
source code, I assure, is worth):

my.interp <- function(x, lag = 1)
{
    r <- unclass(x)  # don't want class-specific subset methods
    i1 <- -1:-lag
    r <- (r[i1] + r[-length(r):-(length(r)-lag+1)])/2
    class(r) <- oldClass(x)
    return(r)
}

Finally,

g1 <- apply(grid1val,1,my.interp)
g2 <- apply(grid2val,2,my.interp)

give the interpolations on gridFinal, provided that all gridFinal
points are within the grid1 and grid2 ones.

If you want the mean from 4 points, you apply once more with lag=3,
cbind/rbind to the result columns/rows o NAs, and you calculate the
mean of the points of the two matrixes.
This is the simplest (and quickest) moving average that you can do.
For more complicated examples, and for 3d, you have to go a little
further, but the principle holds.