diff --git a/DESCRIPTION b/DESCRIPTION
index 78935b854eb76f5ab7dce4fe18216a6c0588ccde..9bd337c40d5957b525a501d1d90fbc03374618b5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -50,6 +50,7 @@ Imports:
     ncdf4,
     parallel,
     plyr,
+    multiApply,
     SpecsVerification (>= 0.5.0)
 Suggests:
     easyVerification,
@@ -61,3 +62,4 @@ LazyData: true
 SystemRequirements: cdo
 Encoding: UTF-8
 RoxygenNote: 5.0.0
+
diff --git a/NAMESPACE b/NAMESPACE
index 3e35710e2d54618b5a2311a1a03f839de76e59f3..dc29243de6bad6473679214674ff935cae723998 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -89,6 +89,7 @@ import(graphics)
 import(mapproj)
 import(maps)
 import(methods)
+import(multiApply)
 import(ncdf4)
 import(parallel)
 import(plyr)
diff --git a/R/Season.R b/R/Season.R
index d9f95a51cd33a77d4a2f29079a21ca3a5d2f49d3..16cfa1639081a153d515a050e35ae81b3813e454 100644
--- a/R/Season.R
+++ b/R/Season.R
@@ -1,96 +1,203 @@
 #'Computes Seasonal Means
 #'
-#'Computes seasonal means on timeseries organized in a array of any number of 
-#'dimensions up to 10 dimensions where the time dimension is one of those 10 
-#'dimensions.
+#'Computes seasonal means (or other operations) on monthly timeseries from n-dimensional arrays with named dimensions
 #'
-#'@param var Array containing the timeseries along one of its dimensions. 
-#'@param posdim Dimension along which to compute seasonal means = Time 
-#'  dimension.
-#'@param monini an integer indicating the first month of the time series: 1 to 
-#'  12.
-#'@param moninf an integer indicating the month when to start the seasonal 
-#'  means: 1 to 12.
-#'@param monsup an integer indicating the month when to stop the seasonal 
-#'  means: 1 to 12.
+#'@param var a numeric n-dimensional array with named dimensions on monthy frequency. 
+#'@param posdim a character indicating the name of the dimension or a integer numeric indicating the position along which to compute seasonal means (or other operations). By default, 'time' dimension is expected.
+#'@param monini an integer indicating the first month of the time series: 1 to 12.
+#'@param moninf an integer indicating the month when to start the seasonal means: 1 to 12.
+#'@param monsup an integer indicating the month when to stop the seasonal means: 1 to 12.
+#'@param operation a character or function indicating the name of a function to be applied in seasonal basins. By default, means are computed. Other operations can be 'sum' for total precipitation.
+#'@param na.rm a logical value indicating whether NA values should be stripped before the computation proceeds.
 #'
-#'@return Array with the same dimensions as var except along the posdim 
-#'  dimension whose length corresponds to the number of seasons. Partial 
-#'  seasons are not accounted for.
-#'@keywords datagen
-#'@author History:\cr
-#'0.1  -  2011-03  (V. Guemas, \email{virginie.guemas@ic3.cat})  -  Original code\cr
-#'1.0  -  2013-09  (N. Manubens, \email{nicolau.manubens@ic3.cat})  -  Formatting to CRAN
+#'@return Array with the same dimensions as var except along the posdim dimension whose length corresponds to the number of seasons. Partial seasons are not accounted for.
+#'
+#'@import multiApply
 #'@examples
-#'# Load sample data as in Load() example:
-#'example(Load)
-#'leadtimes_dimension <- 4
-#'initial_month <- 11
-#'mean_start_month <- 12
-#'mean_stop_month <- 2
-#'season_means_mod <- Season(sampleData$mod, leadtimes_dimension, initial_month,
-#'                           mean_start_month, mean_stop_month)
-#'season_means_obs <- Season(sampleData$obs, leadtimes_dimension, initial_month,
-#'                           mean_start_month, mean_stop_month)
-#'PlotAno(season_means_mod, season_means_obs, startDates, 
-#'        toptitle = paste('winter (DJF) temperatures'), ytitle = c('K'), 
-#'        legends = 'ERSST', biglab = FALSE, fileout = 'tos_season_means.eps')
+#'dat <- 1 : (12 * 5 * 2 * 3 * 2)
+#'dim(dat) <- c(dat = 1, memb = 3, time = 12 * 5, lon = 2, lat = 2) 
+#'res <- Season(var = dat, monini = 1, moninf = 1, monsup = 2)
+#'dat <- 1 : (12 * 2 * 3)
+#'dim(dat) <- c(dat = 2, time = 12, memb = 3)
+#'res <- Season(var = dat, monini = 1, moninf = 1, monsup = 2)
+#'dat <- 1 : (24 * 2 * 3)
+#'dim(dat) <- c(dat = 2, time = 24, memb = 3)
+#'res <- Season(var = dat, monini = 1, moninf = 1, monsup = 2)
+#'dat[c(1, 20, 50)] <- NA
+#'res1 <- Season(var = dat, monini = 1, moninf = 1, monsup = 2, na.rm = FALSE)
+#'res2 <- Season(var = dat, monini = 1, moninf = 1, monsup = 2, na.rm = TRUE)
+#'@import multiApply
 #'@export
-Season <- function(var, posdim = 4, monini, moninf, monsup) {
-  while (monsup < moninf) {
-    monsup <- monsup + 12
-  }
-  #
-  #  Enlarge the size of var to 10 
-  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  # 
-  dimsvar <- dim(var)
-  if (is.null(dimsvar)) {
-    dimsvar <- length(var)
-  }
-  ntime <- dimsvar[posdim]
-  enlvar <- Enlarge(var, 10)
-  outdim <- c(dimsvar, array(1, dim = (10 - length(dimsvar))))
-  u <- IniListDims(outdim, 10)
-  v <- IniListDims(outdim, 10)
-  #
-  #  Compute the seasonal means 
-  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  #
-  ind <- 1:ntime
-  months <- ((ind - 1) + monini - 1) %% 12 + 1
-  years <- ((ind - 1) + monini - 1) %/% 12
+Season <- function(var, posdim = 'time', monini, moninf, monsup, 
+                   operation = mean, na.rm = FALSE) {
+    # Check var
+    if (is.null(var)) {
+        stop("Parameter 'var' cannot be NULL.")
+    }
+    if (!is.numeric(var)) {
+        stop("Parameter 'var' must be a numeric array.")
+    }
+    if (is.null(dim(var))) {
+        dim(var) <- c(length(var))
+        if (is.character(posdim)) {
+            names(dim(var)) <- posdim
+        } else {
+            names(dim(var)) <- 'time'
+            posdim <- 'time'
+        }
+        time_dim <- 1
+    } else {
+        if (is.null(names(dim(var)))) {
+             if (is.numeric(posdim)) {
+                names(dim(var)) <- paste0("D", 1 : length(dim(var)))  
+                names(dim(var))[posdim] <- 'time'
+                time_dim <- posdim
+                posdim <- 'time'
+            } else  {
+                stop("Parameter 'var' must contain dimesnion names.")
+            }
+        } else {
+             if (is.numeric(posdim)) {
+                time_dim <- posdim 
+                posdim <- names(dim(var))[posdim]
+             } else if (is.character(posdim)) {
+                 time_dim <- which(names(dim(var)) == posdim)
+             } else {
+                stop("Parameter 'posdim' must be a integer or character",
+                     "indicating the 'time' dimension.")
+            }
+        }
+    }
+
+    dim_names <- names(dim(var))
+#    series <- apply(var, margins, .Season, 
+#                    monini = monini, moninf = moninf, monsup = monsup,
+#                    operation = operation)
+    series <- Apply(list(var), 
+                    target_dims = posdim, 
+                    fun = .Season, 
+                    monini = monini, moninf = moninf, monsup = monsup,
+                    operation = operation, na.rm = na.rm)$output1
+
+    if (is.null(dim(series))) {
+        dim(series) <- c(time = length(series))
+
+    } else {  #is an array
+      if (length(dim(series)) == 1) {
+        if (!is.null(names(dim(series)))) { 
+          dim(series) <- c(1, dim(series))
+          names(dim(series))[1] <- posdim
+        } else {
+          names(dim(series)) <- posdim
+        }
+      } else {  # more than one dim
+        if (names(dim(series))[1] != "") {
+          dim(series) <- c(1, dim(series))
+          names(dim(series))[1] <- posdim
+        } else {   #regular case
+          names(dim(series))[1] <- posdim
+        }
+      }
+    }
+
+    if (any(dim_names != names(dim(series)))) {
+        pos <- match(dim_names, names(dim(series)))
+        series <- aperm(series, pos)
+        names(dim(series)) <- dim_names
+    }
+    return(series)
+}
 
-  for (jmon in moninf:monsup) {
-    u[[posdim]] <- ind[which(months == ((jmon - 1) %% 12 + 1))]
-    ind0 <- u[[posdim]][1]
-    indf <- u[[posdim]][length(u[[posdim]])]
-    if (indf > (ntime - (monsup - jmon))) {
-      u[[posdim]] <- u[[posdim]][-which(u[[posdim]] == indf)]
-    }
-    if (ind0 < (jmon - moninf + 1)) {
-      u[[posdim]] <- u[[posdim]][-which(u[[posdim]] == ind0)]
+.Season <- function(x, monini, moninf, monsup, operation = mean, na.rm = FALSE) {
+    # Checks:
+    if (!is.numeric(x)) {
+        stop("Parameter 'x' must be a numeric vector.")
+    }
+    if (!is.numeric(monini)) {
+        stop("Parameter 'monini' must be numeric.")
     } 
-    if (jmon == moninf) { 
-      nseas <- length(u[[posdim]])
-      dimsvar[posdim] <- nseas
-      outdim[posdim] <- nseas
-      enlvarout <- array(0, dim = outdim)
-    }
-    v[[posdim]] <- 1:nseas
-    enlvarout[v[[1]], v[[2]], v[[3]], v[[4]], v[[5]], v[[6]], v[[7]], v[[8]], 
-              v[[9]], v[[10]]] <- enlvarout[v[[1]], v[[2]], v[[3]], v[[4]],
-                                            v[[5]], v[[6]], v[[7]], v[[8]], 
-                                            v[[9]], v[[10]]] + enlvar[u[[1]], 
-                                  u[[2]], u[[3]], u[[4]], u[[5]], u[[6]], 
-                                  u[[7]], u[[8]], u[[9]], u[[10]]]
-  }
-  varout <- array(dim = dimsvar)
-  varout[] <- enlvarout
-  varout <- varout / (monsup - moninf + 1)
-  #
-  #  Outputs
-  # ~~~~~~~~~
-  #
-  varout
+    if (length(monini) > 1) {
+        monini <- monini[1] 
+        warning("Parameter 'monini' has length > 1 and only the first ",
+                "element will be used.")
+    }
+    if (monini %% 1 != 0) {
+        stop("Parameter 'monini' must be an integer.")
+    }
+    if (!is.numeric(moninf)) {
+        stop("Parameter 'moninf' must be numeric.")
+    }
+    if (length(moninf) > 1) {
+        moninf <- moninf[1]
+        warning("Parameter 'moninf' has length > 1 and only the first ",
+                "element will be used.")
+    }
+    if (moninf %% 1 != 0) {
+        stop("Parameter 'moninf' must be an integer.")
+    }
+    if (!is.numeric(monsup)) {
+        stop("Parameter 'monsup' must be numeric.")
+    }
+    if (length(monsup) > 1) {
+        monsup <- monsup[1]
+        warning("Parameter 'monsup' has length > 1 and only the first ",
+                "element will be used.")
+    }
+    if (monsup %% 1 != 0) {
+        stop("Parameter 'monsup' must be an integer.")
+    }
+    # Check fun operation
+    if (is.character(operation)) {
+        fun_name <- operation
+        err <- try({operation <- get(operation)}, silent = TRUE)
+        if (!is.function(operation)) {
+            stop("Could not find the function '", fun_name, "'.")
+        }
+     }
+     if (!is.function(operation)) {
+         stop("Parameter 'operation' must be a function or a character string ",
+              "with the name of a function.")
+     }
+
+#    # Correction e.g. 'winter':
+#    while (monsup < moninf) {
+#        monsup <- monsup + 12
+#    }
+    # Correction need if monini is not January: 
+    moninf <- moninf - monini + 1
+    monsup <- monsup - monini + 1
+    moninf <- ifelse(moninf <= 0, moninf + 12, moninf)
+    monsup <- ifelse(monsup <= 0, monsup + 12, monsup)
+
+    while (monsup < moninf) {
+        monsup <- monsup + 12
+    }
+
+    #### Create position index:
+    # Basic index:
+    pos <- moninf : monsup
+    # Extended index for all period:
+    if (length(x) > pos[length(pos)]) {
+        pos2 <- lapply(pos, function(y) {seq(y, length(x), 12)})
+    } else {
+        pos2 <- pos
+    }
+    # Correct if the final season is not complete:
+    maxyear <- min(unlist(lapply(pos2, length)))
+    pos2 <- lapply(pos2, function(y) {y[1 : maxyear]})
+    # Convert to array:
+    pos2 <- unlist(pos2)
+    dim(pos2) <- c(year = maxyear, month = length(pos2)/maxyear)
+    if (na.rm == TRUE) {
+        if (length(which(is.na(x[pos2]))) > 0) {
+             pos2 <- apply(pos2, 1, function(y){y[!is.na(x[y])]})
+             timeseries <- unlist(lapply(pos2, function(y) {operation(x[y])}))
+        } else {
+             timeseries <- apply(pos2, 1, function(y) {operation(x[y])})
+        }
+    } else {
+        timeseries <- apply(pos2, 1, function(y) {operation(x[y])})
+    }
+
+    return(timeseries)
 }
diff --git a/R/Trend.R b/R/Trend.R
index 98a71bf35221fc35b3913331b943cede0f3e386c..d8bfc6014b4be5d19c6224d967f0b18b52e79cad 100644
--- a/R/Trend.R
+++ b/R/Trend.R
@@ -8,26 +8,16 @@
 #'.Trend provides the same functionality but taking a matrix ensemble members 
 #'as input (exp).
 #'
-#'@param var Array of any number of dimensions up to 10.
-#'@param exp M by N matrix of M forecasts from N ensemble members.
-#'@param interval Number of months/years between 2 points along posTR 
-#'  dimension.\cr
-#'  Default = 1.\cr
-#'  The trend would be provided in number of units per month or year.
-#'@param siglev Confidence level for the computation of confidence interval. 
-#'  0.95 by default.
-#'@param conf Whether to compute the confidence levels or not. TRUE by default.
-#'@param posTR Position along which to compute the trend.
 #'@param var An array of any number of dimensions up to 10.
-#'@param exp An M by N matrix representing M forecasts from N ensemble members.
+#'@param posTR An integer indicating the position along which to compute the 
+#'  trend.
 #'@param interval A number of months/years between 2 points along posTR 
 #'  dimension. Set 1 as default.
 #'@param siglev A numeric value indicating the confidence level for the 
 #'  computation of confidence interval. Set 0.95 as default.
 #'@param conf A logical value indicating whether to compute the confidence 
 #'  levels or not. Set TRUE as default.
-#'@param posTR An integer indicating the position along which to compute the 
-#'  trend.
+#'@param exp An M by N matrix representing M forecasts from N ensemble members.
 #'
 #'@return 
 #'\item{$trend}{
diff --git a/man/Season.Rd b/man/Season.Rd
index 08487256eea4cb2072078adbc9cdaf1e7b1158fd..b141883fe18b30b91e9b13c877ad9dfe52fd1fb8 100644
--- a/man/Season.Rd
+++ b/man/Season.Rd
@@ -4,52 +4,42 @@
 \alias{Season}
 \title{Computes Seasonal Means}
 \usage{
-Season(var, posdim = 4, monini, moninf, monsup)
+Season(var, posdim = "time", monini, moninf, monsup, operation = mean,
+  na.rm = FALSE)
 }
 \arguments{
-\item{var}{Array containing the timeseries along one of its dimensions.}
+\item{var}{a numeric n-dimensional array with named dimensions on monthy frequency.}
 
-\item{posdim}{Dimension along which to compute seasonal means = Time 
-dimension.}
+\item{posdim}{a character indicating the name of the dimension or a integer numeric indicating the position along which to compute seasonal means (or other operations). By default, 'time' dimension is expected.}
 
-\item{monini}{an integer indicating the first month of the time series: 1 to 
-12.}
+\item{monini}{an integer indicating the first month of the time series: 1 to 12.}
 
-\item{moninf}{an integer indicating the month when to start the seasonal 
-means: 1 to 12.}
+\item{moninf}{an integer indicating the month when to start the seasonal means: 1 to 12.}
 
-\item{monsup}{an integer indicating the month when to stop the seasonal 
-means: 1 to 12.}
+\item{monsup}{an integer indicating the month when to stop the seasonal means: 1 to 12.}
+
+\item{operation}{a character or function indicating the name of a function to be applied in seasonal basins. By default, means are computed. Other operations can be 'sum' for total precipitation.}
+
+\item{na.rm}{a logical value indicating whether NA values should be stripped before the computation proceeds.}
 }
 \value{
-Array with the same dimensions as var except along the posdim 
- dimension whose length corresponds to the number of seasons. Partial 
- seasons are not accounted for.
+Array with the same dimensions as var except along the posdim dimension whose length corresponds to the number of seasons. Partial seasons are not accounted for.
 }
 \description{
-Computes seasonal means on timeseries organized in a array of any number of 
-dimensions up to 10 dimensions where the time dimension is one of those 10 
-dimensions.
+Computes seasonal means (or other operations) on monthly timeseries from n-dimensional arrays with named dimensions
 }
 \examples{
-# Load sample data as in Load() example:
-example(Load)
-leadtimes_dimension <- 4
-initial_month <- 11
-mean_start_month <- 12
-mean_stop_month <- 2
-season_means_mod <- Season(sampleData$mod, leadtimes_dimension, initial_month,
-                          mean_start_month, mean_stop_month)
-season_means_obs <- Season(sampleData$obs, leadtimes_dimension, initial_month,
-                          mean_start_month, mean_stop_month)
-PlotAno(season_means_mod, season_means_obs, startDates, 
-       toptitle = paste('winter (DJF) temperatures'), ytitle = c('K'), 
-       legends = 'ERSST', biglab = FALSE, fileout = 'tos_season_means.eps')
-}
-\author{
-History:\cr
-0.1  -  2011-03  (V. Guemas, \email{virginie.guemas@ic3.cat})  -  Original code\cr
-1.0  -  2013-09  (N. Manubens, \email{nicolau.manubens@ic3.cat})  -  Formatting to CRAN
+dat <- 1 : (12 * 5 * 2 * 3 * 2)
+dim(dat) <- c(dat = 1, memb = 3, time = 12 * 5, lon = 2, lat = 2) 
+res <- Season(var = dat, monini = 1, moninf = 1, monsup = 2)
+dat <- 1 : (12 * 2 * 3)
+dim(dat) <- c(dat = 2, time = 12, memb = 3)
+res <- Season(var = dat, monini = 1, moninf = 1, monsup = 2)
+dat <- 1 : (24 * 2 * 3)
+dim(dat) <- c(dat = 2, time = 24, memb = 3)
+res <- Season(var = dat, monini = 1, moninf = 1, monsup = 2)
+dat[c(1, 20, 50)] <- NA
+res1 <- Season(var = dat, monini = 1, moninf = 1, monsup = 2, na.rm = FALSE)
+res2 <- Season(var = dat, monini = 1, moninf = 1, monsup = 2, na.rm = TRUE)
 }
-\keyword{datagen}
 
diff --git a/man/Trend.Rd b/man/Trend.Rd
index 4326a1e3b8c6743020222ac9662c8f0aa0ad439d..71d2e1bed2c0b78c96116274cbcc02449b7a0cbc 100644
--- a/man/Trend.Rd
+++ b/man/Trend.Rd
@@ -10,25 +10,10 @@ Trend(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE)
 .Trend(exp, interval = 1, siglev = 0.95, conf = TRUE)
 }
 \arguments{
-\item{var}{Array of any number of dimensions up to 10.}
-
-\item{posTR}{Position along which to compute the trend.}
-
-\item{interval}{Number of months/years between 2 points along posTR 
-dimension.\cr
-Default = 1.\cr
-The trend would be provided in number of units per month or year.}
-
-\item{siglev}{Confidence level for the computation of confidence interval. 
-0.95 by default.}
-
-\item{conf}{Whether to compute the confidence levels or not. TRUE by default.}
-
-\item{exp}{M by N matrix of M forecasts from N ensemble members.}
-
 \item{var}{An array of any number of dimensions up to 10.}
 
-\item{exp}{An M by N matrix representing M forecasts from N ensemble members.}
+\item{posTR}{An integer indicating the position along which to compute the 
+trend.}
 
 \item{interval}{A number of months/years between 2 points along posTR 
 dimension. Set 1 as default.}
@@ -39,8 +24,7 @@ computation of confidence interval. Set 0.95 as default.}
 \item{conf}{A logical value indicating whether to compute the confidence 
 levels or not. Set TRUE as default.}
 
-\item{posTR}{An integer indicating the position along which to compute the 
-trend.}
+\item{exp}{An M by N matrix representing M forecasts from N ensemble members.}
 }
 \value{
 \item{$trend}{
diff --git a/tests/testthat/test-Season.R b/tests/testthat/test-Season.R
new file mode 100644
index 0000000000000000000000000000000000000000..18ce9defff529d03656b0a0fbf4d1839d86178b4
--- /dev/null
+++ b/tests/testthat/test-Season.R
@@ -0,0 +1,118 @@
+context("Generic tests")
+test_that("Sanity checks", {
+  expect_error(
+    Season(var = "A"), "Parameter 'var' must be a numeric array.") 
+  expect_error(
+    Season(var = 1, monini = "A", moninf = "B", monsup = "C"),
+    "Parameter 'monini' must be numeric.")
+  expect_error(
+    Season(var = 1, monini = 1, moninf = "B", monsup = "C"),
+    "Parameter 'moninf' must be numeric.")
+  expect_error(
+    Season(var = 1, monini = 1, moninf = 1, monsup = "C"),
+    "Parameter 'monsup' must be numeric.")
+  expect_error(
+    Season(var = 1, monini = 1.5, moninf = 1, monsup = "C"),
+    "Parameter 'monini' must be an integer.")
+  expect_error(
+    Season(var = 1, monini = 1, moninf = 1.5, monsup = "C"),
+    "Parameter 'moninf' must be an integer.")
+ expect_error(
+    Season(var = 1, monini = 1, moninf = 1, monsup = 1.5),
+    "Parameter 'monsup' must be an integer.")
+  expect_warning(
+    Season(var = 1, monini = c(1, 2), moninf = 1, monsup = 1),
+   "Parameter 'monini' has length > 1 and only the first element will be used.")
+  
+  expect_warning(
+    Season(var = 1, monini = 1, moninf = c(1, 2), monsup = 1),
+   "Parameter 'moninf' has length > 1 and only the first element will be used.")
+  
+  expect_warning(
+    Season(var = 1, monini = 1, moninf = 1, monsup = c(1, 2)),
+   "Parameter 'monsup' has length > 1 and only the first element will be used.")
+
+   expect_error(
+    Season(var = 1, monini = 1, moninf = 1, monsup = 1, operation = "RRR"),
+    "Could not find the function 'RRR'.")
+
+  expect_error(Season(var = NULL), "Parameter 'var' cannot be NULL.")
+   
+  var <- array(1 : 20, dim = c(2, 4, 5))
+  expect_error(Season(var),
+               "Parameter 'var' must contain dimesnion names.") 
+  expect_equal(Season(var, posdim = 3, monini = 1, moninf = 1, monsup = 2),
+               array(5 : 12, dim = c(D1 = 2, D2 = 4, time = 1)))
+
+  var <- array(1 : 20, dim = c(2, 4, 5))
+  names(dim(var)) <- c('x', 'y', 'time')
+  output <- array(1 : 8, dim = c(x = 2, y = 4, time = 1))
+  expect_equal(
+    Season(var, monini = 1, moninf = 1, monsup = 1), output)
+  output <- array(5 : 12, dim = c(x = 2, y = 4, time = 1))
+  expect_equal(
+    Season(var, monini = 1, moninf = 1, monsup = 2), output)
+  output <- array(rep(seq(1.5, 19.5, 2),2), dim = c(x = 1, y = 4, time = 5))
+  expect_equal(
+    Season(var, posdim = 1, monini = 1, moninf = 1, monsup = 2), output)
+
+  var <- 1 : (12 * 2 * 3)
+  dim(var) <- c(dat = 2, time = 12, sdate = 3)
+  output <- array(c(2 : 3, 26 : 27, 50 : 51), 
+                  dim = c(dat = 2, time = 1, sdate = 3))
+  expect_equal(
+    Season(var, posdim = 'time', monini = 1, moninf = 1, monsup = 2),
+    output)
+  expect_equal(
+    Season(var, posdim = 2, monini = 1, moninf = 1, monsup = 2),
+    output)
+  output <- array(13 : 36, dim = c(dat = 2, time = 12, sdate = 1))
+  expect_equal(
+    Season(var, posdim = 3, monini = 1, moninf = 1, monsup = 2),
+    output)
+
+  output <- array(c(11:12, 35:36, 59:60), dim = c(dat = 2, time = 1, sdate = 3))
+  expect_equal(
+    Season(var, posdim = 2, monini = 10, moninf = 2, monsup = 4),
+    output)
+
+  var <- 1 : 10 # caso 1
+  output <- array(1.5, dim = c(time = 1))
+  expect_equal(
+    Season(var, posdim = 'time', monini = 1, moninf = 1, monsup = 2), 
+    output)
+
+  var <- 1 : 10 ; dim(var) <- c(2, 5) #caso 2
+  expect_error(
+     Season(var, posdim = 'time', monini = 1, moninf = 1, monsup = 2),
+     "Parameter 'var' must contain dimesnion names.")
+  output <- array(2 : 3, dim = c(D1 = 2, time = 1))
+  expect_equal(
+     Season(var, posdim = 2, monini = 1, moninf = 1, monsup = 2),
+     output)
+  output <- array(c(4, 6), dim = c(D1 = 2, time = 1))
+  expect_equal(
+     Season(var, posdim = 2, monini = 1, moninf = 1, monsup = 2, operation = sum),
+     output)
+
+  var <- array(1:1200, dim = c(time = 120, lat = 5, lon = 2))
+  output <- array(dim = c(time = 9, lat = 5, lon = 2))
+  output[, 1, 1] <- seq(8.5, 104.5, 12)
+  for (i in 2:5) {
+    output[1, i, 1] <- output[9, i - 1, 1] + 24
+    output[, i, 1] <- seq(output[1, i, 1], output[1, i, 1] + 96, 12)
+  }
+  output[1, 1, 2] <- output[9, 5, 1] + 24
+  output[, 1, 2] <- seq(output[1, 1, 2], output[1, 1, 2] + 96, 12)
+  for (i in 2:5) {
+    output[1, i, 2] <- output[9, i - 1, 2] + 24
+    output[, i, 2] <- seq(output[1, i, 2], output[1, i, 2] + 96, 12)
+  }
+
+  expect_equal(
+     Season(var, posdim = 'time', monini = 11, moninf = 1, monsup = 12),
+     output)
+
+
+})
+