diff --git a/NAMESPACE b/NAMESPACE
index 9f9a0e3481c2eaa1eb04b26776b497b25783564e..efc5c8353c2d2e8a9f20f8e3980faf5d1d2027b1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -85,6 +85,7 @@ importFrom(plyr,.)
 importFrom(plyr,dlply)
 importFrom(reshape2,melt)
 importFrom(s2dv,Reorder)
+importFrom(s2dverification,Subset)
 importFrom(utils,glob2rx)
 importFrom(utils,head)
 importFrom(utils,tail)
diff --git a/NEWS.md b/NEWS.md
index bcdf86dd449a5db10564103f38653464065cd1e4..c8dfc43f5e9d30c271aa6d240e53977051796fc7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
 **Submission date to CRAN: XX-12-2020**
 - New features:
     + PlotPDFsOLE includes parameters to modify legend style
-
+    + CST_RFSlope handless missing values in the temporal dimension and new ncores parameter for parallel computation
 - Fixes:
     + PlotForecastPDF correctly displays terciles labels 
 
diff --git a/R/CST_RFSlope.R b/R/CST_RFSlope.R
index 6bc3dc2d989d373b5f59912b72f7abf7af03fb2c..759820377b10fd520ac8f09c74321698ae5e3c08 100644
--- a/R/CST_RFSlope.R
+++ b/R/CST_RFSlope.R
@@ -15,6 +15,7 @@
 #' over which to compute spectral slopes. If a character array of dimension names is provided, the spectral slopes
 #' will be computed as an average over all elements belonging to those dimensions.
 #' If omitted one of c("ftime", "sdate", "time")  is searched and the first one with more than one element is chosen.
+#' @param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one.
 #' @return CST_RFSlope() returns spectral slopes using the RainFARM convention
 #' (the logarithmic slope of k*|A(k)|^2 where A(k) are the spectral amplitudes).
 #' The returned array has the same dimensions as the \code{exp} element of the input object,
@@ -38,7 +39,7 @@
 #' #[1,] 1.893503 1.893503 1.893503
 #' #[2,] 1.893503 1.893503 1.893503
 #' @export
-CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) {
+CST_RFSlope <- function(data, kmin = 1, time_dim = NULL, ncores = 1) {
 
   slopes <- RFSlope(data$data, kmin, time_dim,
                     lon_dim = "lon", lat_dim = "lat")
@@ -68,11 +69,15 @@ CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) {
 #' with more than one element is chosen.
 #' @param lon_dim Name of lon dimension ("lon" by default).
 #' @param lat_dim Name of lat dimension ("lat" by default).
+#' @param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one.
+#'
 #' @return RFSlope() returns spectral slopes using the RainFARM convention
 #' (the logarithmic slope of k*|A(k)|^2 where A(k) are the spectral amplitudes).
 #' The returned array has the same dimensions as the input array,
 #' minus the dimensions specified by \code{lon_dim}, \code{lat_dim} and \code{time_dim}.
 #' @import multiApply
+#' @import rainfarmr
+#' @importFrom s2dverification Subset
 #' @export
 #' @examples
 #' # Example for the 'reduced' RFSlope function 
@@ -94,8 +99,17 @@ CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) {
 #' #[3,] 1.893503 1.893503 1.893503
 #' #[4,] 1.893503 1.893503 1.893503
 RFSlope <- function(data, kmin = 1, time_dim = NULL,
-                    lon_dim = "lon", lat_dim = "lat") {
-
+                    lon_dim = "lon", lat_dim = "lat", ncores = 1) {
+  if (length(ncores) > 1) {
+    ncores = ncores[1]
+    warning("Parameter 'ncores' has length > 1 and only the first element will be used.")
+  }
+  if (!is.null(ncores)) {
+    ncores <- round(ncores)
+    if (ncores == 0) {
+      ncores = NULL
+    }
+  }
   # Ensure input  grid is square and with even dimensions
   if ( (dim(data)[lon_dim] != dim(data)[lat_dim]) |
        (dim(data)[lon_dim] %% 2 == 1)) {
@@ -139,7 +153,7 @@ RFSlope <- function(data, kmin = 1, time_dim = NULL,
 
   # Repeatedly apply .RFSlope
   result <- Apply(data, c(lon_dim, lat_dim, "rainfarm_samples"),
-                  .RFSlope, kmin)$output1
+                  .RFSlope, kmin, ncores = ncores)$output1
 
   return(slopes = result)
 }
@@ -152,7 +166,11 @@ RFSlope <- function(data, kmin = 1, time_dim = NULL,
 #' @noRd
 
 .RFSlope <- function(pr, kmin) {
-
+  if (any(is.na(pr))) {
+    posna <- unlist(lapply(1:dim(pr)['rainfarm_samples'],
+                         function(x){!is.na(pr[1, 1, x])}))
+    pr <- Subset(pr, 'rainfarm_samples', posna)
+  }
   fxp <- fft2d(pr)
   sx <- fitslope(fxp, kmin = kmin)
   return(sx)
diff --git a/man/CST_RFSlope.Rd b/man/CST_RFSlope.Rd
index 0c4e16710c82a98b1d1245a2321f4675cd5a578e..b76ac93ef98a9ee157a71a664ab86a24824c84c3 100644
--- a/man/CST_RFSlope.Rd
+++ b/man/CST_RFSlope.Rd
@@ -4,7 +4,7 @@
 \alias{CST_RFSlope}
 \title{RainFARM spectral slopes from a CSTools object}
 \usage{
-CST_RFSlope(data, kmin = 1, time_dim = NULL)
+CST_RFSlope(data, kmin = 1, time_dim = NULL, ncores = 1)
 }
 \arguments{
 \item{data}{An object of the class 's2dv_cube', containing the spatial precipitation fields to downscale.
@@ -18,6 +18,8 @@ to average these slopes, which can be specified by parameter \code{time_dim}.}
 over which to compute spectral slopes. If a character array of dimension names is provided, the spectral slopes
 will be computed as an average over all elements belonging to those dimensions.
 If omitted one of c("ftime", "sdate", "time")  is searched and the first one with more than one element is chosen.}
+
+\item{ncores}{is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one.}
 }
 \value{
 CST_RFSlope() returns spectral slopes using the RainFARM convention
diff --git a/man/RFSlope.Rd b/man/RFSlope.Rd
index db3f0e1054ff5a4b93289a453c5d1bd272a4f9e8..5c0c168955e9aef866a0e44c31b302e9e66756ca 100644
--- a/man/RFSlope.Rd
+++ b/man/RFSlope.Rd
@@ -4,7 +4,14 @@
 \alias{RFSlope}
 \title{RainFARM spectral slopes from an array (reduced version)}
 \usage{
-RFSlope(data, kmin = 1, time_dim = NULL, lon_dim = "lon", lat_dim = "lat")
+RFSlope(
+  data,
+  kmin = 1,
+  time_dim = NULL,
+  lon_dim = "lon",
+  lat_dim = "lat",
+  ncores = 1
+)
 }
 \arguments{
 \item{data}{Array containing the spatial precipitation fields to downscale.
@@ -25,6 +32,8 @@ with more than one element is chosen.}
 \item{lon_dim}{Name of lon dimension ("lon" by default).}
 
 \item{lat_dim}{Name of lat dimension ("lat" by default).}
+
+\item{ncores}{is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one.}
 }
 \value{
 RFSlope() returns spectral slopes using the RainFARM convention