diff --git a/NEWS b/NEWS
index 0f184b6b73ae82e93e61d9a27714a5fb5ebe5e6f..ef890a572f3b260109ca5a908b3e03a785f3a5c5 100755
--- a/NEWS
+++ b/NEWS
@@ -5,7 +5,8 @@ NEWS/ChangeLog for hydroPSO
         o 'hydroPSO'    : -) added '...' parameter. It is only used when 'fn' is different from "hydromod". This is only done for 'optim' compatibility.
                           -) fixed small bug related to improper PSO evolution when 'best.update'=="async" (NOT the default option !) 
                              AND (topology!="gbest" | method!="pso")
-                          -) default value of 'IW' was changed from (linearly decreasing) IW.w= c(1.2, 0.4) to 'IW.w=1/(2*log(2))'
+                          -) new 'random.update' parameter for the 'control' variable, in order to allow random update of personal/global best
+                          -) default value of the inertia weight 'IW' was changed from (linearly decreasing) IW.w= c(1.2, 0.4) to 'IW.w=1/(2*log(2))'
                           -) source code was tidy up
                              
 0.1-56	14-Jun-2012
diff --git a/R/PSO_v2012.R b/R/PSO_v2012.R
index 7011b4537c7ff6575981bfa91cfc97a5708e5866..94182f031ae602ab9fb23a9abeb0b24e2e55e6e9 100755
--- a/R/PSO_v2012.R
+++ b/R/PSO_v2012.R
@@ -1387,6 +1387,7 @@ hydroPSO <- function(
 	    Xini.type=c("lhs", "random"),  
 	    Vini.type=c("lhs", "random", "zero"), 
 	    best.update=c("sync", "async"),
+	    random.update= TRUE,
 	    boundary.wall=c("reflecting", "damping", "absorbing", "invisible"),
 	    topology=c("random", "gbest", "lbest", "vonNeumann"), K=3, iter.ini=0, ngbest=4, # only used when 'method=ipso'   
 
@@ -1430,6 +1431,7 @@ hydroPSO <- function(
     Xini.type         <- match.arg(control[["Xini.type"]], con[["Xini.type"]]) 
     Vini.type         <- match.arg(control[["Vini.type"]], con[["Vini.type"]]) 
     best.update       <- match.arg(control[["best.update"]], con[["best.update"]]) 
+    random.update     <- as.logical(con[["rand.update"]])
     boundary.wall     <- match.arg(control[["boundary.wall"]], con[["boundary.wall"]]) 
     topology          <- match.arg(control[["topology"]], con[["topology"]]) 
     K                 <- con[["K"]]      
@@ -1792,11 +1794,11 @@ hydroPSO <- function(
       writeLines(c("Max Iterations    :", maxit), PSOparam.TextFile, sep=" ") 
       writeLines("", PSOparam.TextFile) 
       writeLines(c("Method            :", method), PSOparam.TextFile, sep=" ") 
+      writeLines("", PSOparam.TextFile) 
       if ( method == "ipso" ) {
 	writeLines(c("ngbest           :", ngbest), PSOparam.TextFile, sep=" ") 
 	writeLines("", PSOparam.TextFile)  
       } # IF end
-      writeLines("", PSOparam.TextFile) 
       writeLines(c("Topology          :", topology), PSOparam.TextFile, sep=" ") 
       writeLines("", PSOparam.TextFile) 
       if ( (topology == "lbest") | (topology == "random") ) {
@@ -2134,23 +2136,6 @@ hydroPSO <- function(
 
 	} # ELSE end
 
-      if ( plot ) {
-	ifelse(MinMax == "max", lgof <- max(GoF, na.rm=TRUE), lgof <- min(GoF, na.rm=TRUE)) 
-	colorRamp= colorRampPalette(c("darkred", "red", "orange", "yellow", "green", "darkgreen", "cyan"))
-	XX.Boundaries.current <- computeCurrentXmaxMin(X) 
-	xlim <- range(XX.Boundaries.current)
-	ylim <- range(XX.Boundaries.current)
-	if (iter==1) {
-	   plot(X[,1], X[,2], xlim=X.Boundaries[1,], ylim=X.Boundaries[2,], 
-	        main=paste("Iter= ", iter, ". GoF= ", 
-	        format(lgof, scientific=TRUE, digits=digits), sep=""), 
-	        col=colorRamp(npart), cex=0.5 )
-	} else plot(X[,1], X[,2], xlim=X.Boundaries[1,], ylim=X.Boundaries[2,], 
-	            main=paste("Iter= ", iter, ". GoF= ", 
-	            format(lgof, scientific=TRUE, digits=digits), sep=""), 
-	            col=colorRamp(npart), cex=0.5 )
-	#plotParticles2D(X)
-      } # IF end 
 
       if ( best.update == "sync" ) {
 	    tmp <- sync.update.pgbests(x=X, 
@@ -2192,7 +2177,10 @@ hydroPSO <- function(
       ##########################################################################  
       ###################   Particles Loop (j) - Start  ########################
       ##########################################################################  
-      for (j in 1:npart) {
+      
+      ifelse(random.update, index <- sample(npart), index <- 1:npart)
+        
+      for (j in index) {
       
         if (write2disk) {
         
@@ -2326,7 +2314,25 @@ hydroPSO <- function(
       } # FOR j end: Particles Loop
       ##########################################################################  
       ###################   Particles Loop (j) - End  ##########################
-      ##########################################################################  
+      ########################################################################## 
+       
+      if ( plot ) {
+	ifelse(MinMax == "max", lgof <- max(GoF, na.rm=TRUE), lgof <- min(GoF, na.rm=TRUE)) 
+	colorRamp= colorRampPalette(c("darkred", "red", "orange", "yellow", "green", "darkgreen", "cyan"))
+	XX.Boundaries.current <- computeCurrentXmaxMin(X) 
+	xlim <- range(XX.Boundaries.current)
+	ylim <- range(XX.Boundaries.current)
+	if (iter==1) {
+	   plot(X[,1], X[,2], xlim=X.Boundaries[1,], ylim=X.Boundaries[2,], 
+	        main=paste("Iter= ", iter, ". GoF= ", 
+	        format(lgof, scientific=TRUE, digits=digits), sep=""), 
+	        col=colorRamp(npart), cex=0.5 )
+	} else plot(X[,1], X[,2], xlim=X.Boundaries[1,], ylim=X.Boundaries[2,], 
+	            main=paste("Iter= ", iter, ". GoF= ", 
+	            format(lgof, scientific=TRUE, digits=digits), sep=""), 
+	            col=colorRamp(npart), cex=0.5 )
+	#plotParticles2D(X)
+      } # IF end 
       
       gbest.fit.iter[iter] <- gbest.fit
       
diff --git a/man/hydroPSO.Rd b/man/hydroPSO.Rd
index 5ccd7ccc689002d84424aef482e3d067992aa474..a739a1e416faf56ddb4ea03a04d8b36ad1352892 100755
--- a/man/hydroPSO.Rd
+++ b/man/hydroPSO.Rd
@@ -130,9 +130,13 @@ character, indicates how to initialise the particles' velocities in the swarm. V
 By default \code{Vini.type=lhs}
 }
   \item{best.update}{
-character, indicates how (when) to update the global and local best. Valid values are: \cr
+character, indicates how (when) to update the global/neighbourhood and personal best. Valid values are: \cr
 -)\kbd{sync}: the update is made synchronously, i.e. after computing the position and goodness-of-fit for ALL the particles in the swarm. This is the DEFAULT option\cr
 -)\kbd{async}: the update is made asynchronously, i.e. after computing the position and goodness-of-fit for EACH individual particle in the swarm
+}
+  \item{random.update}{
+OPTIONAL. Only used when \code{best.update='async'}. \cr
+logical, if \code{TRUE} the particles are processed in random order to update their personal best and the global/neighbourhood best. By default \code{random.update=TRUE}
 }
   \item{boundary.wall}{
 character, indicates the type of boundary condition to be applied during optimisation. Valid values are in \code{c('reflecting', 'damping', 'absorbing', 'invisible')} \cr