diff --git a/R/hdist.R b/R/hdist.R
index b0d8bbfeee2517888fc5dc79baa1d43570998eb6..c355fd24977499e09897ca59e0c3eafabe2315d2 100644
--- a/R/hdist.R
+++ b/R/hdist.R
@@ -38,7 +38,7 @@
 #' @examples
 #' data(Oudon)
 #' catchments <- st_geometry(Oudon$obs)
-#' hdist(x = catchments[1:2], y = catchments[3:5], gres = 5, method = "rghosh")
+#' hdist(x = catchments[1:2], y = catchments[3:5], gres = 5, method = "rghosh", parallel=c(FALSE, TRUE), cores=2)
 #' @import sf stars doParallel foreach
 #' @importFrom units set_units drop_units
 #' @useDynLib transfR, .registration = TRUE
@@ -97,12 +97,12 @@ hdist.sfc <- function(x, y, method="rghosh", gres=5, ditself=FALSE, maxsample=2.
       }
     }
     if(verbose) cat("Computing Ghosh distance between catchments\n")
-    gdist <- loop_gdist(xdisc=xdisc, ydisc=ydisc, proj=proj, intersect=FALSE, parallel=c(FALSE, parallel[2]), cores=cores)
+    gdist <- loop_gdist(xdisc=xdisc, ydisc=ydisc, proj=proj, intersect=FALSE, parallel=parallel[2], cores=cores)
     if(ditself){
       if(verbose) cat("Computing Ghosh distance within catchments\n")
       if(!identical(x,y)){
-        gdist <- cbind(gdist,sapply(xdisc,FUN=function(x){call_gdist(pts1=x, pts2=x, proj=proj, rescale=FALSE, diag=TRUE, cores=cores, intersect=FALSE)}))
-        gdist <- rbind(gdist,c(sapply(ydisc,FUN=function(x){call_gdist(pts1=x, pts2=x, proj=proj, rescale=FALSE, diag=TRUE, cores=cores, intersect=FALSE)}),NA))
+        gdist <- cbind(gdist,sapply(xdisc,FUN=function(x){call_gdist(pts1=x, pts2=x, proj=proj, intersect=FALSE, rescale=FALSE, diag=TRUE, parallel=parallel[2], cores=cores)}))
+        gdist <- rbind(gdist,c(sapply(ydisc,FUN=function(x){call_gdist(pts1=x, pts2=x, proj=proj, intersect=FALSE, rescale=FALSE, diag=TRUE, parallel=parallel[2], cores=cores)}),NA))
       }else{
         gdist <- cbind(gdist,diag(gdist))
         gdist <- rbind(gdist,c(diag(gdist),NA))
@@ -112,7 +112,7 @@ hdist.sfc <- function(x, y, method="rghosh", gres=5, ditself=FALSE, maxsample=2.
     if(method=="rghosh2"){
       gc() # Clean memory
       if(verbose) cat("Computing Ghosh distance within the shared areas\n")
-      idist <- loop_gdist(xdisc=xdisc, ydisc=ydisc, proj=proj, intersect=TRUE, parallel=c(FALSE, parallel[2]), cores=cores)
+      idist <- loop_gdist(xdisc=xdisc, ydisc=ydisc, proj=proj, intersect=TRUE, parallel=parallel[2], cores=cores)
       if(verbose) cat("Rescaling Ghosh distance\n")
       gdist <- gdist-idist
     }
@@ -185,13 +185,14 @@ hdist.transfR <- function(x, y, method="rghosh", weightO=0.8, weightC=0.2, ...){
   return(mdist)
 }
 
-call_gdist <- function(pts1, pts2, proj, rescale, diag, cores, intersect){
+call_gdist <- function(pts1, pts2, proj, intersect, rescale, diag, parallel, cores){
+  if(!parallel) cores = 1
   if(intersect){
     # Compute the Ghosh distance matrix for common points only
     iseq <- st_equals(pts1, pts2, sparse = FALSE)
     pts_inter <- pts1[apply(iseq, MARGIN = 1, FUN = any)]
     if(length(pts_inter)>0){
-      call_gdist(pts1=pts_inter, pts2=pts_inter, proj=proj, rescale=FALSE, diag=FALSE, cores=cores, intersect=FALSE)
+      call_gdist(pts1=pts_inter, pts2=pts_inter, proj=proj, intersect=FALSE, rescale=FALSE, diag=FALSE, parallel=parallel, cores=cores)
     }else{
       return(0)
     }
@@ -203,28 +204,12 @@ call_gdist <- function(pts1, pts2, proj, rescale, diag, cores, intersect){
 }
 
 loop_gdist <- function(xdisc, ydisc, proj, intersect, parallel, cores){
-  if(length(parallel)==1) parallel[1:2] <- parallel
-  if(parallel[1]){ # Never used by the package, FORTRAN/openMP parallelization is preferred
-    if(!identical(xdisc,ydisc)){
-      res <- foreach::"%dopar%"(foreach::foreach(i=1:length(xdisc),.packages="sf"),
-                                sapply(ydisc,FUN=function(x){call_gdist(pts1=x, pts2=xdisc[[i]], proj=proj, rescale=FALSE, diag=FALSE, cores=1, intersect=intersect)}))
-      gdist <- matrix(unlist(res), nrow = length(xdisc), ncol = length(ydisc), byrow = TRUE)
-    }else{
-      res <- foreach::"%dopar%"(foreach::foreach(i=1:length(xdisc),.packages="sf"),{
-        tmp <- sapply(ydisc[i:length(xdisc)],FUN=function(x){call_gdist(pts1=x, pts2=xdisc[[i]], proj=proj, rescale=FALSE, diag=FALSE, cores=1, intersect=intersect)})
-        c(rep(NA,i-1),tmp)
-      })
-      gdist <- matrix(unlist(res), nrow = length(xdisc), ncol = length(ydisc), byrow = TRUE)
-      for(i in 1:length(xdisc)) gdist[,i] <- gdist[i,]
-    }
+  if(!parallel) cores = 1
+  gdist <- matrix(NA,nrow=length(xdisc),ncol=length(ydisc))
+  if(!identical(xdisc,ydisc)){
+    for(i in 1:length(xdisc)) gdist[i,] <- sapply(ydisc,FUN=function(x){call_gdist(pts1=x, pts2=xdisc[[i]], proj=proj, intersect=intersect, rescale=FALSE, diag=FALSE, parallel=parallel, cores=cores)})
   }else{
-    if(!parallel[2]) cores = 1
-    gdist <- matrix(NA,nrow=length(xdisc),ncol=length(ydisc))
-    if(!identical(xdisc,ydisc)){
-      for(i in 1:length(xdisc)) gdist[i,] <- sapply(ydisc,FUN=function(x){call_gdist(pts1=x, pts2=xdisc[[i]], proj=proj, rescale=FALSE, diag=FALSE, cores=cores, intersect=intersect)})
-    }else{
-      for(i in 1:length(xdisc)) gdist[i,i:length(xdisc)] <- gdist[i:length(xdisc),i] <- sapply(ydisc[i:length(xdisc)],FUN=function(x){call_gdist(pts1=x, pts2=xdisc[[i]], proj=proj, rescale=FALSE, diag=FALSE, cores=cores, intersect=intersect)})
-    }
+    for(i in 1:length(xdisc)) gdist[i,i:length(xdisc)] <- gdist[i:length(xdisc),i] <- sapply(ydisc[i:length(xdisc)],FUN=function(x){call_gdist(pts1=x, pts2=xdisc[[i]], proj=proj, intersect=intersect, rescale=FALSE, diag=FALSE, parallel=parallel, cores=cores)})
   }
   units(gdist) <- units(st_distance(xdisc[[1]][1],ydisc[[1]][1])) #units are lost by some operations
   return(gdist)
diff --git a/man/hdist.Rd b/man/hdist.Rd
index 69f384a5c3720c4ce8b22b4fc95e07a650e2d593..19c208b0b6d0fb1fd64c742408c2a9375d7fa611 100644
--- a/man/hdist.Rd
+++ b/man/hdist.Rd
@@ -51,9 +51,9 @@ It will add one row and one column in the distance matrix. Only used if method i
 \item{proj}{logical indicating if spatial layer are using a projection. If TRUE, euclidean
 distance is used. If FALSE, the great-circle distance is used}
 
-\item{parallel}{logical indicating if the computation should be parallelised. Could be a vector
+\item{parallel}{logical indicating whether the computation should be parallelised. Could be a vector
 of length two to distinguish between parallelisation of sampling and Ghosh distance
-(as sampling can be memory intensive)}
+(because sampling over large areas can be memory intensive)}
 
 \item{cores}{the number of cores to use for parallel execution if \code{parallel} is TRUE.
 If not specified, the number of cores is set to the value of \code{parallel::detectCores()}}