diff --git a/R/PSO_v2012.R b/R/PSO_v2012.R
index 0948dfca38fe2d29ce19f69e43e5499f911c3bbb..a6023b40216fe932c5d3f7d8783f1827b2846cab 100755
--- a/R/PSO_v2012.R
+++ b/R/PSO_v2012.R
@@ -2002,7 +2002,6 @@ hydroPSO <- function(
 
     if (topology != "random") {
       nc <- K  
-      #nc <- npart
       if (trunc(K/2) != ceiling(K/2)) {
         N   <- (K-1)/2
       } else N  <- K/2
@@ -2020,8 +2019,6 @@ hydroPSO <- function(
 	  X.neighbours[i,j+N+NN] <- neigh.index
 	} # FOR end
       } # FOR end                      
-      #rownames(X.neighbours) <- paste("Part", 1:npart, sep="")
-      #colnames(X.neighbours) <- paste("Neigh", 1:nc, sep="")
     } # IF end 
 
     LocalBest.fit <- rep(fn.worst.value, npart)
@@ -2655,7 +2652,7 @@ hydroPSO <- function(
       NormSwarmRadius <- swarm.radius/swarm.diameter
 
       if ( (verbose) & ( iter/REPORT == floor(iter/REPORT) ) ) 
-	   message( "iter:", format(iter, width=6, justify="right"), 
+	   message( "iter:", format(iter, width=nchar(maxit), justify="right"), 
 		    "  Gbest:", formatC( gbest.fit, format="E", digits=digits, flag=" "), 
 		    "  Gbest_rate:", format( round(gbest.fit.rate*100, 2), width=6, nsmall=2, justify="left"), "%",
 		    "  Iter_best_fit:", formatC(pbest.fit.iter, format="E", digits=digits, flag=" "),               
@@ -2917,8 +2914,6 @@ hydroPSO <- function(
 
       # Writing the file 'X.neighbours.txt' 
       fname <- paste(file.path(drty.out), "/", "Particles_Neighbours.txt", sep="") 
-      print(X.neighbours)	
-      print(paste("Neigh", 1:npart, sep=""))
       ifelse(topology == "lbest", nc <- K, nc <- npart)
       write.table(X.neighbours, file=fname, col.names=paste("Neigh", 1:nc, sep=""), row.names=paste("Part", 1:npart, sep=""), sep="  ", na="", quote=FALSE) 
 
diff --git a/R/plot_2parOF.R b/R/plot_2parOF.R
index fd20afa683c21ae5b411144a8617ac84b657f19e..31f58875127cc0efd4a5c3637c8cdb636c86b4c4 100755
--- a/R/plot_2parOF.R
+++ b/R/plot_2parOF.R
@@ -89,11 +89,6 @@ plot_2parOF <- function(params,
     ifelse(MinMax=="max", decreasing<-FALSE, decreasing<-TRUE)
     p <- p[order(p[, gof.name], decreasing = decreasing), ]
     
-    # Define a transformation
-    norm_trans <- function(){
-      trans_new('norm', function(x) pnorm(x), function(x) qnorm(x))
-    }
-  
     if (type=="sp") {
       colnames(p) <- c("x", "y", gof.name)
       
diff --git a/R/plot_NparOF.R b/R/plot_NparOF.R
index 4bfa5714b88adb35d7cdddfad7f720c33c466df6..7ab0c3b2c0b4448d6b4692c7280e483338a8d9a8 100755
--- a/R/plot_NparOF.R
+++ b/R/plot_NparOF.R
@@ -7,15 +7,15 @@
 ################################################################################
 #                             'plot_NparOF'                                    #
 ################################################################################
-# Author : Mauricio Zambrano Bigarini                                          #
+# Author : Mauricio Zambrano Bigiarini                                         #
 # Started: Nov 30th, 2010                                                      #   
 # Updates: 17-Jan-2011 ; 28-Jan-2011 ; 09-Mar-2011                             #
-#          17-Feb-2012 ; 21-Feb-2012 ; 09-Mar-2012 ; 23-Mar-2012               #    
+#          17-Feb-2012 ; 21-Feb-2012 ; 09-Mar-2012 ; 23-Mar-2012 ; 19-Nov-2012 #    
 ################################################################################
 # Purpose: For 'n' user-defined parameters, it produces 'sum(1:(npar-1))'      #
-#         'plot_2parOF' plots, with the  values of the objective funtion in    #
+#         'plot_2parOF' plots, with the  values of the objective function in   #
 #         a 2D box,  where the boundaries of each parameter are used as axis.  #
-#         The 'sum(1:(npar-1)) plots corresponds to all the posible            #
+#         The 'sum(1:(npar-1)) plots corresponds to all the possible           #
 #         combinations of 2 parameters among all the 'n' parameters provided   #
 ################################################################################
 # nrows  : numeric, with the amount of rows to be used in the plotting window. 
@@ -26,6 +26,7 @@ plot_NparOF <- function(params,
                         gofs,
                         param.names=colnames(params),
                         MinMax=c("min", "max"),
+                        beh.thr=NA, 
                         nrows="auto",
                         gof.name="GoF", 
                         main=paste(gof.name, "Surface"),
@@ -38,6 +39,10 @@ plot_NparOF <- function(params,
                         ) {
 
     
+  ##############################################################################
+  # 1)                            Checkings                                    #
+  ##############################################################################
+  
     # Checking 'params'
     if (missing(params)) 
       stop("Missing argument: 'params' must be provided !!" )
@@ -57,7 +62,7 @@ plot_NparOF <- function(params,
     # Number of parameters that will be analysed
     npar <- length(param.names)
 
-    # creating the varaible that will store the position of the selected parameters within 'params'
+    # creating the variable that will store the position of the selected parameters within 'params'
     par.pos <- numeric(npar)
 
     # Checking 'param.names'
@@ -66,7 +71,32 @@ plot_NparOF <- function(params,
         stop("Invalid argument: The field '", param.names[i], "' doesn't exist in 'params'")
       
       par.pos[i] <- which(colnames(params) == param.names[i])
-    } # FOR end  
+    } # FOR end
+  
+    
+  ##############################################################################
+  # 2)                            Computations                                 #
+  ##############################################################################
+  
+    # Filtering out those parameter sets above/below a certain threshold
+    if (!is.na(beh.thr)) {  
+       # Checking 'beh.thr'
+       mx <- max(gofs, na.rm=TRUE)
+       if (beh.thr > mx)
+         stop("Invalid argument: 'beh.thr' must be lower than ", mx ,"!!")
+    
+      # Computing the row index of the behavioural parameter sets
+      ifelse(MinMax=="min", beh.row.index <- which(gofs <= beh.thr), 
+                            beh.row.index <- which(gofs >= beh.thr) )
+    
+      # Removing non-behavioural parameter sets & gofs
+      params <- params[beh.row.index, ]
+      gofs   <- gofs[beh.row.index]
+   
+      # Amount of behavioural parameter sets 
+      nbeh <- nrow(params)
+      if (verbose) message( "[ Number of behavioural parameter sets: ", nbeh, " ]" )
+    } # IF end  
 
     # If the user didn't provide 'GOFcuts', the 5 quantiles are used
     if (length(GOFcuts) == 1){
@@ -74,6 +104,10 @@ plot_NparOF <- function(params,
         GOFcuts <- unique( quantile( as.numeric(gofs), 
                            probs=c(0, 0.25, 0.5, 0.75, 0.9, 0.95, 1), na.rm=TRUE) )
     } # IF end
+    
+  ##############################################################################
+  # 3)                            Plotting                                     #
+  ##############################################################################  
    
     # Number of plots that will be drawn   
     nplots <- sum(1:(npar-1))
diff --git a/R/plot_params.R b/R/plot_params.R
index f5cc5512cab0644c8981b767a300ee79ca15b531..592e56c720739df37ab0e1be79457e1e17ec8264 100755
--- a/R/plot_params.R
+++ b/R/plot_params.R
@@ -17,7 +17,7 @@
 # corresponding objective function value (usually for plotting the 
 # efficiencies of different parameter sets)
 
-# params : data.frame, whose colums represents the behavioural parameter 
+# params : data.frame, whose columns represents the behavioural parameter 
 #          sets and the goodness of fit of each one of them 
 # of.name: character, with the name of the column that contains the values 
 #          of the objective function for each parameter set
@@ -156,9 +156,9 @@ plot_params.default <- function(params,
   } # IF end
   
   ##############################################################################
-  # 2)                            Plotting                                    #
+  # 3)                            Plotting                                     #
   ##############################################################################  
-  # Ploting ALL the PARAMETER SETS
+  # Plotting ALL the PARAMETER SETS
   if (verbose) message( "                                        ")  
   if (verbose) message( "[            Plotting ...              ]")  
   
@@ -181,7 +181,7 @@ plot_params.default <- function(params,
   old.par <- par(no.readonly=TRUE)
   if (!do.png) on.exit(par(old.par)) 
 
-  # Plotting the distribution of all the parameter sampled with LH, not only the behaviourals
+  # Plotting the distribution of all the parameter sampled with LH, not only the behavioural ones
   MinMax.colour <- "coral"
   par(mfrow=c(lnr, ncol))
   par(mar=c(5,4.5,1,2)+0.1) # Default: par(mar=c(5,4,4,2)+0.1)
diff --git a/R/read_params.R b/R/read_params.R
index 379e6dac6e0ff0fedfc357b764e1d530ff5ce44c..6be6310ba9ccf54d1077f4da04b7de4bf1b2cfb0 100755
--- a/R/read_params.R
+++ b/R/read_params.R
@@ -31,7 +31,7 @@
 #          the file format: header is set to TRUE if and only if the first row 
 #          contains one fewer field than the number of columns.
 # param.cols: numerical vector, with the column position in 'param.fname' that 
-#             represent the paramter values used during the calibrations
+#             represent the parameter values used during the calibrations
 # param.names: character, with the names to be given to each one of the parameters
 #              stored in the columns represented by 'param.cols' (usually, the 
 #              most sensitive parameters in 'param.fname' )
@@ -41,7 +41,7 @@
 # plot   : logical, indicating if dotty plots of the parameter sets and its
 #          objective function have to be plotted or not
 # beh.thr: OPTIONAL, only used when 'plot=TRUE'. \cr
-#          numeric, with a behavioural threshold to be used for ploting 
+#          numeric, with a behavioural threshold to be used for plotting 
 #          a horizontal line 
 # beh.col: OPTIONAL, only used when 'plot=TRUE'. \cr
 #          color to be used for plotting the horizontal line on 'beh.thr'
diff --git a/man/ReadPlot_params.Rd b/man/ReadPlot_params.Rd
index 4330db4dba23ec80dc8670364c2f751e730d8c46..7856f7d981d89c8af12af79e70f71336c9a9ace1 100755
--- a/man/ReadPlot_params.Rd
+++ b/man/ReadPlot_params.Rd
@@ -123,7 +123,9 @@ Valid values are in: \code{c('min', 'max')}
 }
   \item{beh.thr}{
 OPTIONAL \cr
-numeric, value for drawing a horizontal line for separating behavioural from non behavioural parameter sets
+numeric, threshold value used for selecting parameter sets that have to be used in the analysis (\sQuote{behavioural parameters}, using the GLUE terminology) \cr
+If \code{MinMax='min'}, only parameter sets with a goodness-of-fit value (given by \code{gofs}) less than or equal to \code{beh.thr} will be considered for the subsequent analysis. \cr
+If \code{MinMax='max'}, only parameter sets with a goodness-of-fit value (given by \code{gofs}) greater than or equal to \code{beh.thr} will be considered for the subsequent analysis
 }
   \item{beh.col}{
 OPTIONAL. Only used when \code{plot=TRUE} \cr
diff --git a/man/plot_NparOF.Rd b/man/plot_NparOF.Rd
index 26ba684d9c212f2cfa925f7e8d0a8dadda464156..b2ff4747e0bbd4f66007a836cb6535b1a90a7514 100755
--- a/man/plot_NparOF.Rd
+++ b/man/plot_NparOF.Rd
@@ -16,10 +16,10 @@ The \kbd{sum(1:(npar-1))} plots corresponds to all the possible combinations of
 }
 \usage{
 plot_NparOF(params, gofs, param.names=colnames(params), MinMax=c("min", "max"), 
-           nrows="auto",  gof.name="GoF", main=paste(gof.name, "Surface"), 
-           GOFcuts="auto", colorRamp= colorRampPalette(c("darkred", "red", 
-           "orange", "yellow", "green", "darkgreen", "cyan")), points.cex=0.7, 
-           alpha=1, axis.rot=c(0, 0), verbose=TRUE)
+         beh.thr=NA, nrows="auto", gof.name="GoF", main=paste(gof.name, "Surface"), 
+         GOFcuts="auto", colorRamp= colorRampPalette(c("darkred", "red", 
+         "orange", "yellow", "green", "darkgreen", "cyan")), points.cex=0.7, 
+         alpha=1, axis.rot=c(0, 0), verbose=TRUE)
 }
 %- maybe also 'usage' for other objects documented here.
 \arguments{
@@ -34,6 +34,12 @@ character, names for the parameters in \code{params} that have to be plotted (\c
 }
   \item{MinMax}{
 character, indicates whether the optimum value in \code{gofs} corresponds to the minimum or maximum of the objective function. Valid values are in: \code{c('min', 'max')}
+}
+  \item{beh.thr}{
+OPTIONAL \cr
+numeric, threshold value used for selecting parameter sets that have to be used in the analysis (\sQuote{behavioural parameters}, using the GLUE terminology) \cr
+If \code{MinMax='min'}, only parameter sets with a goodness-of-fit value (given by \code{gofs}) less than or equal to \code{beh.thr} will be considered for the subsequent analysis. \cr
+If \code{MinMax='max'}, only parameter sets with a goodness-of-fit value (given by \code{gofs}) greater than or equal to \code{beh.thr} will be considered for the subsequent analysis
 }
   \item{nrows}{
 numeric, number of rows to be used in the plotting window \cr
@@ -88,7 +94,7 @@ Mauricio Zambrano-Bigiarini, \email{mzb.devel@gmail.com}
 %% ~Make other sections like Warning with \section{Warning }{....} ~
 
 \seealso{
-\code{\link{plot_2parOF}}, \code{\link{read_results}}, \code{\link{plot_results}}, \code{\link{plot_GofPerParticle}}, \code{\link{plot_ParamsPerIter}}
+\code{\link{plot_2parOF}}, \code{\link{read_results}}, \code{\link{plot_results}}, \code{\link{plot_GofPerParticle}}, , \code{\link{plot_params}}, \code{\link{plot_ParamsPerIter}}
 }
 \examples{
 # Number of dimensions to be optimised
@@ -103,9 +109,9 @@ setwd("~")
 set.seed(100)
 
 # Running PSO with the 'rosenbrock' test function, writing the results to text files
-hydroPSO(fn="rosenbrock",
-        lower=rep(-30, nparam), upper=rep(30, nparam),    
-        control=list(MinMax="min", npart=2*nparam, write2disk=TRUE)
+hydroPSO(fn=rosenbrock,
+        lower=rep(-10, nparam), upper=rep(10, nparam),    
+        control=list(write2disk=TRUE)
         ) # hydroPSO
   
 # reading the 'Particles.txt' output file of hydroPSO