Problem with universal kriging using gstat

data.grid <- expand.grid(x=seq(xLL,xUL,xN), y=seq(xLL,xUL,xN))
names(data.grid) = c("X","Y")
gridded(data.grid) <- ~X+Y

memory.size()

memory.size(TRUE)

memory.limit(size=2048)

round(memory.limit()/1048576.0, 2)

options(scipen=3)

options(digits=15)
rm(list=ls())
require(sp)

require(gstat)

require(maptools)

require(foreign)

cellsize <- 10

#                 1           2               3             4
VarName=c("logit(SOM)","logit(TN)","logit(OLSENP)","logit(EXTK)")
out=paste("Variable","No","Model","C0","C","Spatial","Range","R2","p","SSE",sep=",")
setwd("E:\\Rcode\\EJSS\\Hongtong\\village")
data <- read.dbf("libaocun_gs.dbf")
names(data)

class(data)

names(data)[14] <- "X"
names(data)[15] <- "Y"

data$X <- data$X  + 1750
data$Y <- data$Y +  650
 
coordinates(data) <- ~X+Y
proj4string(data) <- CRS("+proj=tmerc +lat_0=0 +lon_0=111 +k=1.0 +x_0=500000 +y_0=0 +ellps=krass

#read in a shape file of the boundary
aSHAPE <- readShapePoly("libao_bd.shp", IDvar="BSM", proj4string=CRS(proj4string(data))) 

#generate the grid coordinates (in meters)
xLL <- round(slot(aSHAPE,"bbox")[1]) - 100
xUR <- round(slot(aSHAPE,"bbox")[3]) + 100
yLL <- round(slot(aSHAPE,"bbox")[2]) - 100
yUR <- round(slot(aSHAPE,"bbox")[4]) + 100
xN <- round((xUR-xLL)/cellsize)
yN <- round((yUR-yLL)/cellsize)
data.grid = SpatialGrid(GridTopology(c(xLL,yLL),c(cellsize,cellsize),c(xN,yN)))
proj4string(data.grid) <- CRS(proj4string(aSHAPE))


#data logit transformation
SOM_MAX      <- max(data$SOM)*1.1
SOM_MIN      <- min(data$SOM)*0.9
SOM_DIFF     <- SOM_MAX - SOM_MIN
data$SOMt    <- log(((data$SOM-SOM_MIN)/SOM_DIFF)/(1-(data$SOM-SOM_MIN)/SOM_DIFF))

TN_MAX      <- max(data$TN)*1.1
TN_MIN      <- min(data$TN)*0.9
TN_DIFF     <- TN_MAX - TN_MIN
data$TNt    <- log(((data$TN-TN_MIN)/TN_DIFF)/(1-(data$TN-TN_MIN)/TN_DIFF))

OLSENP_MAX      <- max(data$OLSENP)*1.1
OLSENP_MIN      <- min(data$OLSENP)*0.9
OLSENP_DIFF     <- OLSENP_MAX - OLSENP_MIN
data$OLSENPt    <-

EXTK_MAX      <- max(data$EXTK)*1.1
EXTK_MIN      <- min(data$EXTK)*0.9
EXTK_DIFF     <- EXTK_MAX - EXTK_MIN
data$EXTKt    <- log(((data$EXTK-EXTK_MIN)/EXTK_DIFF)/(1-(data$EXTK-EXTK_MIN)/EXTK_DIFF))

#SOM
   # plot
   X11(width=8, height=6.5)
   par(mfrow=c(1,1))
   spplot(aSHAPE["BSM"], scales=list(draw=TRUE, cex=0.7), xlab="Easting (m)", ylab="Northing (m)",

   X11(width=8, height=6.5)
   par(mfrow=c(1,2))
   hist(data$SOMt, main="", xlab=VarName[1], n=25)
   boxplot(data$SOMt, xlab=VarName[1])

   #estimate semivariogram model
   soil.v <- variogram(SOMt~X+Y+X*Y+X^2+Y^2, data=data, cutoff=1050, width=90)
   vmodel <- vgm(0.5, "Sph", 1000, 0.1)
   v.sph <- fit.variogram(object=soil.v, model=vmodel, fit.sills=TRUE, fit.ranges=TRUE,

   v.sph

   X11(width=8, height=6.5)
   par(mfrow=c(1,1))
   plot(soil.v, model=v.sph, col="red", cex=1, lwd=2.0, main=VarName[1], xlab="Separation distance

#Univeral kriging
data.gstat <- gstat(id="SOMt", formula=SOMt~X+Y+X*Y+X^2+Y^2, data=data, nmax=15, model=v.sph)

#Predicting surface
data.grid <- predict(data.gstat, data.grid)