diff --git a/.Rbuildignore b/.Rbuildignore
index 0954c2a0092c62f36e7508d2d58568d39b20cdf1..d77251231585568f613aaf37eb80927a65d6740e 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,6 +1,8 @@
-.git
-.gitignore
-.tar.gz
-.pdf
-./.nc
+.*\.git$
+.*\.gitignore$
+.*\.tar.gz$
+.*\.pdf$
+./.nc$
+.*^(?!data)\.RData$
+.*\.gitlab-ci.yml$
 
diff --git a/.gitignore b/.gitignore
index b26d556a841a100900ec7e4824aded4dd3e3b82a..35a09175a4a4e587947985914ec082148bc51880 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,4 +15,6 @@ master_pull.txt
 *.ps
 Rplots.pdf
 .nfs*
+*.RData
+!data/*.RData
 
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000000000000000000000000000000000000..d7d860551ebf3b2eff156461c6a26540e23f0709
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,10 @@
+stages:
+  - build
+build:
+  stage: build
+  script:
+    - module load R
+    - R CMD build --resave-data .
+    - R CMD check --as-cran --no-manual ClimProjDiags_*.tar.gz
+    - R -e 'covr::package_coverage()'
+
diff --git a/DESCRIPTION b/DESCRIPTION
index 4fa1039c05e2cdb47b6a282416dad5ae18cccfe7..1776287eeb07b6c166b7d515f1514438791d974a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,14 +1,22 @@
 Package: ClimProjDiags
 Title: Set of Tools to Compute Various Climate Indices
-Version: 0.0.2
+Version: 1.0.0
 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)
+    person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "ctb"),
+    person("An-Chi", "Ho", , "an.ho@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)
 Imports:
     multiApply (>= 2.0.0),
     climdex.pcic,
@@ -20,8 +28,9 @@ 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,
+    testthat,
     rmarkdown
-VignetteBuilder: knitr 
+VignetteBuilder: knitr
diff --git a/NAMESPACE b/NAMESPACE
index 8c0c72b49711d1083ace038843337edfa80dfba2..2b91ef634043ccf05154bf1eac41ac4b16ff94b7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,12 +5,12 @@ export(Climdex)
 export(CombineIndices)
 export(DTRIndicator)
 export(DTRRef)
-export(DailyAno)
 export(Extremes)
 export(SeasonSelect)
 export(SelBox)
 export(Subset)
 export(Threshold)
+export(TimeAnomaly)
 export(WaveDuration)
 export(WeightedMean)
 import(PCICt)
diff --git a/R/DailyAno.R b/R/DailyAno.R
deleted file mode 100644
index d3df4cf3aba4f7d6900e1bcbaab056925a92dbe8..0000000000000000000000000000000000000000
--- a/R/DailyAno.R
+++ /dev/null
@@ -1,85 +0,0 @@
-#'Daily anomalies
-#'
-#'@description This function computes daily anomalies from a vector containing the daily time series.
-#'
-#'@param data A vector of daily data.
-#'@param jdays A vector of the corresponding day of the year. This vector must be the same length as parameter \code{data}.
-#'@param dates If \code{jdays} is not supplied, a vector of dates corresponding to the observations in \code{data} with defined calendar attributes.
-#'@param calendar A character indicating the calendar type.
-#'@param na.rm A logical indicating whether missing values should be removed. If \code{na.rm} is FALSE an NA value in any of the arguments will cause a value of NA to be returned, otherwise (TRUE by default) NA values are ignored.
-#'
-#'@return A vector of daily anomalies of the same length as parameter \code{data}.
-#'
-#'@examples
-#'# Time series in a vector example:
-#'data <- 1:10
-#'jdays <- c(rep(1, 5), rep(2, 5))
-#'daily_anomaly <- DailyAno(data = data, jdays = jdays, na.rm = TRUE)
-#'print(daily_anomaly)
-#'@export
-DailyAno <- function(data, jdays = NULL, dates = NULL, calendar = NULL, na.rm = TRUE) {
-  if (is.null(data)) {
-    stop("Parameters 'data' cannot be NULL.")
-  }
-  if (!is.vector(data)) {
-    stop("Parameters 'data' and 'jdays' must be a vector.")
-  }
-  if (is.null(jdays) & is.null(dates)) {
-    stop("At least one of the parameters 'jdays' or 'dates' must be supplied.")
-  }
-  if (is.null(jdays)) {
-    if (is.null(calendar)) {
-      calendar <- attributes(dates)$calendar
-      if (is.null(calendar)) {
-        calendar <- attributes(dates)$variables$time$calendar
-      }
-      if (is.null(calendar)) {
-        stop("The attribute 'calendar' must be present in the parameter 'dates' or specified in parameter 'calendar'.")
-      }
-      # (end identifying)
-    }  
-    if (!any(class(dates) %in% c('POSIXct'))) {
-      dates <- try( {
-        if (is.character(dates)) {
-          as.POSIXct(dates, format = "%Y%m%d")
-        } else {
-          as.POSIXct(dates)
-        }
-      })
-      if ('try-error' %in% class(dates) | sum(is.na(dates)) == length(dates)) {
-        stop("Dates provided in parameter 'dates' must be of class 'POSIXct' or convertable to 'POSIXct'.")
-      }
-    }
-    dates <- as.PCICt(dates, cal = calendar)
-    dates = as.character(dates)
-    jdays <- as.numeric(strftime(dates, format = "%j"))
-    if (calendar == "gregorian" | calendar == "standard" | calendar == "proleptic_gregorian") {
-      year <- as.numeric(strftime(dates, format = "%Y"))
-      if (length(unique(year)) > 1) {
-        pos <- ((year / 100) %% 1 == 0) + ((year / 4) %% 1 == 0) + ((year / 400) %% 1 == 0)
-        pos <- which(pos == 0 | pos == 2 | pos == 4)
-        if (length(pos) > 0) {
-          pos <- pos[which(jdays[pos] > 59)] 
-          jdays[pos] <- jdays[pos] + 1
-        }
-      }
-    }
-  }
-  if (length(data) != length(jdays)) {
-    stop("Parameters 'data' and 'jdays' must have the same lenght.")
-  }
-  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.")
-  }
-  climatology <- tapply(data, INDEX = jdays, FUN = mean, na.rm = na.rm)
-  anomalies <- c()
-  for (i in 1 : length(data)) {
-    index <-  which(names(climatology) == jdays[i])
-    anomalies[i] <- data[i] - climatology[index]
-  }
-  return(anomalies)
-}
\ No newline at end of file
diff --git a/R/Extremes.R b/R/Extremes.R
index 46ecbaf80d1621060a64bf6659f9bd4fb4af93f2..c565c08ad98b64ccdeecefea86174220ffccdae9 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -36,7 +36,7 @@
 #'attr(data, 'Variables')$dat1$time <- time
 #'
 #'a <- Extremes(data, threshold = threshold, op = ">", min.length = 6, spells.can.span.years = TRUE, 
-#'              max.missing.days = 5, ncores = NULL)
+#'             max.missing.days = 5, ncores = NULL)
 #'str(a)
 #'@export
 Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.years = TRUE, max.missing.days = 5, 
diff --git a/R/TimeAnomaly.R b/R/TimeAnomaly.R
new file mode 100644
index 0000000000000000000000000000000000000000..4d3e99b39f8e3fd9f1dddda4003da5428a5f883e
--- /dev/null
+++ b/R/TimeAnomaly.R
@@ -0,0 +1,133 @@
+#'Time anomalies
+#'
+#'@description This function computes temporal anomalies from a vector containing time series.
+#'@param data A numeric n-dimensional array of data including one time dimension.
+#'@param dates A vector of dates if \code{jdays} is not supplied. If NULL (by default), the 'time' attributes of parameter 'data' are considered.
+#'@param storefreq The time frequency of input data. Can be 'daily' (default) or 'monthly'.
+#'@param jdays A vector of the corresponding day of the year. This vector must be the same length as the time dimension of the parameter \code{data}.
+#'@param timedim An integer number indicating the position of the time dimension in the parameter \code{data}. If NULL (by default), use the dimension called 'time' or 'dates' in parameter \code{data}.
+#'@param na.rm A logical indicating whether missing values should be removed. If \code{na.rm} is FALSE an NA value in any of the arguments will cause a value of NA to be returned, otherwise (TRUE by default) NA values are ignored.
+#'@param ncores The number of cores to be used when computing the index.
+#'@return A named array of time anomalies of the same dimension as parameter \code{data}.
+#'
+#'@import multiApply 
+#'@import PCICt
+#'@examples
+#'# Two-year data array:
+#'data <- rnorm(2 * 3 * 731 * 2, 2, 1)
+#'dim(data) <- c(lon = 2, lat = 3, time = 365 + 366, model = 2)
+#'time <- seq(ISOdate(1903,1,1), ISOdate(1904,12,31), "days")
+#'metadata <- list(time = list(standard_name = 'time', long_name = 'time',
+#'                             calendar = 'standard',
+#'                             units = 'days since 1970-01-01 00:00:00', prec = 'double',
+#'                             dim = list(list(name ='time', unlim = FALSE))))
+#'attr(time, 'variables') <- metadata
+#'attr(data, 'Variables')$dat1$time <- time
+#'attr(data, 'Variables')$dat2$time <- time
+#'attr(data, 'Variables')$common[[2]]$dim[[3]]$len <- length(time)
+#'attr(data, 'Variables')$common[[2]]$dim[[3]]$vals <- time
+#'
+#'a <- TimeAnomaly(data, timedim = 3, na.rm = TRUE, ncores = NULL)
+#'str(a)
+#'@export
+TimeAnomaly <- function(data, dates = NULL, storefreq = 'daily', jdays = NULL, 
+                        timedim = NULL, na.rm = TRUE, ncores = NULL) {
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
+  }
+  if (is.null(dim(data))) {
+    dim(data) <- c(time = length(data))
+    timedim <- 1
+    warning(paste0("Parameter 'data' is a vector. Used as a time series ",
+                   "and parameter 'timedim' is set to be 1.")
+    )
+  }
+ if (is.null(jdays)) {
+   if (is.null(dates)) {
+     dates <- attr(data, 'Variables')$common$time
+     if (is.null(dates)) {
+       dates <- attr(data, 'Variables')$dat1$time
+     }
+   }
+     if (is.null(dates)) {
+       stop(paste0("No dates provided in parameter 'dates' nor as attribute ",
+                   "of parameter 'data'. At least one of the parameters ",
+                   "'jdays' or 'dates' must be supplied.")
+       )
+     }
+   
+    if (!any(class(dates) %in% c('POSIXct'))) {
+      dates <- try( {
+        if (is.character(dates)) {
+          as.POSIXct(dates, format = "%Y%m%d")
+        } else {
+          as.POSIXct(dates)
+        }
+      })
+      if ('try-error' %in% class(dates) | sum(is.na(dates)) == length(dates)) {
+        stop(paste0("Dates provided in parameter 'dates' must be of class ",
+                    "'POSIXct' or convertable to 'POSIXct'.")
+        )
+      }
+    }
+
+   if (storefreq == 'daily') {
+      jdays <- as.numeric(strftime(dates, format = "%j"))
+   } else if (storefreq == 'monthly') {
+      jdays <- as.numeric(strftime(dates, format = "%m"))
+   } else stop("Parameter 'storefreq' should be either 'daily' or 'monthly'.")
+
+
+    year <- as.numeric(strftime(dates, format = "%Y"))
+    if (length(unique(year)) > 1) {
+        pos <- ((year / 100) %% 1 == 0) + ((year / 4) %% 1 == 0) +
+                ((year / 400) %% 1 == 0)
+        pos <- which(pos == 0 | pos == 2 | pos == 4)
+      if (length(pos) > 0) {
+          pos <- pos[which(jdays[pos] > 59)] 
+          jdays[pos] <- jdays[pos] + 1
+      }
+    }
+ }
+
+  if (is.null(timedim)) {
+    timedim <- which(names(dim(data)) == c("dates") | 
+                     names(dim(data)) == c("time"))
+  }
+  if (length(timedim) == 0) {
+    stop(paste0("The time dimension is not found in data. Please specify ",
+                "with parameter 'timedim'.")
+    )
+  }
+
+  if (!is.logical(na.rm)) {
+    stop("Parameter 'na.rm' must be logical.")
+  }
+  if (length(na.rm) > 1) {
+    na.rm = na.rm[1]
+    warning(paste0("Parameter 'na.rm' has length > 1 and only the first ",
+                   "element will be used.")
+    )
+  }
+
+ margin <- c(1 : length(dim(data)))[-c(timedim)]
+
+ TimeAnomaly <- Apply(data, margins = margin, fun = .TimeAnomaly, jdays = jdays, 
+                      na.rm = na.rm, ncores = ncores)
+ TimeAnomaly <- TimeAnomaly$output1
+ names(dim(TimeAnomaly))[1] <- names(dim(data))[timedim]
+ if(is.null(names(dim(TimeAnomaly))[1])) {
+   names(dim(TimeAnomaly))[1] <- 'time'
+ }
+ return(TimeAnomaly)
+}
+
+.TimeAnomaly <- function(data, jdays, na.rm = na.rm) {
+  climatology <- tapply(data, INDEX = jdays, FUN = mean, na.rm = na.rm)
+  anomalies <- c()
+  for (i in 1 : length(data)) {
+    index <-  which(names(climatology) == jdays[i])
+    anomalies[i] <- data[i] - climatology[index]
+  }
+  return(anomalies)
+}
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
deleted file mode 100644
index dcd53c37459889b5610d15c5524ac93e2345a6fe..0000000000000000000000000000000000000000
--- a/man/DailyAno.Rd
+++ /dev/null
@@ -1,33 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/DailyAno.R
-\name{DailyAno}
-\alias{DailyAno}
-\title{Daily anomalies}
-\usage{
-DailyAno(data, jdays = NULL, dates = NULL, calendar = NULL,
-  na.rm = TRUE)
-}
-\arguments{
-\item{data}{A vector of daily data.}
-
-\item{jdays}{A vector of the corresponding day of the year. This vector must be the same length as parameter \code{data}.}
-
-\item{dates}{If \code{jdays} is not supplied, a vector of dates corresponding to the observations in \code{data} with defined calendar attributes.}
-
-\item{calendar}{A character indicating the calendar type.}
-
-\item{na.rm}{A logical indicating whether missing values should be removed. If \code{na.rm} is FALSE an NA value in any of the arguments will cause a value of NA to be returned, otherwise (TRUE by default) NA values are ignored.}
-}
-\value{
-A vector of daily anomalies of the same length as parameter \code{data}.
-}
-\description{
-This function computes daily anomalies from a vector containing the daily time series.
-}
-\examples{
-# Time series in a vector example:
-data <- 1:10
-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..fbad43c1151813a8d4d347d96c98189fc1bae42b 100644
--- a/man/Extremes.Rd
+++ b/man/Extremes.Rd
@@ -55,6 +55,7 @@ attr(time, "variables") <- metadata
 attr(data, 'Variables')$dat1$time <- time
 
 a <- Extremes(data, threshold = threshold, op = ">", min.length = 6, spells.can.span.years = TRUE, 
-             max.missing.days = 5, ncores = NULL)
+            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..bd7c05a9dbe1770ff6043e9fcab6c9f44d981b04 100644
--- a/man/Threshold.Rd
+++ b/man/Threshold.Rd
@@ -40,3 +40,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/TimeAnomaly.Rd b/man/TimeAnomaly.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..b3e3827f6d05902e210f480086130cc49f27fb77
--- /dev/null
+++ b/man/TimeAnomaly.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/TimeAnomaly.R
+\name{TimeAnomaly}
+\alias{TimeAnomaly}
+\title{Time anomalies}
+\usage{
+TimeAnomaly(data, dates = NULL, storefreq = "daily", jdays = NULL,
+  timedim = NULL, na.rm = TRUE, ncores = NULL)
+}
+\arguments{
+\item{data}{A numeric n-dimensional array of data including one time dimension.}
+
+\item{dates}{A vector of dates if \code{jdays} is not supplied. If NULL (by default), the 'time' attributes of parameter 'data' are considered.}
+
+\item{storefreq}{The time frequency of input data. Can be 'daily' (default) or 'monthly'.}
+
+\item{jdays}{A vector of the corresponding day of the year. This vector must be the same length as the time dimension of the parameter \code{data}.}
+
+\item{timedim}{An integer number indicating the position of the time dimension in the parameter \code{data}. If NULL (by default), use the dimension called 'time' or 'dates' in parameter \code{data}.}
+
+\item{na.rm}{A logical indicating whether missing values should be removed. If \code{na.rm} is FALSE an NA value in any of the arguments will cause a value of NA to be returned, otherwise (TRUE by default) NA values are ignored.}
+
+\item{ncores}{The number of cores to be used when computing the index.}
+}
+\value{
+A named array of time anomalies of the same dimension as parameter \code{data}.
+}
+\description{
+This function computes temporal anomalies from a vector containing time series.
+}
+\examples{
+# Two-year data array:
+data <- rnorm(2 * 3 * 731 * 2, 2, 1)
+dim(data) <- c(lon = 2, lat = 3, time = 365 + 366, model = 2)
+time <- seq(ISOdate(1903,1,1), ISOdate(1904,12,31), "days")
+metadata <- list(time = list(standard_name = 'time', long_name = 'time',
+                            calendar = 'standard',
+                            units = 'days since 1970-01-01 00:00:00', prec = 'double',
+                            dim = list(list(name ='time', unlim = FALSE))))
+attr(time, 'variables') <- metadata
+attr(data, 'Variables')$dat1$time <- time
+attr(data, 'Variables')$dat2$time <- time
+attr(data, 'Variables')$common[[2]]$dim[[3]]$len <- length(time)
+attr(data, 'Variables')$common[[2]]$dim[[3]]$vals <- time
+
+a <- TimeAnomaly(data, timedim = 3, na.rm = TRUE, 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)
 }
+
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000000000000000000000000000000000000..1f16a317626155bd0f65e8a398805c97612e5bfc
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,4 @@
+library(testthat)
+library(ClimProjDiags)
+
+test_check("ClimProjDiags")
diff --git a/tests/testthat/test-TimeAnomaly.R b/tests/testthat/test-TimeAnomaly.R
new file mode 100644
index 0000000000000000000000000000000000000000..0247ac43ceb3d56014a3b744f2b96cb183d0f010
--- /dev/null
+++ b/tests/testthat/test-TimeAnomaly.R
@@ -0,0 +1,157 @@
+context("Generic tests")
+test_that("Sanity checks", {
+
+  expect_error(
+    TimeAnomaly(data = NULL),
+    "Parameter 'data' cannot be NULL."
+  )
+
+  data <- 1:10
+  dim(data) <- c(2, 5)
+  expect_error(  
+    TimeAnomaly(data),
+    paste0("No dates provided in parameter 'dates' nor as attribute of ",
+           "parameter 'data'. At least one of the parameters 'jdays' or ",
+           "'dates' must be supplied.")
+  )
+
+  x <- c("1jan1960", "2jan1960", "31mar1960", "30jul1960")
+  expect_error(
+    TimeAnomaly(data = 1:8, dates = c(x, x)),
+    paste0("Dates provided in parameter 'dates' must be of class 'POSIXct' ",
+           "or convertable to 'POSIXct'.")
+  )
+
+  data <- 1 : 100
+  dim(data) <- c(day = 50, model = 2)
+  dates <- as.PCICt(paste0(rep(1, 50), "-", rep(c(1, 2), 25), "-",
+                    sort(rep(1901 : 1925, 2))), cal = "standard",
+                    format = "%d-%m-%Y")
+  expect_error(
+    TimeAnomaly(data = data, dates = dates, storefreq = 'seasonal'),
+    "Parameter 'storefreq' should be either 'daily' or 'monthly'."
+  )
+  expect_error(
+    TimeAnomaly(data = data, dates = dates),
+    paste0("The time dimension is not found in data. Please specify with ",
+           "parameter 'timedim'.")
+  )
+  expect_error(
+    TimeAnomaly(data = data, dates = dates, timedim = 1, na.rm = 1),
+    "Parameter 'na.rm' must be logical."
+  )
+
+
+  expect_warning(
+    TimeAnomaly(data = 1:10, jdays = rep(1:2,5)),
+    paste0("Parameter 'data' is a vector. Used as a time series and parameter ",
+           "'timedim' is set to be 1.")
+  )
+
+  expect_warning(
+    TimeAnomaly(data = data, dates = dates, timedim = 1, na.rm = c(TRUE, TRUE)),
+    "Parameter 'na.rm' has length > 1 and only the first element will be used."
+  )
+
+  tmp <- as.array(c(-4, -4, -2, -2, 0, 0, 2, 2, 4, 4))
+  names(dim(tmp)) <- c('time')
+  expect_equal(
+    TimeAnomaly(data = 1:10, jdays = rep(1:2, 5)), 
+    tmp
+  )
+
+  data <- 1 : 100
+  dim(data) <- c(time = length(data))
+  jj <- rep(1 : 2, 50)
+  output <- rep(-49 : 49, each = 2)
+  output <- output[which(output %% 2 != 0)]
+  dim(output) <- c(time = length(output))
+  expect_equal(
+    TimeAnomaly(data = data, jdays = jj), 
+    output
+  )
+
+  data <- rnorm(730 * 3 * 2 * 1, 2, 1)
+  dim(data) <- c(time = 730, memb = 3, lon = 2, lat = 1)
+  time <- seq.Date(as.Date("1990-01-01", format = "%Y-%d-%m"),
+          as.Date("1991-31-12", format = "%Y-%d-%m"), 1)
+  expect_equal(
+    dim(TimeAnomaly(data, dates = time, timedim = 1, na.rm = TRUE, ncores = NULL)),
+    dim(data)
+  )
+  data[1,1,1,1] <- NA
+  expect_equal(
+    which(is.na(TimeAnomaly(data, dates = time, timedim = 1, na.rm = TRUE,
+          ncores = NULL))),
+    c(1)
+  )
+  expect_equal(
+    which(is.na(TimeAnomaly(data, dates = time, timedim = 1, na.rm = FALSE,
+          ncores = NULL))),
+    c(1, 366)
+  )
+
+  data <- rnorm(800 * 3 * 6 * 2, 25, 8)
+  dim(data) <- c(dates = 800, lat = 3, lon = 6, model = 2)
+  dates <- seq(ISOdate(2000, 1, 1), ISOdate(2002, 03, 10), "days")
+  metadata <- list(time = list(standard_name = 'time', long_name = 'time',
+                               calendar = 'standard',
+                               units = 'days since 1970-01-01 00:00:00', 
+                               prec = 'double',
+                               dim = list(list(name = 'time', unlim = FALSE))))
+  attr(dates, 'variables') <- metadata
+  attr(data, 'Variables')$dat1$time <- dates
+  attr(data, 'Variables')$dat2$time <- dates
+  attr(data, 'Variables')$common[[2]]$dim[[3]]$len = length(dates)
+  attr(data, 'Variables')$common[[2]]$dim[[3]]$vals <- dates
+  expect_equal(
+    dim(TimeAnomaly(data)),
+    c(dates = 800, lat = 3, lon = 6, model = 2)
+  )
+
+  data <- 1 : 24
+  dim(data) <- c(time = 24, lat = 1, lon = 1)
+  dates <- seq(ISOdate(2000, 1, 1), ISOdate(2001, 12, 31), "months")
+  metadata <- list(time = list(standard_name = 'time', long_name = 'time',
+                               calendar = 'standard',
+                               units = 'days since 1970-01-01 00:00:00',
+                               prec = 'double',
+                               dim = list(list(name ='time', unlim = FALSE))))
+  attr(dates, 'variables') <- metadata
+  attr(data, 'Variables')$dat1$time <- dates
+  attr(data, 'Variables')$dat2$time <- dates
+  attr(data, 'Variables')$common[[2]]$dim[[3]]$len = length(dates)
+  attr(data, 'Variables')$common[[2]]$dim[[3]]$vals <- dates
+  tmp <- sort(rep(c(6, -6), 12))
+  dim(tmp) <- c(time = 24, lat = 1, lon = 1)
+  expect_equal(
+    TimeAnomaly(data, storefreq = 'monthly'),
+    tmp
+  )
+
+  data <- 1:(1 * 1 * 731 * 1)
+  dim(data) <- c(lon = 1, lat = 1, time = 365 + 366, model = 1)
+  time <- seq(ISOdate(1903,1,1), ISOdate(1904,12,31), "days")
+  metadata <- list(time = list(standard_name = 'time', long_name = 'time',
+                               calendar = 'standard',
+                               units = 'days since 1970-01-01 00:00:00',
+                               prec = 'double',
+                               dim = list(list(name = 'time', unlim = FALSE))))
+  attr(time, 'variables') <- metadata
+  attr(data, 'Variables')$dat1$time <- time
+  attr(data, 'Variables')$dat2$time <- time
+  attr(data, 'Variables')$common[[2]]$dim[[3]]$len = length(time)
+  attr(data, 'Variables')$common[[2]]$dim[[3]]$vals <- time
+  aa <- rep(-182.5, 59)
+  bb <- rep(-183, 365-59)
+  cc <- rep(182.5, 59)
+  dd <- 0
+  ee <- rep(183, 366-60)
+  tmp <- c(aa, bb, cc, dd, ee)
+  dim(tmp) <- c(time = 731, lon = 1, lat = 1, model = 1)
+  expect_equal(
+    TimeAnomaly(data),
+    tmp
+  )
+
+})
diff --git a/vignettes/extreme_indices.Rmd b/vignettes/extreme_indices.Rmd
index 57eb34f29c044c8a390614e50bc43ca05bfa0d39..4ec5a6d4e10c01530c9dd6a445f21ac77b82156a 100644
--- a/vignettes/extreme_indices.Rmd
+++ b/vignettes/extreme_indices.Rmd
@@ -149,13 +149,11 @@ In order to evaluate the future projections, it is necessary to compute the inde
 The next steps should be followed:
 
 
-To remove seasonality effects, the anomaly is computed for each day and gridpoint by applying the `DailyAno` function. The name of the first dimensions is defined as 'time' dimension. 
+To remove seasonality effects, the anomaly is computed for each day and gridpoint by applying the `TimeAnomaly` function. 
 
 
 ```r
-anomaly_data <- apply(tmax_historical, c(1,2,4,5), DailyAno, dates = attributes(tmax_historical)$Variables$dat1$time)
-
-names(dim(anomaly_data))[1] <- "time"
+anomaly_data <- TimeAnomaly(tmax_historical)
 ```