diff --git a/R/bk.grid.R b/R/bk.grid.R
index 8fb3404cbc1bf903d0902b8515036be1725cf969..72172f335d7165d5dbcae646efc683a325171f1f 100644
--- a/R/bk.grid.R
+++ b/R/bk.grid.R
@@ -18,7 +18,8 @@ bk.grid <- function(point.obj,
                     mode=3,
                     duplicate = "error",
                     dupfun = NULL,
-                    method="gqr")
+                    method="gqr",
+ 		    get.lm=FALSE)
 {
     tmp   <- check.gridparams(angle,xsw,xne,ysw,yne,
                               dx,dy,nx,ny)
@@ -78,6 +79,15 @@ bk.grid <- function(point.obj,
         
     snbbit<-rep(0,1+n*nz)
     snbbit[1]<-1
+    if(get.lm){
+      mu <- double(ntrend*nz)
+      lambda <- double(n*nz)
+      lambd0 <- double(nz)
+    } else {
+      mu <- double(ntrend)
+      lambda <- double(n)
+      lambd0 <- double(1)
+    }
     
 
     ans<-.C("bk_grid",
@@ -121,10 +131,12 @@ bk.grid <- function(point.obj,
 	    nsmax   = as.integer(nsmax),
             lwork   = as.integer(lwork),
             mode    = as.integer(mode),
-            lambda  = double(n*nz),
-            lambd0  = double(nz),
+            mu      = as.double(mu),
+            lambda  = as.double(lambda),
+            lambd0  = double(lambd0),
             bits    = as.integer(c(integer(nz),snbbit)),
             ierr    = integer(1),
+            get.lm  = as.integer(get.lm),
             glsmth  = as.integer(method),
 #           ,.Package= "baykrig"
             )
@@ -134,10 +146,13 @@ bk.grid <- function(point.obj,
                  z=matrix(ans$zg,nx,ny),
                  var=matrix(ans$varg,nx,ny),
                  done=matrix(ans$dog, nx, ny),
-                 snb=matrix(ans$bits[(nz+2):(nz+n*nz+1)],nrow=n,ncol=nz,byrow=F),
-                 lambda=matrix(ans$lambda,nrow=n,ncol=nz,byrow=FALSE),
-                 lambda0=matrix(ans$lambd0,nx,ny)
+                 snb=matrix(ans$bits[(nz+2):(nz+n*nz+1)],nrow=n,ncol=nz,byrow=F)
                  )
+    if(get.lm){
+      retval$lambda <- matrix(ans$lambda,nrow=n,ncol=nz,byrow=FALSE)
+      retval$lambda0 <- matrix(ans$lambd0,nx,ny)
+      retval$mu <- matrix(ans$mu,nrow=n,ncol=ntrend,byrow=FALSE)
+    }
     retval$z[retval$done<=0] <- NA
     retval$var[retval$done<=0] <- NA
     retval$lambda0[retval$done<=0] <- NA
diff --git a/devel/fresh/.Rhistory b/devel/fresh/.Rhistory
index 2eb5a65d113951b44da4f4982d34ba3e71065d61..ca4c650ac85912b04168f071600f36283c174d85 100644
--- a/devel/fresh/.Rhistory
+++ b/devel/fresh/.Rhistory
@@ -1,32 +1,32 @@
-leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=100, ny=100, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean",method="ols")
+ls()
 library(baykrig)
 library(baykrig)
-leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=5, ny=5, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
 library(baykrig)
-leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=5, ny=5, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
-leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=50, ny=50, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
-leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=150, ny=50, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
 library(baykrig)
-leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=10, ny=10, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
-plot(leman.bk)
-leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=100, ny=100, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
-plot(leman.bk)
-leman.bk.2<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=200, ny=200, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
 library(baykrig)
-leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=10, ny=10, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
-plot(leman.bk)
+leman.bk.pt <- bk.points(point = leman.88.pt , at = "cadpbm", prior=leman.prior, var.mod.obj = leman.88.vmsph, x=c(530,550),y=c(142,147), trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
 library(baykrig)
-leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=50, ny=50, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
-plot(leman.bk)
-plot(leman.bk,show.snb=T)
-plot(leman.bk,show.snb=T,bias=T)
-dev.print(file="bias.ps")
+leman.bk.pt <- bk.points(point = leman.88.pt , at = "cadpbm", prior=leman.prior.1, var.mod.obj = leman.88.vmsph, x=c(530,550),y=c(142,147), trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
+example(prior)
+example(emp.prior)
+ls("package:baykrig")
+example(empirical.prior)
+leman.bk.pt <- bk.points(point = leman.88.pt , at = "cadpbm", prior=leman.prior.1, var.mod.obj = leman.88.vmsph, x=c(530),y=c(142), trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
+leman.bk.pt
+str(leman.bk.pt)
 q()
 library(baykrig)
-plot(leman.bk,show.snb=T,bias=T)
-ls()
+example(bk.grid)
 library(baykrig)
-leman.bk.pt <- bk.points(point = leman.88.pt , at = "cadpbm", prior=leman.prior, var.mod.obj = leman.88.vmsph, x=c(530,550),y=c(142,147), trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
-ls()
+example(empirical.prior)
+save.image()
+help(bk.grid)
 library(baykrig)
-leman.bk.pt <- bk.points(point = leman.88.pt , at = "cadpbm", prior=leman.prior, var.mod.obj = leman.88.vmsph, x=c(530,550),y=c(142,147), trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
+library(baykrig)
+leman.bk1 <- bk.grid(point = leman.88.pt, at = "cadpbm", prior=leman.prior.1,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=100, ny=100, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean",get.lm=T)
+library(baykrig)
+leman.bk1 <- bk.grid(point = leman.88.pt, at = "cadpbm", prior=leman.prior.1,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=100, ny=100, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean",get.lm=F)
+str(leman.bk1)
+str(leman.bk1$var)
+(leman.bk1$var)
+q()
diff --git a/src/bk_grid.c b/src/bk_grid.c
index 3c53cc0fb94498cf8f247380343adc8f0163c71f..640eff12005e05a8937a79aa2533283b6b8107b2 100644
--- a/src/bk_grid.c
+++ b/src/bk_grid.c
@@ -1,6 +1,6 @@
 #include "bk_grid.h" 
 
-void F77_NAME(bk_grid)(double *xsw,
+void bk_grid__(double *xsw,
 	     double *ysw,
 	     double *xne,
 	     double *yne,
@@ -40,10 +40,12 @@ void F77_NAME(bk_grid)(double *xsw,
 	     int *nsmax,
 	     int *lwork,
 	     int *mode,
+	       double *mu,
 	     double *lambda,
 	     double *lambd0,
 	     int *searchnb,
 	     int *ierr,
+	       int *retlm,
 	     int *glsmth){
 /* simple Fortran wrapper */
 bk_grid(xsw,
@@ -86,10 +88,12 @@ bk_grid(xsw,
 	     nsmax,
 	     lwork,
 	     mode,
+	mu,
 	     lambda,
 	     lambd0,
 	     searchnb,
 	     ierr,
+	retlm,
 	     glsmth);
 }
 
@@ -133,10 +137,12 @@ void bk_grid(double *xsw,
 	     int *nsmax,
 	     int *lwork,
 	     int *mode,
+	     double *mu,
 	     double *lambda,
 	     double *lambd0,
 	     int *searchnb,
 	     int *ierr,
+	     int *retlm,
 	     int *glsmth){
     
     int nz=(*nx)*(*ny), ldc0=(*n),ldphwk=(*ntrend),
@@ -148,7 +154,6 @@ void bk_grid(double *xsw,
 	*fwork, *fwrk2, *f0work, *dist, *kwork,
 	*rhswork, *fpwork, *fpfwork, *fpf0wrk, *chlup, 
 	*cminv, *work, *ferr, *berr,
-	*mu, 
 	cov0;
     int *indsnb, *indsnw, *indsrt, *ipiv, *ipvt, *iwork;
 
@@ -183,7 +188,6 @@ void bk_grid(double *xsw,
     ferr   =Calloc((size_t)(*n),double);
     berr   =Calloc((size_t)(*n),double);
     iwork  =Calloc((size_t)(3*(*n)),int);
-    mu     =Calloc((size_t)(*ntrend),double);
 #else
     c0vec  =(double *) R_alloc((*n),sizeof(double));
     muwrk  =(double *) R_alloc((*ntrend),sizeof(double));
@@ -213,7 +217,6 @@ void bk_grid(double *xsw,
     ferr   =(double *) R_alloc((*n),sizeof(double)); 
     berr   =(double *) R_alloc((*n),sizeof(double)); 
     iwork  =(int *) R_alloc(3*(*n),sizeof(int)); 
-    mu     =(double *) R_alloc((*ntrend),sizeof(double));
 #endif
 
     F77_CALL(bkgrid)(xsw,
@@ -301,10 +304,10 @@ void bk_grid(double *xsw,
 		     lambd0,
 		     searchnb,
 		     ierr,
+		     retlm,
 		     glsmth);
 
 #ifndef TRANSIENT
-    Free(mu);
     Free(iwork);
     Free(berr);
     Free(ferr);
@@ -338,7 +341,6 @@ void bk_grid(double *xsw,
 
 #else
     /*
-    free(mu);
     free(iwork);
     free(berr);
     free(ferr);
diff --git a/src/bk_grid.h b/src/bk_grid.h
index 0ea014e73027447deb33786b133a24f4b9675a13..5a584ddd4b157b16777446c75f79e314493e7156 100644
--- a/src/bk_grid.h
+++ b/src/bk_grid.h
@@ -43,13 +43,15 @@ void bk_grid(double *xsw,
 	     int *nsmax,
 	     int *lwork,
 	     int *mode,
+	     double *mu,
 	     double *lambda,
 	     double *lambd0,
 	     int *searchnb,
 	     int *ierr,
+	     int *retlm,
 	     int *glsmth);
 
-void F77_NAME(bk_grid)(double *xsw,
+void bk_grid__(double *xsw,
 	     double *ysw,
 	     double *xne,
 	     double *yne,
@@ -89,10 +91,12 @@ void F77_NAME(bk_grid)(double *xsw,
 	     int *nsmax,
 	     int *lwork,
 	     int *mode,
+	       double *mu,
 	     double *lambda,
 	     double *lambd0,
 	     int *searchnb,
 	     int *ierr,
+	       int *retlm,
 	     int *glsmth);
 	     
 void F77_NAME(bkgrid)(double *xsw,
@@ -180,6 +184,7 @@ void F77_NAME(bkgrid)(double *xsw,
 		      double *lambd0,
 		      int *searchnb,
 		      int *ierr,
+		      int *retlm,
 		      int *glsmth);
     
 
diff --git a/src/bkgrid.f b/src/bkgrid.f
index 908cc737edfe36a2b1f1e153d2c4952394a47218..1a804e7a9e258c8a59b14abb0c7ddcd839ab4641 100644
--- a/src/bkgrid.f
+++ b/src/bkgrid.f
@@ -83,9 +83,10 @@
      .                  LAMBD0,
      .                  BITS,
      .                  IERR,
+     .                  RETLM,
      .                  GLSMTH)
       IMPLICIT NONE
-      INTEGER NX,NY,NZ,N,COVTYPE,TREND,NTREND,LDZG,
+      INTEGER NX,NY,NZ,N,COVTYPE,TREND,NTREND,LDZG,RETLM,
      .        NSEARCH,NSMIN,NSMAX,MODE,IERR,INDSNB(*),INDSNW(*),
      .        INDSRT(*),IPIV(*),IPT,EXTRAP,DOG(NX,*),LDKWRK,EXTCOV,
      .        LDCOV,LDFWRK,LDMPR,LDPHPR,LDPHWK,LDCVBT,LWORK,LDLMBD,
@@ -94,7 +95,7 @@
      .                 XG(*),YG(*),ZG(LDZG,*),VARG(LDZG,*),
      .                 LON(*),LAT(*),Z(*),COVMAT(LDCOV,*),C0VEC(*),COV0,
      .                 RSEARCH,FWORK(LDFWRK,*),F0WORK(*),
-     .                 KWORK(LDKWRK,*),RHSWORK(*),MU(*),
+     .                 KWORK(LDKWRK,*),RHSWORK(*),MU(NTREND,*),
      .                 LAMBDA(LDLMBD,*),FWRK2(LDFWRK,*),
      .                 COVPAR(*),
      .                 FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*),
@@ -255,10 +256,16 @@ c     loop over all points and pass them to KRIGE:
       if (bits(1+nz).ne.0) then
          usesbbt=1 
       end if
-      PCNT=0
+      if (retlm.eq.1) then
+         PCNT=0
+      else
+         PCNT=1
+      end if
       DO 20 I=1,NX
          DO 10 J=1,NY
-            PCNT=PCNT+1
+            if (retlm.eq.1) then
+               PCNT=PCNT+1
+            end if
             BITS(I+NY*(J-1)) = PCNT
             X0=XG(I)
             Y0=YG(J)
@@ -335,7 +342,7 @@ c           the main work is now done by BK:
      .              IPIV,
      .              IWORK,
      .              MODE,
-     .              MU,
+     .              MU(1,pcnt),
      .              NTREND,
      .              Z0,
      .              NA0,
@@ -359,12 +366,16 @@ c           extract results for this tile
                   DOG(I,J)=-1
                   ZG(I,J)=0
                   VARG(I,J)=0
-                  LAMBD0(I,J)=0
+                  if (retlm.eq.1) then
+                     LAMBD0(I,J)=0
+                  end if
 c                  write(*,*)"x"
                ELSE
                   ZG(I,J)=Z0
                   VARG(I,J)=VAR0
-                  LAMBD0(I,J)=L0
+                  if (retlm.eq.1) then
+                     LAMBD0(I,J)=L0
+                  end if
 c                  write(*,*)"o"
 c                  VARG(I,J)=IWORK(1)*1.0D0
 c                  VARG(I,J)=lambd0