diff --git a/R/Analogs.R b/R/Analogs.R
index a69a66dfc80b0943bffcfc3085da5d25145edcff..f0ebd610721f0a7ddcabfbec495436bb4d12ef26 100644
--- a/R/Analogs.R
+++ b/R/Analogs.R
@@ -66,8 +66,8 @@
 #'@param loocv_window a logical vector only to be used if 'obs' does not have the dimension
 #''window'. It indicates whether to apply leave-one-out cross-validation in the creation 
 #'of the window. It is recommended to be set to TRUE. Default to TRUE.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
-#'
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@return An 's2dv' object. The element 'data' contains the dowscaled field, 'lat' the
 #'downscaled latitudes, and 'lon' the downscaled longitudes. If fun_analog is set to NULL
 #'(default), the output array in 'data' also contains the dimension 'analog' with the best 
@@ -87,7 +87,7 @@
 #'@export
 CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analog = NULL, lat_dim = "lat", 
                         lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", member_dim = "member",
-                        region = NULL, return_indices = FALSE, loocv_window = TRUE, ncores = 1) {
+                        region = NULL, return_indices = FALSE, loocv_window = TRUE, ncores = NULL) {
 
   # input exp and obs must be s2dv_cube objects
   if (!inherits(exp,'s2dv_cube')) {
@@ -204,7 +204,8 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo
 #'@param loocv_window a logical vector only to be used if 'obs' does not have the dimension
 #''window'. It indicates whether to apply leave-one-out cross-validation in the creation 
 #'of the window. It is recommended to be set to TRUE. Default to TRUE.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@import multiApply
 #'@import CSTools
 #'@importFrom s2dv InsertDim CDORemap
@@ -232,7 +233,7 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp,
                     obs2_lats = NULL, obs2_lons = NULL, nanalogs = 3, fun_analog = NULL, 
                     lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", 
                     member_dim = "member", region = NULL, return_indices = FALSE, 
-                    loocv_window = TRUE, ncores = 1) {
+                    loocv_window = TRUE, ncores = NULL) {
   #-----------------------------------
   # Checkings
   #-----------------------------------
@@ -336,7 +337,14 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp,
       stop("Missing time dimension in 'obs2', or does not match the parameter 'time_dim'")
     }
   }
-
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
   # Select a function to apply to the analogs selected for a given observation
   if (!is.null(fun_analog)) {
     stopifnot(fun_analog %in% c("mean", "wmean", "max", "min", "median"))
@@ -370,22 +378,22 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp,
 
   obs_interpolated <- Interpolation(exp = obs_train, lats = obs_train_lats, lons = obs_train_lons, 
                                     target_grid = grid_exp, lat_dim = lat_dim, lon_dim = lon_dim, 
-                                    method_remap = "conservative", region = region)
+                                    method_remap = "conservative", region = region, ncores = ncores)
   # If after interpolating 'obs' data the coordinates do not match, the exp data is interpolated to
   # the same grid to force the matching
   if (!.check_coords(lat1 = obs_interpolated$lat, lat2 = exp_lats, lon1 = obs_interpolated$lon, lon2 = exp_lons)) {
     exp_interpolated <- Interpolation(exp = exp, lats = exp_lats, lons = exp_lons, target_grid = grid_exp, 
                                       lat_dim = lat_dim, lon_dim = lon_dim, method_remap = "conservative", 
-                                      region = region)$data
+                                      region = region, ncores = ncores)$data
   } else {
     exp_interpolated <- exp
   }
 
   # Create window if user does not have it in the training observations
   if ( !("window" %in% names(dim(obs_interpolated$data))) ) {
-    obs_train_interpolated <- generate_window(obj = obs_interpolated$data, sdate_dim = sdate_dim, 
-                                              time_dim = time_dim, loocv = loocv_window)
-    obs_hres <- generate_window(obj = obs, sdate_dim = sdate_dim, time_dim = time_dim, loocv = loocv_window)
+    obs_train_interpolated <- .generate_window(obj = obs_interpolated$data, sdate_dim = sdate_dim, 
+                                              time_dim = time_dim, loocv = loocv_window, ncores = ncores)
+    obs_hres <- .generate_window(obj = obs, sdate_dim = sdate_dim, time_dim = time_dim, loocv = loocv_window, ncores = ncores)
   }
 
   #-----------------------------------
@@ -411,7 +419,7 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp,
       ntimes <- dim(res.ind)[names(dim(res.ind)) == time_dim]
       res.ind <- Apply(res.ind, target_dims = c("ind", sdate_dim), function(x)   
         sapply(1:nsdates, function(s) seq(ntimes * nsdates)[ - (ntimes * (s - 1) + 1:ntimes)][x[, s]]),
-        output_dims = c("ind", sdate_dim))$output1
+        output_dims = c("ind", sdate_dim), ncores = ncores)$output1
     }
 
     # restore ensemble dimension in observations if it existed originally
@@ -484,7 +492,7 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp,
 # Add the dimension window to an array that contains, at least, the start date and time
 # dimensions
 # object has at least dimensions sdate and time
-generate_window <- function(obj, sdate_dim, time_dim, loocv, size = NULL) {
+.generate_window <- function(obj, sdate_dim, time_dim, loocv, size = NULL, ncores = NULL) {
 
   rsdates <- 1:dim(obj)[names(dim(obj)) == sdate_dim]
   ntimes <- dim(obj)[names(dim(obj)) == time_dim] 
@@ -497,11 +505,11 @@ generate_window <- function(obj, sdate_dim, time_dim, loocv, size = NULL) {
     if (loocv) {
         obj_window <- Apply(obj, target_dims = c(time_dim, sdate_dim),
           fun = function(x) sapply(rsdates, function(s) as.vector(x[ rtimes, -s])), 
-          output_dims = c('window', sdate_dim))$output1
+          output_dims = c('window', sdate_dim), ncores = ncores)$output1
     # Generate window without cross-validation        
     } else {
         obj_window <- Apply(obj, target_dims = c(time_dim, sdate_dim),
-          fun = as.vector, output_dims = 'window')$output1
+          fun = as.vector, output_dims = 'window', ncores = ncores)$output1
     }
   }
   # Generate window of the size specified by the user. Only applied with CV
@@ -513,28 +521,20 @@ generate_window <- function(obj, sdate_dim, time_dim, loocv, size = NULL) {
    
     # Concatenate data from previous, target and posterior months
     obj_new <- Apply(obj, target_dims = list(c("time", "smonth")),
-      fun = as.vector, output_dims = "time")$output1
+      fun = as.vector, output_dims = "time", ncores = ncores )$output1
 
     if (loocv) {
       obj_window <- Apply(list(obj_new, rtimes, rsdates), target_dims = list(c(time_dim, sdate_dim), NULL, NULL),
         fun = function(x, t, s) as.vector(x[(ntimes + t - size):(ntimes + t + size), -s]), 
-        output_dims = 'window')$output1  
+        output_dims = 'window', ncores = ncores)$output1  
       names(dim(obj_window))[(length(names(dim(obj_window))) - 1):length(names(dim(obj_window)))] <- c(time_dim, sdate_dim) 
     } else {
       obj_window <- Apply(obj_new, target_dims = c(time_dim, sdate_dim),
         fun = function(x) sapply(rtimes, function(t) as.vector(x[(ntimes + t - size):(ntimes + t + size), ])),
-        output_dims = c('window', time_dim))$output1
+        output_dims = c('window', time_dim), ncores = ncores)$output1
 
     }
   } 
 
   return(obj_window)
 }
-
-
-
-
-
-
-
-
diff --git a/R/Intbc.R b/R/Intbc.R
index 27de3735dbe88f0a60f8892406e287b55f40a0f7..1cb558d5050c4a5f29062820246e9fe374de80d4 100644
--- a/R/Intbc.R
+++ b/R/Intbc.R
@@ -49,7 +49,8 @@
 #'to the left border, while lonmax refers to the right border. latmin indicates the lower
 #'border, whereas latmax indicates the upper border. If set to NULL (default), the function
 #'takes the first and last elements of the latitudes and longitudes.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@return An 's2dv' object. The element 'data' contains the dowscaled field, 'lat' the
 #'downscaled latitudes, and 'lon' the downscaled longitudes.
 #'@examples
@@ -68,7 +69,7 @@
 
 CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, points = NULL, 
                       method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", 
-                      sdate_dim = "sdate", member_dim = "member", region = NULL, ncores = 1,...)
+                      sdate_dim = "sdate", member_dim = "member", region = NULL, ncores = NULL, ...)
 {
   if (!inherits(exp,'s2dv_cube')) {
     stop("Parameter 'exp' must be of the class 's2dv_cube'")
@@ -83,7 +84,7 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point
                int_method = int_method, bc_method = bc_method, points = points, 
                source_file = exp$attrs$source_files[1], method_point_interp = method_point_interp, 
                lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, member_dim = member_dim, 
-               region = region, ncores = ncores,...)
+               region = region, ncores = ncores, ...)
 
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   exp$data <- res$data
@@ -162,8 +163,8 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point
 #'to the left border, while lonmax refers to the right border. latmin indicates the lower
 #'border, whereas latmax indicates the upper border. If set to NULL (default), the function
 #'takes the first and last elements of the latitudes and longitudes.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
-#'
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@import CSTools 
 #'
 #'@seealso \code{\link[CSTools]{BiasCorrection}}
@@ -186,7 +187,7 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point
 #'@export
 Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid, bc_method, int_method = NULL,
                   points = NULL, method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", 
-                  time_dim = "time", member_dim = "member", source_file = NULL, region = NULL, ncores = 1, ...) {
+                  time_dim = "time", member_dim = "member", source_file = NULL, region = NULL, ncores = NULL, ...) {
 
   if (!inherits(bc_method, 'character')) {
     stop("Parameter 'bc_method' must be of the class 'character'")
@@ -250,11 +251,18 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid,
             "'obs_lats' and 'obs_lons'.")
     region <- c(obs_lons[1], obs_lons[length(obs_lons)], obs_lats[1], obs_lats[length(obs_lats)])
   }
-
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
   exp_interpolated <- Interpolation(exp = exp, lats = exp_lats, lons = exp_lons, target_grid = target_grid, 
                                     method_remap = int_method, points = points, source_file = source_file,
                                     lat_dim = lat_dim, lon_dim = lon_dim, method_point_interp = method_point_interp, 
-                                    region = region)
+                                    region = region, ncores = ncores)
  
   # If after interpolating 'exp' data the coordinates do not match, the obs data is interpolated to
   # the same grid to force the matching
@@ -263,7 +271,7 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid,
     obs_interpolated <- Interpolation(exp = obs, lats = obs_lats, lons = obs_lons, target_grid = target_grid, 
                                       method_remap = int_method, points = points, source_file = source_file,
                                       lat_dim = lat_dim, lon_dim = lon_dim, 
-                                      method_point_interp = method_point_interp, region = region)
+                                      method_point_interp = method_point_interp, region = region, ncores = ncores)
     obs_ref <- obs_interpolated$data
   } else {
     obs_ref <- obs
@@ -289,32 +297,32 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid,
     }
 
     res <- BiasCorrection(exp = exp_interpolated$data, obs = obs_ref, memb_dim = member_dim,
-                          sdate_dim = sdate_dim, ...) 
+                          sdate_dim = sdate_dim, ncores = ncores, ...) 
   } 
   else if (bc_method == 'cal' | bc_method == 'calibration') {
     if (dim(exp_interpolated$data)[member_dim] == 1) {
       stop('Calibration must not be used with only one ensemble member.')
     }
     res <- Calibration(exp = exp_interpolated$data, obs = obs_ref, memb_dim = member_dim, 
-                       sdate_dim = sdate_dim, ...)
+                       sdate_dim = sdate_dim, ncores = ncores, ...)
   }
   else if (bc_method == 'qm' | bc_method == 'quantile_mapping') {
    
     res <- QuantileMapping(exp = exp_interpolated$data, obs = obs_ref, na.rm = TRUE, 
-                           memb_dim = member_dim, sdate_dim = sdate_dim, ...)
+                           memb_dim = member_dim, sdate_dim = sdate_dim, ncores = ncores, ...)
   }
   else if (bc_method == 'dbc' | bc_method == 'dynamical_bias') {
     # the temporal dimension must be only one dimension called "time"
     if (all(c(time_dim, sdate_dim) %in% names(dim(exp_interpolated$data)))) {
       exp_interpolated$data <- Apply(exp_interpolated$data, target_dims = c(time_dim, sdate_dim), 
-        fun = as.vector, output_dims = "time")$output1
+        fun = as.vector, output_dims = "time", ncores = ncores)$output1
     } 
     if (all(c(time_dim, sdate_dim) %in% names(dim(obs_ref)))) {
       obs_ref <- Apply(obs_ref, target_dims = c(time_dim, sdate_dim), fun = as.vector, 
-        output_dims = "time")$output1
+        output_dims = "time", ncores = ncores)$output1
     }
     # REMEMBER to add na.rm = T in colMeans in .proxiesattractor 
-    res <- DynBiasCorrection(exp = exp_interpolated$data, obs = obs_ref, ...)
+    res <- DynBiasCorrection(exp = exp_interpolated$data, obs = obs_ref, ncores = ncores, ...)
   } 
  
   # Return a list of three elements
@@ -322,7 +330,3 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid,
 
   return(res)
 }
-
-
-
-
diff --git a/R/Interpolation.R b/R/Interpolation.R
index 1599bf3be04dae19b40ef65466b03bfd913e1d2f..0a53afaf66f865a7cdaf3129a711e2119e6b78b0 100644
--- a/R/Interpolation.R
+++ b/R/Interpolation.R
@@ -37,7 +37,8 @@
 #'@param method_point_interp a character vector indicating the interpolation method to
 #'interpolate model gridded data into the point locations. Accepted methods are "nearest",
 #'"bilinear", "9point", "invdist4nn", "NE", "NW", "SE", "SW".
-#'
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@seealso \code{\link[s2dverification]{CDORemap}}
 #'
 #'@return An s2dv object containing the dowscaled field.
@@ -52,7 +53,7 @@
 #'@export
 CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_grid = NULL, 
                               lat_dim = "lat", lon_dim = "lon", region = NULL, 
-                              method_point_interp = NULL)
+                              method_point_interp = NULL, ncores = NULL)
 {
   if (!inherits(exp,'s2dv_cube')) {
     stop("Parameter 'exp' must be of the class 's2dv_cube'")
@@ -70,7 +71,7 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr
   res <- Interpolation(exp = exp$data, lats = exp$coords[[lat_dim]], lons = exp$coords[[lon_dim]], 
                        source_file = exp$attrs$source_files[1], points = points,
                        method_remap = method_remap, target_grid = target_grid, lat_dim = lat_dim,
-                       lon_dim = lon_dim, region = region, method_point_interp = method_point_interp)
+                       lon_dim = lon_dim, region = region, method_point_interp = method_point_interp, ncores = ncores)
 
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   exp$data <- res$data
@@ -129,6 +130,8 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr
 #'interpolate model gridded data into the point locations. Accepted methods are "nearest",
 #'"bilinear", "9point", "invdist4nn", "NE", "NW", "SE", "SW". Only needed if the downscaling 
 #'is to a point location.
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@import multiApply
 #'@import plyr 
 #'@importFrom s2dv CDORemap
@@ -147,7 +150,7 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr
 #'@export
 Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, method_remap = NULL,  
                           target_grid = NULL, lat_dim = "lat", lon_dim = "lon", region = NULL, 
-                          method_point_interp = NULL) 
+                          method_point_interp = NULL, ncores = NULL) 
 {
   if (!is.null(method_remap)) {
     if (!inherits(method_remap, 'character')) {
@@ -232,7 +235,14 @@ Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, me
            "be passed to CDO, and it must be a grid recognised by CDO or a NetCDF file.")
     }
   }
-
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
   #----------------------------------
   # Limits of the region defined by the model data
   #----------------------------------
@@ -256,12 +266,27 @@ Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, me
   # Map regrid with CDO
   #----------------------------------
   if (is.null(points)) {
-    res <- CDORemap(data_array = exp, 
-                    lats = lats, 
-                    lons = lons, 
-                    grid = target_grid, 
-                    method = method_remap, 
-                    crop = region)
+    
+    .KnownLonNames <- s2dv:::.KnownLonNames
+    .KnownLatNames <- s2dv:::.KnownLatNames
+    .warning <- s2dv:::.warning
+    
+    res <- multiApply::Apply(data = list(data_array = exp), 
+                             target_dims = c(lon_dim,lat_dim),
+                             fun = CDORemap, 
+                             lats = lats, 
+                             lons = lons,
+                             grid = target_grid, 
+                             method = method_remap,
+                             crop = region, 
+                             ncores = ncores)
+
+    # res <- CDORemap(data_array = exp,
+    #                 lats = lats,
+    #                 lons = lons,
+    #                 grid = target_grid,
+    #                 method = method_remap,
+    #                 crop = region)
 
     # Return a list
     res <- list(data = res$data_array, obs = NULL, lon = res$lons, lat = res$lats)
@@ -271,16 +296,16 @@ Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, me
   #----------------------------------
   } else {  
     # First create interpolation weights, depending on the chosen method
-    weights <- create_interp_weights(ncfile = source_file, locids = 1:unique(lengths(points)),
+    weights <- .create_interp_weights(ncfile = source_file, locids = 1:unique(lengths(points)),
                                      lats = points[[lat_dim]], lons = points[[lon_dim]], 
                                      method = method_point_interp, region = list(lat_min = region[3], 
                                      lat_max = region[4], lon_min = region[1], lon_max = region[2])) 
 
     # Select coarse-scale data to be interpolated
-    model_data_gridpoints <- get_model_data(weights.df = weights, mdata = exp)
+    model_data_gridpoints <- .get_model_data(weights.df = weights, mdata = exp, ncores = ncores)
 
     # Interpolate model data to point locations
-    res <- interpolate_data(model_data_gridpoints, weights)
+    res <- .interpolate_data(model_data_gridpoints, weights, ncores = ncores)
   
     # Return a list
     res <- list(data = res, obs = NULL, lon = points[[lon_dim]], lat = points[[lat_dim]])
@@ -293,7 +318,7 @@ Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, me
 # Compute weights for interpolation at several (lat,lon) positions
 # We assume that grid boxes are centered in the grid point.
 #======================
-create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
+.create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
                                   method = c("nearest", "bilinear", "9point", "invdist4nn", "NE", 
                                              "NW", "SE", "SW")) 
 { 
@@ -307,7 +332,7 @@ create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
   #----------------
   # Read grid description and compute (i,j) of requested locations (including decimals)
   #----------------
-  griddes <- get_griddes(nc_cropped1)
+  griddes <- .get_griddes(nc_cropped1)
 
   if (is.null(griddes$yinc)) {
     system(paste0('rm ', nc_cropped1))
@@ -318,7 +343,7 @@ create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
   # If latitudes are decreasingly ordered, revert them
   if (griddes$yinc < 0) {
     system(paste0('cdo invertlat ', nc_cropped1, ' ', nc_cropped2))
-    griddes <- get_griddes(nc_cropped2)
+    griddes <- .get_griddes(nc_cropped2)
   }
   # remove temporary files
   system(paste0('rm ', nc_cropped1))
@@ -328,7 +353,7 @@ create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
     stop("'griddes' not found in the NetCDF source files")
   } 
   
-  gridpoints <- latlon2ij(griddes, lats, lons)
+  gridpoints <- .latlon2ij(griddes, lats, lons)
 
   #----------------
   # Compute the weights according to the selected method
@@ -573,7 +598,7 @@ create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
 # Works only for 'lonlat' and 'gaussian' grids.
 # Grids are supposed to cover whole globe.
 #======================
-latlon2ij <- function(griddes, lats, lons) {
+.latlon2ij <- function(griddes, lats, lons) {
   #------------
   # Check input params
   #------------
@@ -623,7 +648,7 @@ latlon2ij <- function(griddes, lats, lons) {
 #======================
 # Use cdo griddes to obtain grid information
 #======================
-get_griddes <- function(ncfile) {
+.get_griddes <- function(ncfile) {
   tmp <- system(paste0("cdo griddes ", ncfile,
                        " 2>/dev/null | egrep 'gridtype|xsize|ysize|xfirst|xinc|yfirst|yinc'"), intern = T)
   arr <- do.call(rbind, strsplit(tmp,"\\s+= ", perl = T))
@@ -631,7 +656,7 @@ get_griddes <- function(ncfile) {
   names(griddes) <- arr[,1]
 
   if(griddes$gridtype == "gaussian") {
-    griddes$yvals <- get_lats(ncfile)
+    griddes$yvals <- .get_lats(ncfile)
   }
 
   # Convert some fields to numeric. Ensures all fields are present.
@@ -645,7 +670,7 @@ get_griddes <- function(ncfile) {
 #======================
 # Use nco to obtain latitudes. Latitudes shall be named "lat" or "latitude".
 #======================
-get_lats <- function(ncfile) {
+.get_lats <- function(ncfile) {
 
   tmp <- system(paste0('ncks -H -s "%f " -v latitude ',ncfile),intern=T)
 
@@ -663,7 +688,7 @@ get_lats <- function(ncfile) {
 # Uses StartR. All ... parameters go to Start (i.e. specify dat, var, 
 # sdate, time, ensemble, num_procs, etc)
 #======================
-get_model_data <- function(weights.df, mdata) {
+.get_model_data <- function(weights.df, mdata, ncores = NULL) {
  
   #-----------------
   # Get data for all combinations of i and j.
@@ -705,7 +730,9 @@ get_model_data <- function(weights.df, mdata) {
   #-----------------
   # Retrieve with multiApply
   #-----------------
-  sub_mdata <- Apply(mdata, target_dims = list(c(latdim, londim)), fun = function(x) {laply(1:length(is),function(k) { x[js[k],is[k]] }) })$output1
+  sub_mdata <- Apply(mdata, target_dims = list(c(latdim, londim)), 
+                     fun = function(x) {laply(1:length(is),function(k) { x[js[k],is[k]] }) }, 
+                     ncores = ncores)$output1
   names(dim(sub_mdata))[1] <- "gridpoint"
 
   #-----------------
@@ -718,7 +745,7 @@ get_model_data <- function(weights.df, mdata) {
 # Multiply the grid-point series by the weights, 
 # to obtain the desired interpolations
 #======================
-interpolate_data <- function(model_data, weights.df) {
+.interpolate_data <- function(model_data, weights.df, ncores) {
   #-----------------
   # Multiply each gridpoint matrix by its corresponding weight
   #-----------------
@@ -732,9 +759,8 @@ interpolate_data <- function(model_data, weights.df) {
   #interp_data <- apply(weighted_data, -gpdim, function(x) { rowsum(x, weights.df$locid) }) 
   #names(dim(interp_data))[1] <- "location"
   interp_data <- Apply(weighted_data, target_dims = gpdim, fun = function(x) { 
-    rowsum(x, weights.df$locid)}, output_dims = c("location", "aux"))$output1
+    rowsum(x, weights.df$locid)}, output_dims = c("location", "aux"),
+    ncores = ncores)$output1
   
   return(interp_data)
 }
-
-
diff --git a/R/Intlr.R b/R/Intlr.R
index 5ac11c3d7a68bb0e909b35d5adec722e06a239f7..b4b8a75e4ae20d086c2bb09010ae7c9978a41d12 100644
--- a/R/Intlr.R
+++ b/R/Intlr.R
@@ -394,7 +394,7 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
   # Predictand: observations
   else if (lr_method == '4nn') {
     predictor <- .find_nn(coar = exp, lats_hres = obs_lats, lons_hres = obs_lons, lats_coar = exp_lats, 
-                          lons_coar = exp_lons, lat_dim = lat_dim, lon_dim = lon_dim, nn = 4, ncores = ncores)
+                         lons_coar = exp_lons, lat_dim = lat_dim, lon_dim = lon_dim, nn = 4, ncores = ncores)
     
     if (is.null(points)) {
       if (!is.null(target_grid)) {
diff --git a/R/LogisticReg.R b/R/LogisticReg.R
index d58e8ab0eeaf47fc088964f04d0ca87aa8e6d710..f569610f885941b6b9a9b8b85fa213c89dac1214 100644
--- a/R/LogisticReg.R
+++ b/R/LogisticReg.R
@@ -62,8 +62,8 @@
 #'takes the first and last elements of the latitudes and longitudes.
 #'@param loocv a logical vector indicating whether to perform leave-one-out cross-validation
 #'in the fitting of the logistic regression. Default to FALSE.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
-#'
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@import multiApply
 #'@import nnet
 #'@importFrom laply plyr  
@@ -91,7 +91,7 @@
 CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_method = "ens_mean",
                   probs_cat = c(1/3,2/3), return_most_likely_cat = FALSE, points = NULL, 
                   method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", 
-                  member_dim = "member", region = NULL, loocv = FALSE, ncores = 1) {
+                  member_dim = "member", region = NULL, loocv = FALSE, ncores = NULL) {
 
   if (!inherits(exp,'s2dv_cube')) {
     stop("Parameter 'exp' must be of the class 's2dv_cube'")
@@ -198,8 +198,8 @@ CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_me
 #'takes the first and last elements of the latitudes and longitudes.
 #'@param loocv a logical vector indicating whether to perform leave-one-out cross-validation
 #'in the fitting of the logistic regression. Default to FALSE.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
-#'
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@import multiApply
 #'@import nnet
 #'@importFrom laply plyr  
@@ -227,7 +227,7 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
                   int_method = NULL, log_reg_method = "ens_mean", probs_cat = c(1/3,2/3), 
                   return_most_likely_cat = FALSE, points = NULL, method_point_interp = NULL, 
                   lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", member_dim = "member", 
-                  source_file = NULL, region = NULL, loocv = FALSE, ncores = 1) {
+                  source_file = NULL, region = NULL, loocv = FALSE, ncores = NULL) {
 
   #-----------------------------------
   # Checkings
@@ -303,7 +303,14 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
             "'obs_lats' and 'obs_lons'.")
     region <- c(obs_lons[1], obs_lons[length(obs_lons)], obs_lats[1], obs_lats[length(obs_lats)])
   }
-
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
   # the code is not yet prepared to handle members in the observations
   restore_ens <- FALSE
   if (member_dim %in% names(dim(obs))) {
@@ -319,18 +326,18 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
   exp_interpolated <- Interpolation(exp = exp, lats = exp_lats, lons = exp_lons, target_grid = target_grid,
                                     method_remap = int_method, points = points, source_file = source_file,
                                     lat_dim = lat_dim, lon_dim = lon_dim, method_point_interp = method_point_interp,
-                                    region = region)
+                                    region = region, ncores = ncores)
 
   # compute ensemble mean anomalies
   if (log_reg_method == "ens_mean") {
-    predictor <- get_ens_mean_anom(obj_ens = exp_interpolated$data, member_dim = member_dim, sdate_dim = sdate_dim)
+    predictor <- .get_ens_mean_anom(obj_ens = exp_interpolated$data, member_dim = member_dim, sdate_dim = sdate_dim, ncores = ncores)
 
     target_dims_predictor <- sdate_dim
   } 
   else if (log_reg_method == "ens_mean_sd")  {
-    ens_mean_anom <- get_ens_mean_anom(obj_ens = exp_interpolated$data, member_dim = member_dim, 
-                                       sdate_dim = sdate_dim)
-    ens_sd <- get_ens_sd(obj_ens = exp_interpolated$data, member_dim = member_dim)
+    ens_mean_anom <- .get_ens_mean_anom(obj_ens = exp_interpolated$data, member_dim = member_dim, 
+                                       sdate_dim = sdate_dim, ncores = ncores)
+    ens_sd <- .get_ens_sd(obj_ens = exp_interpolated$data, member_dim = member_dim, ncores = ncores)
 
     #merge two arrays into one array of predictors
     predictor <- abind(ens_mean_anom, ens_sd, along = 1/2)
@@ -338,7 +345,7 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
      
     target_dims_predictor <- c(sdate_dim, "pred")
   } else if (log_reg_method == "sorted_members") {   
-    predictor <- sort_members(obj_ens = exp_interpolated$data, member_dim = member_dim)
+    predictor <- .sort_members(obj_ens = exp_interpolated$data, member_dim = member_dim, ncores = ncores)
 
     target_dims_predictor <- c(sdate_dim, member_dim)
   } else {
@@ -352,7 +359,7 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
     obs_interpolated <- Interpolation(exp = obs, lats = obs_lats, lons = obs_lons, target_grid = target_grid,
                                       method_remap = int_method, points = points, source_file = source_file,
                                       lat_dim = lat_dim, lon_dim = lon_dim,
-                                      method_point_interp = method_point_interp, region = region)
+                                      method_point_interp = method_point_interp, region = region, ncores = ncores)
     obs_ref <- obs_interpolated$data
   } else {
     obs_ref <- obs
@@ -362,19 +369,19 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
   obs_cat <- Apply(obs_ref, target_dims = sdate_dim, function(x) {
                    terc <- convert2prob(as.vector(x), prob = probs_cat) 
                    apply(terc, 1, function(r) which (r == 1))},          
-               output_dims = sdate_dim)$output1
+               output_dims = sdate_dim, ncores = ncores)$output1
 
   res <- Apply(list(predictor, obs_cat), target_dims = list(target_dims_predictor, sdate_dim), fun = function(x, y)
-               .log_reg(x = x, y = y, loocv = loocv,probs_cat=probs_cat), output_dims = c(sdate_dim, "category"))$output1
+               .log_reg(x = x, y = y, loocv = loocv), output_dims = c(sdate_dim, "category"), ncores = ncores)$output1
 
   if (return_most_likely_cat) {
-    res <- Apply(res, target_dims = c(sdate_dim, "category"), most_likely_category,
-                 output_dims = sdate_dim)$output1
+    res <- Apply(res, target_dims = c(sdate_dim, "category"), .most_likely_category,
+                 output_dims = sdate_dim, ncores = ncores)$output1
   }
 
   # restore ensemble dimension in observations if it existed originally
   if (restore_ens) {
-    obs_ref <- s2dv::InsertDim(obs_ref, posdim = 1, lendim = 1, name = member_dim)
+    obs_ref <- s2dv::InsertDim(obs_ref, posdim = 1, lendim = 1, name = member_dim, ncores = ncores)
   }
 
   res <- list(data = res, obs = obs_ref, lon = exp_interpolated$lon, lat = exp_interpolated$lat)
@@ -382,7 +389,7 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
   return(res)
 }
 
-most_likely_category <- function(data) {
+.most_likely_category <- function(data) {
 # data, expected dims: start date, category (in this order)
 
   if (all(is.na(data))) {
@@ -393,33 +400,33 @@ most_likely_category <- function(data) {
   return(mlc)
 }
 
-sort_members <- function(obj_ens, member_dim) {
+.sort_members <- function(obj_ens, member_dim, ncores = NULL) {
 
-  sorted <- Apply(obj_ens, target_dims = member_dim, sort, decreasing = TRUE, na.last = TRUE)$output1  
+  sorted <- Apply(obj_ens, target_dims = member_dim, sort, decreasing = TRUE, na.last = TRUE, ncores = ncores)$output1  
 
   return(sorted)
 }
 
-get_ens_sd <- function(obj_ens, member_dim) {
+.get_ens_sd <- function(obj_ens, member_dim, ncores = NULL) {
 
   # compute ensemble spread
-  ens_sd <- Apply(obj_ens, target_dims = member_dim, sd, na.rm = TRUE)$output1
+  ens_sd <- Apply(obj_ens, target_dims = member_dim, sd, na.rm = TRUE, ncores = ncores)$output1
 
   return(ens_sd)
 }
 
-get_ens_mean_anom <- function(obj_ens, member_dim, sdate_dim) {
+.get_ens_mean_anom <- function(obj_ens, member_dim, sdate_dim, ncores = NULL) {
 
   require(s2dv)
 
   # compute climatology
-  clim <- Apply(obj_ens, target_dims = c(member_dim, sdate_dim), mean)$output1
+  clim <- Apply(obj_ens, target_dims = c(member_dim, sdate_dim), mean, ncores = ncores)$output1
 
   # compute ensemble mean
-  ens_mean <- Apply(obj_ens, target_dims = member_dim, mean, na.rm = TRUE)$output1
+  ens_mean <- Apply(obj_ens, target_dims = member_dim, mean, na.rm = TRUE, ncores = ncores)$output1
 
   # compute ensemble mean anomalies
-  anom <- Ano(ens_mean, clim)
+  anom <- Ano(ens_mean, clim, ncores = ncores)
 
   return(anom)
 }
@@ -437,7 +444,7 @@ get_ens_mean_anom <- function(obj_ens, member_dim, sdate_dim) {
 
   } else {
     # training
-    lm1 <- train_lr(df = tmp_df, loocv = loocv)
+    lm1 <- .train_lr(df = tmp_df, loocv = loocv)
 
     # prediction
     res <- pred_lr(lm1 = lm1, df = tmp_df, loocv = loocv,probs_cat=probs_cat)
@@ -448,7 +455,7 @@ get_ens_mean_anom <- function(obj_ens, member_dim, sdate_dim) {
 #-----------------------------------
 # Function to train the logistic regressions.
 #-----------------------------------
-train_lr <- function(df, loocv) {
+.train_lr <- function(df, loocv) {
 
   require(nnet)
   
@@ -534,5 +541,3 @@ pred_lr <- function(df, lm1, loocv,probs_cat) {
 
   return(pred_vals)
 }
-
-