diff --git a/DESCRIPTION b/DESCRIPTION
index 1b8fadd2ce80a939636e20eff34188fcce5d4ba1..5fccc0ab3185dce2361eac6d0bb4cbce5caafd95 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: gmGeostats
-Version: 0.10-8-9001
-Date: 2021-07-15
+Version: 0.10.9-9001
+Date: 2021-10-01
 Title: Geostatistics for Compositional Analysis
 Authors@R: c(person(given = "Raimon", 
         family = "Tolosana-Delgado", 
diff --git a/NAMESPACE b/NAMESPACE
index feecde33e161434c24f0320222223cfc181784fc..1d5250eb6a6ad0f087a249ad96533cfa0f8e699f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -102,6 +102,7 @@ S3method(precision,accuracy)
 S3method(predict,LMCAnisCompo)
 S3method(predict,genDiag)
 S3method(predict,gmCgram)
+S3method(predict,gmSpatialModel)
 S3method(print,mask)
 S3method(setMask,DataFrameStack)
 S3method(setMask,GridTopology)
@@ -236,6 +237,7 @@ exportClasses(gmValidationStrategy)
 exportMethods(Predict)
 exportMethods(as.gstat)
 exportMethods(logratioVariogram)
+exportMethods(predict)
 exportMethods(stackDim)
 exportMethods(variogram)
 import(RColorBrewer)
diff --git a/NEWS.md b/NEWS.md
index 590be6e55bed0882f987d037abd4c13196cf66bc..446da0c0eb1559d9c38b2934fe394fe253fafcfe 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,13 @@
+# gmGeostats 0.10.9-9001
+
+* (2021-09-30) minor bug corrected on the way fit_lmc passes its arguments to surrogate gstat and gmGeostats functions
+* (2021-10-01) management of `predict()` methods for S3 and S4 objects, wth the help of the internal function `Predict()` (not recommended to use)
+
+# gmGeostats 0.10.9
+
+(release of polished 0.10.8.900x dev versions)
+* minor correction on C source code for compatibility with clang13
+
 # gmGeostats 0.10.8-9001
 
 * (2021-07-22) conversion methods between variogram models: `as.gmCgram()` methods for "variogramModel" and "variogramModelList" objects (from package "gstat")
diff --git a/R/gmSpatialModel.R b/R/gmSpatialModel.R
index 89bd9e1e1f5b685dc30b6c2508eb6892f1d81bb7..27b5bed1cc436c4358819350d0aa088b11386293 100644
--- a/R/gmSpatialModel.R
+++ b/R/gmSpatialModel.R
@@ -283,6 +283,8 @@ make.gmCompositionalMPSSpatialModel = function(
 
 ### exporter to gstat -----
 as.gstat.gmSpatialModel <- function(object, ...){
+  # extra arguments
+  lldots = list(...)
   # data elements
   coords = sp::coordinates(object)
   X = compositions::rmult(object@data, V= gsi.getV(object@data), orig=gsi.orig(object@data))
@@ -306,10 +308,17 @@ as.gstat.gmSpatialModel <- function(object, ...){
     if(any(is.infinite(beta))) beta = NULL
     # neighbourhood
     ng = object@parameters
+    # manage manual changes of parameters given in dots...
+    for(nm in names(ng)){
+      tk = grep(nm, names(lldots))
+      if(length(tk)>1) ng[[tk]] = lldots[[tk]] 
+    }
+    # convert
     if(!is(ng, "gmKrigingNeighbourhood")) stop("as.gstat: object@parameters must be of class 'gmGaussianMethodParameters'!")
     res = compo2gstatLR(coords=coords, compo=compo, V=Vinv, lrvgLMC=lrvgLMC, 
                         nscore=FALSE, formulaterm = formulaterm, prefix=prefix, beta=beta, 
-                        nmax=ng$nmax, nmin=ng$nmin, omax=ng$omax, maxdist=ng$maxdist, force=ng$force)
+                        nmax=ng$nmax, nmin=ng$nmin, omax=ng$omax, maxdist=ng$maxdist, force=ng$force, 
+                        ...)
     return(res)
   }else{
     # non-compositional case
@@ -320,9 +329,15 @@ as.gstat.gmSpatialModel <- function(object, ...){
     if(any(is.infinite(beta))) beta = NULL
     # neighbourhood
     ng = object@parameters
+    # manage manual changes of parameters given in dots...
+    for(nm in names(ng)){
+      tk = grep(nm, names(lldots))
+      if(length(tk)>1) ng[[tk]] = lldots[[tk]] 
+    }
     res = rmult2gstat(coords=coords, data=X, V=V, vgLMC=vgLMC, 
                     nscore=FALSE, formulaterm = formulaterm, prefix=prefix, beta=beta, 
-                    nmax=ng$nmax, nmin=ng$nmin, omax=ng$omax, maxdist=ng$maxdist, force=ng$force)    
+                    nmax=ng$nmax, nmin=ng$nmin, omax=ng$omax, maxdist=ng$maxdist, force=ng$force,
+                    ...)    
   }
 }
 
@@ -400,26 +415,38 @@ as.gmSpatialModel.gstat = function(object, V=NULL, ...){
 #' @name predict_gmSpatialModel
 NULL
 
+
+
 #' @rdname predict_gmSpatialModel
-setMethod("predict", signature(object="gmSpatialModel"),
-          function(object, newdata, pars, ...){
-            Predict(object, newdata, pars, ...)
-          })
-          
+#' @export
+#' @method predict gmSpatialModel
+predict.gmSpatialModel <- function(object, newdata, pars=NULL, ...){
+  if(is.null(pars)){
+    return(Predict(object, newdata, ...))
+  }else{
+    return(Predict(object, newdata, pars, ...))
+  }
+}
+
+#' @rdname predict_gmSpatialModel
+#' @export
+setMethod("predict", signature(object="gmSpatialModel"), definition = predict.gmSpatialModel)
+
 
 
 #' @rdname predict_gmSpatialModel
 #' @include gmAnisotropy.R
 #' @include preparations.R
-#'  @export
-setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="missing"),
-          function(object, newdata, pars, ...){
-            if(is.null(object$pars)) object$pars = KrigingNeighbourhood()
-            predict(object, newdata, pars = object$pars, ...)
+#' @export
+setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY"),
+          function(object, newdata, ...){
+            if(is.null(object@parameters)) object@parameters = KrigingNeighbourhood() 
+            Predict(object, newdata, pars = object@parameters , ...)
           }
 )
 
 
+
 #' @rdname predict_gmSpatialModel
 #' @include gmSpatialMethodParameters.R
 #' @export
diff --git a/R/gstatCompatibility.R b/R/gstatCompatibility.R
index c7bfed67a37734311d10fd50dbeafbec713bc9ed..8f141edee7f75ddf23824464a1195cbc429626a9 100644
--- a/R/gstatCompatibility.R
+++ b/R/gstatCompatibility.R
@@ -30,10 +30,15 @@ fit_lmc  <- function(v, ...) UseMethod("fit_lmc", v)
 #' @describeIn fit_lmc wrapper around gstat::fit.lmc method
 #' @param g spatial data object, containing the original data
 #' @param model LMC or variogram model to fit
+#' @param fit.ranges logical, should ranges be modified? (default=FALSE)
+#' @param fit.lmc logical, should the nugget and partial sill matrices be modified (default=TRUE) 
+#' @param correct.diagonal positive value slightly larger than 1, for multiplying the direct variogram 
+#' models and reduce the risk of numerically negative eigenvalues
 #' @export
 #' @method fit_lmc gstatVariogram
-fit_lmc.gstatVariogram <- function(v, g, model,...){
-  res = gstat::fit.lmc(as.gstatVariogram(v, ...), as.gstat(g, ...), as.variogramModel(model, ...))
+fit_lmc.gstatVariogram <- function(v, g, model, fit.ranges = FALSE, fit.lmc = !fit.ranges, correct.diagonal = 1.0, ...){
+  res = gstat::fit.lmc(as.gstatVariogram(v, ...), as.gstat(g, ...), as.variogramModel(model, ...),
+                       fit.ranges = fit.ranges, fit.lmc = fit.lmc, correct.diagonal=correct.diagonal)
   class(res$model) = c("variogramModelList", "list")
   return(res)
 }
diff --git a/src/gmGeostats.c b/src/gmGeostats.c
index 4ad58861cc36daabc25dd2e0b7536e68369be246..c2c4cf3acf4957ca93415fe75c1d22506e4fc397 100644
--- a/src/gmGeostats.c
+++ b/src/gmGeostats.c
@@ -6,10 +6,10 @@
  */
 // attention: comment this if not compiling
 #include <stdio.h>
-#include <Rinternals.h>
 #ifdef _OPENMP
 #include <omp.h>
 #endif
+#include <Rinternals.h>
 
 #define inR   // attention: this must be uncommented if not compiling
 
diff --git a/src/gmGeostats.o b/src/gmGeostats.o
index e09d652bcedc29a6557f12fb8edec9a882e43e69..c55dc18ec61857d88b19575da5ab7536fc49a5ed 100644
Binary files a/src/gmGeostats.o and b/src/gmGeostats.o differ
diff --git a/src/gmGeostats.so b/src/gmGeostats.so
index 4403983c462507fa93aaf58f0b1f1aa083135264..68c8127550c28b8d5a6a138659a5eecfbfb561cf 100755
Binary files a/src/gmGeostats.so and b/src/gmGeostats.so differ