diff --git a/DESCRIPTION b/DESCRIPTION
index f8e69b11f56e9726be4263064cfa5fb96ee411a7..d7691f7166a63a8c3f07280ad393e0da61963ffe 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,14 +1,21 @@
 Package: ClimProjDiags
 Title: Set of Tools to Compute Various Climate Indices
-Version: 0.0.3
+Version: 0.0.4
 Authors@R: c(
     person("BSC-CNS", role = c("aut", "cph")),
     person("Alasdair", "Hunter", , "alasdair.hunter@bsc.es", role = "aut"),   
     person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre")),
     person("Nicolau", "Manubens", , "nicolau.manubens@bsc.es", role = "ctb"),
     person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "ctb"))
-Description: Set of tools to compute metrics and indices for climate analysis. The package provides functions to compute extreme indices, evaluate the agreement between models and combine theses models into an ensemble. Multi-model time series of climate indices can be computed either after averaging the 2-D fields from different models provided they share a common grid or by combining time series computed on the model native grid. Indices can be assigned weights and/or combined to construct new indices.
-Depends: R (>= 3.2.0)
+Description: Set of tools to compute metrics and indices for climate analysis.
+    The package provides functions to compute extreme indices, evaluate the
+    agreement between models and combine theses models into an ensemble. Multi-model
+    time series of climate indices can be computed either after averaging the 2-D
+    fields from different models provided they share a common grid or by combining
+    time series computed on the model native grid. Indices can be assigned weights
+    and/or combined to construct new indices.
+Depends:
+    R (>= 3.2.0)
 Imports:
     multiApply (>= 2.0.0),
     climdex.pcic,
@@ -20,8 +27,8 @@ URL: https://earth.bsc.es/gitlab/es/ClimProjDiags
 BugReports: https://earth.bsc.es/gitlab/es/ClimProjDiags/issues
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 6.0.1.9000
-Suggests: 
+RoxygenNote: 5.0.0
+Suggests:
     knitr,
     rmarkdown
-VignetteBuilder: knitr 
+VignetteBuilder: knitr
diff --git a/R/Threshold.R b/R/Threshold.R
index aa43b3a0f989c05c11b1c0db84bd4a986eed0b71..311163a07747f93b9b97ab677b70583831372ffa 100644
--- a/R/Threshold.R
+++ b/R/Threshold.R
@@ -8,6 +8,7 @@
 #'@param base.range The years used for computing the threshold.
 #'@param qtiles Numeric vector with values between 0 and 1 indicating the quantiles to be computed.
 #'@param ncores The number of cores to be used when computing the threshold.
+#'@param na.rm A logical value. If TRUE, any NA and NaN's are removed before the quantiles are computed (default as FALSE).
 #'
 #'@return An array with similar dimensions as the \code{data} input, but without 'time' dimension, and a new 'jdays' dimension.
 #'
@@ -29,7 +30,7 @@
 #'a <- Threshold(data, dates = NULL, base.range = NULL, qtiles = 0.9, ncores = NULL)
 #'str(a)
 #'@export
-Threshold <- function(data, dates = NULL, calendar = NULL, base.range = NULL, qtiles = 0.9, ncores = NULL) {
+Threshold <- function(data, dates = NULL, calendar = NULL, base.range = NULL, qtiles = 0.9, ncores = NULL, na.rm = FALSE) {
   if (is.null(data)) {
     stop("Parameter 'data' cannot be NULL.")
   }
@@ -93,6 +94,14 @@ Threshold <- function(data, dates = NULL, calendar = NULL, base.range = NULL, qt
   if (stop_error) {
     stop("Parameter 'dates' must be of the same length as the 'time' dimension of the parameter 'data'.")
   }
+  if (!is.logical(na.rm)) {
+    stop("Parameter 'na.rm' must be logical.")
+  }
+  if (length(na.rm) > 1) {
+    na.rm <- na.rm[1]
+    warning("Parameter 'na.rm' has length > 1 and only the first ",
+            "element will be used.")
+  }
   dates <- as.PCICt(dates, cal = calendar)
   dates = as.character(dates)
   jdays <- as.numeric(strftime(dates, format = "%j"))
@@ -170,7 +179,7 @@ Threshold <- function(data, dates = NULL, calendar = NULL, base.range = NULL, qt
   if (length(dim(data)) > 1) {
     result <- Apply(data = data, margins = margins, 
                     fun = .Threshold, indices = jdays, qtiles = qtiles, 
-                    ncores = ncores)
+                    ncores = ncores, na.rm = na.rm)
     names(dim(result$output1)) <- c("jdays", dim_names[-time_dim])
   } else {
     result <- list() 
@@ -179,6 +188,6 @@ Threshold <- function(data, dates = NULL, calendar = NULL, base.range = NULL, qt
   }
   return(result$output1)
 }
-.Threshold <- function(data, indices,  qtiles) {
-  tapply(X = data, INDEX = indices, FUN = quantile, probs = qtiles)
+.Threshold <- function(data, indices,  qtiles, na.rm = FALSE) {
+  tapply(X = data, INDEX = indices, FUN = quantile, probs = qtiles, na.rm = na.rm)
 } 
diff --git a/man/AnoAgree.Rd b/man/AnoAgree.Rd
index 63c1450039013f94bee9c4e66d29c0923fcdab84..c929352d71c27790644570ed71b937a22a7072d9 100644
--- a/man/AnoAgree.Rd
+++ b/man/AnoAgree.Rd
@@ -34,3 +34,4 @@ a <- rnorm(6)
 agree <- AnoAgree(ano = a, membersdim = 1, na.rm = TRUE, ncores = NULL)
 print(agree)
 }
+
diff --git a/man/Climdex.Rd b/man/Climdex.Rd
index 12ebac5e16cb1c3092843fac2f42372948c189de..3ff6bb76494faed3037cdf06037f308802799bf5 100644
--- a/man/Climdex.Rd
+++ b/man/Climdex.Rd
@@ -63,3 +63,4 @@ str(thres)
 clim <- Climdex(data, metric = "t90p", threshold = thres)
 str(clim)
 }
+
diff --git a/man/CombineIndices.Rd b/man/CombineIndices.Rd
index 6d032cddcba684a8375572aa6c259da2c9d5e982..95abc5ed4f5564527b93bd9646221b7c4b1a8d3d 100644
--- a/man/CombineIndices.Rd
+++ b/man/CombineIndices.Rd
@@ -33,3 +33,4 @@ dim(b) <- c(lon = 2, lat = 3, mod = 4)
 comb_ind <- CombineIndices(indices = list(a, b), weights = c(2, 1), operation = "add")
 print(comb_ind)
 }
+
diff --git a/man/DTRIndicator.Rd b/man/DTRIndicator.Rd
index f46d34d75cec99190a4249cffdb3a61f48e7a291..33e7f507882ba4583c5796645b05991a1d93af22 100644
--- a/man/DTRIndicator.Rd
+++ b/man/DTRIndicator.Rd
@@ -60,3 +60,4 @@ aa <- DTRIndicator(tmax, tmin, ref = a, by.seasons = FALSE, ncores = NULL)
 str(aa)
 dim(aa$indicator)
 }
+
diff --git a/man/DTRRef.Rd b/man/DTRRef.Rd
index 1239eea8a778a097d22d892887aab5c327f5b158..6e32e31a0085e87305b3e30c601dec7fe28fa6bc 100644
--- a/man/DTRRef.Rd
+++ b/man/DTRRef.Rd
@@ -68,3 +68,4 @@ dim(tmin) <- c(2, 3, 365)
 a <- DTRRef(tmax, tmin, by.seasons = FALSE, dates = time,  timedim = 3, ncores = NULL)
 str(a)
 }
+
diff --git a/man/DailyAno.Rd b/man/DailyAno.Rd
index dcd53c37459889b5610d15c5524ac93e2345a6fe..d4033f79e505fbd3df46a6304be677212c811f1b 100644
--- a/man/DailyAno.Rd
+++ b/man/DailyAno.Rd
@@ -31,3 +31,4 @@ jdays <- c(rep(1, 5), rep(2, 5))
 daily_anomaly <- DailyAno(data = data, jdays = jdays, na.rm = TRUE)
 print(daily_anomaly)
 }
+
diff --git a/man/Extremes.Rd b/man/Extremes.Rd
index 4ce49cf700146e6412fbab5ee1fc15956441c219..0939fa67275433b18b9f62df84150911ea72e352 100644
--- a/man/Extremes.Rd
+++ b/man/Extremes.Rd
@@ -58,3 +58,4 @@ a <- Extremes(data, threshold = threshold, op = ">", min.length = 6, spells.can.
              max.missing.days = 5, ncores = NULL)
 str(a)
 }
+
diff --git a/man/SeasonSelect.Rd b/man/SeasonSelect.Rd
index 33066eed0e4ce2cd5b1a1abbe92930b43854975c..a71fd242c9a71d0647fa10fad5fc1a440823af70 100644
--- a/man/SeasonSelect.Rd
+++ b/man/SeasonSelect.Rd
@@ -45,3 +45,4 @@ attr(data, 'Variables')$common[[2]]$dim[[3]]$vals <- time
 a <- SeasonSelect(data = data, season = 'JJA')
 str(a)
 }
+
diff --git a/man/SelBox.Rd b/man/SelBox.Rd
index 38e3547bf6ae6d6e699e78821aaefb1311896535..07e92b20e28a41e872059f311378f8e1f90b7573 100644
--- a/man/SelBox.Rd
+++ b/man/SelBox.Rd
@@ -43,3 +43,4 @@ a <- SelBox(data = data, lon = lon, lat = lat, region = c(2, 20, 1, 5),
            londim = 1, latdim = 2, mask = NULL)
 str(a)
 }
+
diff --git a/man/Subset.Rd b/man/Subset.Rd
index 5d24310077de9c6af93a3c10bf1c26fcb6f33240..634cfe2a60d07b8d31eeafa4885fca01afd71b42 100644
--- a/man/Subset.Rd
+++ b/man/Subset.Rd
@@ -31,3 +31,4 @@ data_subset <- Subset(data, c('time', 'model'),
 dim(data_subset)
 
 }
+
diff --git a/man/Threshold.Rd b/man/Threshold.Rd
index bcd630091cfeb6620c42eceb9902b0fa6a4e3aa5..cec12e68bc83ea93767f034c6aa417b5076e81b0 100644
--- a/man/Threshold.Rd
+++ b/man/Threshold.Rd
@@ -5,7 +5,7 @@
 \title{Daily thresholds based on quantiles for n-dimensional arrays}
 \usage{
 Threshold(data, dates = NULL, calendar = NULL, base.range = NULL,
-  qtiles = 0.9, ncores = NULL)
+  qtiles = 0.9, ncores = NULL, na.rm = FALSE)
 }
 \arguments{
 \item{data}{A numeric n-dimensional array containing daily data.}
@@ -19,6 +19,8 @@ Threshold(data, dates = NULL, calendar = NULL, base.range = NULL,
 \item{qtiles}{Numeric vector with values between 0 and 1 indicating the quantiles to be computed.}
 
 \item{ncores}{The number of cores to be used when computing the threshold.}
+
+\item{na.rm}{A logical value. If TRUE, any NA and NaN's are removed before the quantiles are computed (default as FALSE).}
 }
 \value{
 An array with similar dimensions as the \code{data} input, but without 'time' dimension, and a new 'jdays' dimension.
@@ -40,3 +42,4 @@ attr(data, 'Variables')$dat1$time <- time
 a <- Threshold(data, dates = NULL, base.range = NULL, qtiles = 0.9, ncores = NULL)
 str(a)
 }
+
diff --git a/man/WaveDuration.Rd b/man/WaveDuration.Rd
index ba5ace2ddab9ff827fbfc6b0a8f35900391f9bda..b0682523c96b1d604c09c6c397b8040523c7e4e7 100644
--- a/man/WaveDuration.Rd
+++ b/man/WaveDuration.Rd
@@ -48,3 +48,4 @@ threshold <- rep(40, 31)
 a <- WaveDuration(data, threshold, op = ">", spell.length = 6, by.seasons = TRUE, ncores = NULL)
 str(a)
 }
+
diff --git a/man/WeightedMean.Rd b/man/WeightedMean.Rd
index 0141c70d5019781985d2fcdddf40ef79a608aa83..23f31c6c415b83ecd829fba0cc27c3c82d750504 100644
--- a/man/WeightedMean.Rd
+++ b/man/WeightedMean.Rd
@@ -60,3 +60,4 @@ a <- WeightedMean(data = data, lon = lon, lat = lat, region = NULL,
                  mask = NULL, londim = 1, latdim = 2)
 str(a)
 }
+