diff --git a/DESCRIPTION b/DESCRIPTION
index c96b4a3..4409e7c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -23,8 +23,8 @@ Suggests:
     visualizeR
 Type: Package
 Title: A climate4R package for general climate data manipulation and transformation
-Version: 2.0.2
-Date: 2020-09-02
+Version: 2.1.2
+Date: 2021-07-07
 Authors@R: c(person("Bedia", "Joaquín", email = "bediaj@unican.es", role = c("aut","cre"), comment = c(ORCID = "0000-0001-6219-4312")),
              person("Jorge", "Baño Medina", email = "bmedina@ifca.unican.es", role = "aut"),
              person("Ana", "Casanueva", email = "ana.casanueva@unican.es", role = "aut"),
diff --git a/NAMESPACE b/NAMESPACE
index 2312530..db7236c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -21,7 +21,9 @@ export(dataSplit)
 export(detrendGrid)
 export(easyVeri.signif)
 export(easyVeri2grid)
+export(fillGrid)
 export(fillGridDates)
+export(fillGridSpatial)
 export(filterGrid)
 export(filterNA)
 export(get2DmatCoordinates)
diff --git a/NEWS b/NEWS
index 73ff8b6..1da7bef 100644
--- a/NEWS
+++ b/NEWS
@@ -120,9 +120,20 @@ This release of transformeR does not ensure backwards compatibility. Some deprec
  * Force nearest-neighbour application among non-overlapping domains in `interpGrid`
  * Other minor documentation updates
 
-## 2.0.2 (in devel)
+## 2.1.0 (17 Mar 2021)
+ * Implementation of Lamb Weather Typing for southern hemisphere
+ * Inclusion of type 27 (unclassified) in Lamb WT
+ * New functionality for spatial aggregation of rectangular domains
+ * Bug fix in partially masked latitudinal weighted domain averages.
+ * New option to use masking in interpGrid
  * Bug fix in subsetVar attributes
- * Internal falg in bindGrid.member to skip temporal checks
- * Minor changes and documentation updates
+ * Internal flag in bindGrid.member to skip temporal checks
+ * Other minor changes and documentation updates
 
+## 2.1.1 (31 May 2021)
+ * New dateutils-family helper functions for consistent date-time class handling (including time-zone defs)
+ * Consistent treatment of multiple time resolutions in fillGridDates 
 
+## 2.1.2 (07 Jul 2021)
+ * Allow for non-rescaled EOF recovery in EOF2clim
+ * Ensure regular grids in upscaleGrid outputs 
\ No newline at end of file
diff --git a/R/EOF2clim.R b/R/EOF2clim.R
index bd2a1d8..3b56b89 100644
--- a/R/EOF2clim.R
+++ b/R/EOF2clim.R
@@ -2,14 +2,29 @@
 #' @description Converts a \code{prinComp} output to a climatology of EOFs. 
 #' @param prinCompObj PCA object
 #' @param ind.var index position of the target variable
-#' @param member index position of the member whose EOFs are to be displayed
-#' @param n.eofs Number of EOFs to be displayed (from 1 to n.eofs)
+#' @param member Integer. index position of the member whose EOFs are to be displayed.
+#'  Default to 1 (ignored if no members available)
+#' @param n.eofs Integer. Number of EOFs to be displayed (from 1 to \code{n.eofs}). Default to one -first- EOF.
+#' @param rescale Logical flag. Default to \code{TRUE}, which rescales the EOF to the original input variable units.
 #' @return A climatological grid. Note that EOFs are treated as members
 #' @importFrom magrittr %>% 
 #' @export
 #' @author J Bedia
+#' @examples \dontrun{
+#' require(climate4R.datasets)
+#' data("NCEP_Iberia_hus850", "NCEP_Iberia_psl", "NCEP_Iberia_ta850")
+#' multigrid <- makeMultiGrid(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
+#' # In this example, we retain the PCs explaining the 99\% of the variance
+#' pca <- prinComp(multigrid, v.exp = c(.95,0.90,.90), keep.orig = FALSE)
+#' require(visualizeR)
+#' # Original EOF
+#' EOF2clim(pca, ind.var = "psl") %>% spatialPlot(backdrop.theme = "coastline")
+#' # Rescaled EOF
+#' EOF2clim(pca, ind.var = "psl", rescale = FALSE) %>% spatialPlot(backdrop.theme = "coastline")
+#' }
 
-EOF2clim <- function(prinCompObj, ind.var, member, n.eofs) {
+
+EOF2clim <- function(prinCompObj, ind.var, member = 1L, n.eofs = 1L, rescale = TRUE) {
   varNames <- attributes(prinCompObj)$names
   levs <- attributes(prinCompObj)$level
   x <- attributes(prinCompObj)$xCoords
@@ -22,7 +37,10 @@ EOF2clim <- function(prinCompObj, ind.var, member, n.eofs) {
   eofs <- prinCompObj[[ind.var]][[member]]$EOFs[, 1:n.eofs, drop = FALSE]
   prinCompObj <- NULL
   # Rescale EOFs
-  aux <- (eofs * sigma + mu) %>% t() %>% mat2Dto3Darray(x, y) %>% list()
+  if (isTRUE(rescale)) {
+    eofs <- (eofs * sigma + mu)
+  }
+  aux <- t(eofs) %>% mat2Dto3Darray(x, y) %>% list()
   # Recover grids structure (EOFS are treated as members, time = 1 like a climatology
   Data <- do.call("abind", c(aux, along = -1)) %>% unname()
   attr(Data, "dimensions") <- c("time", "member", "lat", "lon")
diff --git a/R/aggregateGrid.R b/R/aggregateGrid.R
index a7189dd..e2d727b 100644
--- a/R/aggregateGrid.R
+++ b/R/aggregateGrid.R
@@ -1,6 +1,6 @@
 #     aggregateGrid.R Flexible grid aggregation along selected dimensions
 #
-#     Copyright (C) 2017 Santander Meteorology Group (http://www.meteo.unican.es)
+#     Copyright (C) 2021 Santander Meteorology Group (http://www.meteo.unican.es)
 #
 #     This program is free software: you can redistribute it and/or modify
 #     it under the terms of the GNU General Public License as published by
@@ -20,16 +20,19 @@
 #' @description Aggregates a grid along the target dimensions using user-defined functions.
 #' @param grid a grid or multigrid to be aggregated.
 #' @param aggr.d Daily aggregation function (for sub-daily data only). A list indicating the name of the
-#'  aggregation function in first place, and other optional arguments to be passed to the aggregation function. See the examples.
+#'  aggregation function in first place, and other optional arguments to be passed to the aggregation function.
+#'  To be on the safe side, the function in \code{FUN} should be always indicated as a character string. See the examples.
 #' @param aggr.m Same as \code{aggr.d}, but indicating the monthly aggregation function. 
 #' @param aggr.y Same as \code{aggr.d}, but indicating the annual aggregation function. 
 #' @param aggr.s Same as \code{aggr.d}, but indicating the seasonal aggregation function. The season can be indicated
 #' as shown in this example:  aggr.s = list(FUN = list("mean", na.rm = TRUE), season = c(12,1,2)) 
 #' @param aggr.mem Same as \code{aggr.d}, but indicating the function for computing the member aggregation.
-#' @param aggr.lat Same as \code{aggr.d}, indicating the aggregation function to be applied along latitude.
+#' @param aggr.spatial Same as \code{aggr.d}, but indicating the aggregation function in case of rectangular domains to be aggregated
+#' as a unique time series grid (or multimember time series grid.)
+#' @param aggr.lat Same as \code{aggr.d}, indicating the aggregation function to be applied along latitude only.
 #' @param weight.by.lat Logical. Should latitudinal averages be weighted by the cosine of latitude?.
-#' Default to \code{TRUE}. Ignored if no \code{aggr.lat} function is indicated, or a function different from \code{"mean"}
-#' is applied.
+#' Default to \code{TRUE}. Ignored if no \code{aggr.lat} or \code{aggr.spatial} function is indicated,
+#'  or a function different from \code{"mean"} is applied.
 #' @param aggr.lon Same as \code{aggr.lat}, but for longitude.
 #' @param aggr.loc Same as \code{aggr.d}, indicating the aggregation function to be applied along the loc dimension.
 #' @template templateParallelParams
@@ -95,6 +98,7 @@ aggregateGrid <- function(grid,
                           aggr.m = list(FUN = NULL),
                           aggr.y = list(FUN = NULL),
                           aggr.s = list(FUN = NULL, season = NULL),
+                          aggr.spatial = list(FUN = NULL),
                           aggr.lat = list(FUN = NULL),
                           weight.by.lat = TRUE,
                           aggr.lon = list(FUN = NULL),
@@ -124,6 +128,11 @@ aggregateGrid <- function(grid,
                 climatology(clim.fun = aggr.s$FUN)  
         }) %>% bindGrid(dimension = "time")
     }
+    if (!is.null(aggr.spatial$FUN)) {
+        grid <- spatialAggregation(grid, aggr.fun = aggr.spatial,
+                                   weight.by.lat = weight.by.lat,
+                                   parallel, max.ncores, ncores)
+    }
     if (!is.null(aggr.lat$FUN)) {
         grid <- latAggregation(grid, aggr.lat, weight.by.lat, parallel, max.ncores, ncores)
     }
@@ -275,6 +284,7 @@ timeAggregation <- function(grid, aggr.type = c("DD","MM","YY"), aggr.fun, paral
 
 
 latAggregation <- function(grid, aggr.fun, weight.by.lat, parallel, max.ncores, ncores) {
+    stopifnot(is.logical(weight.by.lat))
     dimNames <- getDim(grid)
     if (!"lat" %in% dimNames) {
         message("There is not lat dimension: 'aggr.lat' option was ignored.")
@@ -380,3 +390,48 @@ locAggregation <- function(grid, aggr.fun, weight.by.lat, parallel, max.ncores,
       }
       return(grid)
 }
+
+#' @title Spatial Aggregation
+#' @description Spatial aggregation for rectangular domains. 
+#' @param grid Input grid. 
+#' @param aggr.fun Aggregation function
+#' @param weight.by.lat Logical flag
+#' @return A spatially averaged time series grid (possibly multimember)
+#' @keywords internal
+#' @author J Bedia
+
+spatialAggregation <- function(grid, aggr.fun, weight.by.lat, parallel, max.ncores, ncores) {
+  stopifnot(is.logical(weight.by.lat))
+  dimNames <- getDim(grid)
+  if (!any(c("lon", "lat") %in% dimNames)) {
+    message("\'lon\' and/or \'lat\' dimensions are missing in the input grid: 'aggr.spatial' option was ignored.")
+  } else {
+    if (isTRUE(weight.by.lat)) {
+      lat.weights <- latWeighting(grid)
+      weight.matrix <- matrix(lat.weights, nrow = getShape(grid, "lat"), ncol = getShape(grid, "lon"))
+      if (aggr.fun[["FUN"]] == "mean") {
+        message("Calculating areal weights...")
+        aggr.fun <- list(FUN = "weighted.mean", w = weight.matrix, na.rm = TRUE)
+      } else {
+        message("Spatial weighting skipped: It only applies to \'mean\' aggregation function")
+      }
+    }
+    parallel.pars <- parallelCheck(parallel, max.ncores, ncores)
+    mar <- grep("lat|lon", dimNames, invert = TRUE)
+    aggr.fun[["MARGIN"]] <- mar
+    aggr.fun[["X"]] <- grid$Data
+    out <- if (parallel.pars$hasparallel) {
+      message("[", Sys.time(), "] - Aggregating lat dimension in parallel...")
+      on.exit(parallel::stopCluster(parallel.pars$cl))
+      aggr.fun[["cl"]] <- parallel.pars$cl
+      do.call("parApply", aggr.fun)
+    } else {
+      message("[", Sys.time(), "] - Aggregating spatially...")
+      do.call("apply", aggr.fun)
+    }
+    grid$Data <- out
+    attr(grid$Data, "dimensions") <- dimNames[mar]
+    message("[", Sys.time(), "] - Done.")
+  }
+  return(grid)
+}
diff --git a/R/binaryGrid.R b/R/binaryGrid.R
index 4020111..1faa4b7 100644
--- a/R/binaryGrid.R
+++ b/R/binaryGrid.R
@@ -72,7 +72,7 @@ binaryGrid <- function(x,
         if (isRegular(x)) {
           xx <- suppressWarnings(array3Dto2Dmat(redim(subsetGrid(x,members = j), member = FALSE)$Data))
         } else {
-          xx <- x$Data[j,,]
+          xx <- x$Data[j,,] %>% as.matrix()
         }
         s <- matrix(runif(nrow(xx)*ncol(xx),min = 0,max = 1),nrow = nrow(xx), ncol = ncol(xx))
         xbin <- (xx > s)*1
diff --git a/R/bindGrid.R b/R/bindGrid.R
index c66fb4e..da011ed 100644
--- a/R/bindGrid.R
+++ b/R/bindGrid.R
@@ -1,6 +1,6 @@
 #     bindGrid.R Grid binding along user-defined dimension
 #
-#     Copyright (C) 2020 Santander Meteorology Group (http://www.meteo.unican.es)
+#     Copyright (C) 2021 Santander Meteorology Group (http://www.meteo.unican.es)
 #
 #     This program is free software: you can redistribute it and/or modify
 #     it under the terms of the GNU General Public License as published by
@@ -354,51 +354,51 @@ sortDim.spatial <- function(grid, dimension = c("y", "x")) {
 #' @author M De Felice, J Bedia
 
 bindGrid.time <- function(..., tol) {
-  grid.list <- list(...)
-  if (length(grid.list) == 1) {
-    grid.list <- unlist(grid.list, recursive = FALSE)
-  }
-  if (length(grid.list) < 2) {
-    stop("The input must be a list of at least two grids")
-  }
-  grid.list <- lapply(grid.list, function(i){
-    loc <- !isRegular(i)
-    redim(i, loc = loc, var = TRUE)
-  })
-    for (i in 2:length(grid.list)) {
+   grid.list <- list(...)
+   if (length(grid.list) == 1) {
+      grid.list <- unlist(grid.list, recursive = FALSE)
+   }
+   if (length(grid.list) < 2) {
+      stop("The input must be a list of at least two grids")
+   }
+   grid.list <- lapply(grid.list, function(i){
+      loc <- !isRegular(i)
+      redim(i, loc = loc, var = TRUE)
+   })
+   for (i in 2:length(grid.list)) {
       # Spatial test
       if (!isTRUE(all.equal(grid.list[[1]]$xyCoords, grid.list[[i]]$xyCoords, check.attributes = FALSE, tolerance = tol))) {
-        stop("Input data is not spatially consistent")
+         stop("Input data is not spatially consistent")
       }
       # Member
       if (getShape(grid.list[[1]])[match('member', getDim(grid.list[[1]]))] != getShape(grid.list[[i]])[match('member', getDim(grid.list[[i]]))]) {
-        stop("Member dimension is not spatially consistent")
+         stop("Member dimension is not spatially consistent")
       }
-    }
-    ref <- grid.list[[1]]
-    dimNames <- getDim(ref) 
-    dim.bind <- grep("^time", dimNames)
-    data.list <- lapply(grid.list, FUN = "[[", "Data")
-    ref[["Data"]] <- unname(do.call("abind", c(data.list, along = dim.bind)))
-    data.list <- NULL
-    start.list <- lapply(grid.list, FUN = function(x) {
+   }
+   ref <- grid.list[[1]]
+   dimNames <- getDim(ref) 
+   dim.bind <- grep("^time", dimNames)
+   data.list <- lapply(grid.list, FUN = "[[", "Data")
+   ref[["Data"]] <- unname(do.call("abind", c(data.list, along = dim.bind)))
+   data.list <- NULL
+   start.list <- lapply(grid.list, FUN = function(x) {
       getRefDates(x)
-    })
-    end.list <- lapply(grid.list, FUN = function(x) {
+   })
+   end.list <- lapply(grid.list, FUN = function(x) {
       getRefDates(x, "end")
-    })
-    grid.list <- NULL
-    refdates <- list(start = do.call(c, start.list),
-                     end = do.call(c, end.list))
-    attr(ref[["Data"]], "dimensions") <- dimNames
-    n.vars <- getShape(ref, "var")
-    if (n.vars > 1) refdates <- rep(list(refdates), n.vars)
-    ref[["Dates"]] <- refdates
-    ref <- tryCatch({sortDim.time(ref)}, error = function(err) {
+   })
+   grid.list <- NULL
+   refdates <- list(start = do.call("c", start.list),
+                    end = do.call("c", end.list))
+   attr(ref[["Data"]], "dimensions") <- dimNames
+   n.vars <- getShape(ref, "var")
+   if (n.vars > 1) refdates <- rep(list(refdates), n.vars)
+   ref[["Dates"]] <- refdates
+   ref <- tryCatch({sortDim.time(ref)}, error = function(err) {
       warning("time dimension could not be sorted!")
       ref})
-    ref <- redim(ref, drop = TRUE)
-    return(ref)
+   ref <- redim(ref, drop = TRUE)
+   return(ref)
 }
 #end
 
diff --git a/R/fillGrid.R b/R/fillGrid.R
new file mode 100644
index 0000000..48ff528
--- /dev/null
+++ b/R/fillGrid.R
@@ -0,0 +1,127 @@
+#' @title Fill missing dates and extend the latitude-longitude domain of a given grid.
+#' @description fill with NA missing dates in grids and station datasets
+#' @param grid grid or station data. The lonLim and latLim arguments are only valid for gridded data.
+#' @param tz Optional. Time zone. See Details in \code{\link{fillGridDates}}.
+#' @param lonLim Optional. A vector with the minimum and maximum longitude boundaries to be filled with NAs. Default to lonLim = c(-180,180).
+#' @param latLim Optional. A vector with the minimum and maximum latitude boundaries to be filled with NAs. Default to latLim = c(-90,90).
+#' @return A grid filled with NAs in the previously missing date positions and/or in the latitude-longitude domain indicated.
+#' @author J. Baño-Medina
+#' @export
+#' @examples \donttest{
+#' require(climate4R.datasets) 
+#' require(visualizeR)
+#' grid <- get(data("NCEP_Iberia_psl"))
+#' spatialPlot(climatology(grid), backdrop.theme = "coastline")
+#' grid <- fillGrid(grid, tz = NULL, lonLim = c(-180,180), latLim = c(-90,90))
+#' spatialPlot(climatology(grid), backdrop.theme = "coastline")
+#' }
+fillGrid <- function(grid, tz = "", lonLim = c(-180,180), latLim = c(-90,90)) {
+  if (!is.null(tz)) grid <- fillGridDates(tz = tz)
+  if (!is.null(lonLim)) grid <- fillGridSpatial(grid, lonLim = lonLim, latLim = latLim)
+  return(grid)
+}
+
+
+#' @title Fill missing dates
+#' @description fill with NA missing dates in grids and station datasets
+#' @param grid grid or station data
+#' @param lonLim Optional. A vector with the minimum and maximum longitude boundaries to be filled with NAs. Default to lonLim = c(-180,180).
+#' @param latLim Optional. A vector with the minimum and maximum latitude boundaries to be filled with NAs. Default to latLim = c(-90,90).
+#' @return A grid filled with NAs in the previously missing date positions
+#' @author J. Baño-Medina
+#' @export
+
+fillGridSpatial <- function(grid, lonLim = c(-180,180), latLim = c(-90,90)) {
+  resX <- attr(grid$xyCoords,"resX")
+  resY <- attr(grid$xyCoords,"resY")
+  lonLimGrid <- c(grid$xyCoords$x[1],grid$xyCoords$x[length(grid$xyCoords$x)])
+  latLimGrid <- c(grid$xyCoords$y[1],grid$xyCoords$y[length(grid$xyCoords$y)])
+  
+  la <- rev(seq(lonLimGrid[1],lonLim[1],-resX)[-1])
+  a <- array(dim = c(getShape(grid,"time"),length(grid$xyCoords$y),length(la)))
+  lb <- seq(lonLimGrid[2],lonLim[2],resX)[-1]
+  b <- array(dim = c(getShape(grid,"time"),length(grid$xyCoords$y),length(lb)))
+  lonLimFinal <- c(la,grid$xyCoords$x,lb) %>% unique()
+  
+  lc <- rev(seq(latLimGrid[1],latLim[1],-resY)[-1])
+  c <- array(dim = c(getShape(grid,"time"),length(lc),length(lonLimFinal)))
+  ld <- seq(latLimGrid[2],latLim[2],resY)[-1]
+  d <- array(dim = c(getShape(grid,"time"),length(ld),length(lonLimFinal)))
+  latLimFinal <- c(lc,grid$xyCoords$y,ld) %>% unique()
+  
+  grid %<>% redim()
+  grid <- lapply(1:getShape(grid, "member"), FUN = function(z) {
+    grid %<>% subsetGrid(members = z) %>% redim(member = FALSE)
+    dimNames <- attr(grid$Data,"dimensions")
+    ab <- do.call("abind",list(a,grid$Data,b))
+    grid$Data <- do.call("abind",list(c,ab,d,"along" = 2))
+    attr(grid$Data,"dimensions") <- dimNames
+    attr(grid$Data,"dimnames") <- NULL
+    
+    grid$xyCoords$x <- lonLimFinal
+    grid$xyCoords$y <- latLimFinal
+    return(grid)
+  }) %>% bindGrid(dimension = "member")
+}
+
+
+
+#' #' @title Fill missing dates
+#' #' @description fill with NA missing dates in grids and station datasets
+#' #' @param grid grid or station data
+#' #' @param tz Optional. Time zone. See Details.
+#' #' @details The function attempts to recover the time zone of the input grid when this is correctly defined.
+#' #' Otherwise, the function will leave it as unknown. See \code{\link{timezones}} for more details.
+#' #' @return A grid filled with NAs in the previously missing date positions
+#' #' @author M Iturbide
+#' #' @export
+#' 
+#' 
+#' fillGridDates <- function(grid, tz = "") {
+#'   station <- ("loc" %in% getDim(grid)) 
+#'   grid <- redim(grid, runtime = TRUE, var = TRUE)
+#'   start <- getRefDates(grid)
+#'   end <- getRefDates(grid, which = "end")
+#'   day.step <- as.numeric(names(which.max(table(difftime(c(start, NA), c(NA, start), units = "days")))))
+#'   message("Time difference of ", day.step, " days")
+#'   formato <- "%Y-%m-%d %H:%M:%S"
+#'   if (day.step >= 1) formato <- "%Y-%m-%d"
+#'   tz <- attr(start[1], "tzone")
+#'   usetz <- TRUE
+#'   if (is.null(tz)) {
+#'     tz <- ""
+#'     usetz <- FALSE
+#'     warning("Undefined time zone")
+#'   }
+#'   start <- as.POSIXlt(start, format = formato, tz = tz)
+#'   end <- as.POSIXlt(end, format = formato, tz = tz)
+#'   xs <- as.POSIXlt(as.character(seq.POSIXt(start[1], start[length(start)],
+#'                                            by = day.step*24*60*60)),
+#'                    format = formato, tz = tz)
+#'   xe <- as.POSIXlt(as.character(seq.POSIXt(end[1], end[length(end)],
+#'                                            by = day.step*24*60*60)),
+#'                    format = formato, tz = tz)
+#'   end <- NULL
+#'   test <- data.frame("date" = start, "wh" = TRUE)
+#'   start <- NULL
+#'   result <- merge(data.frame("date" = xs), test,
+#'                   by.y = "date", by.x = "date", all.x = TRUE)
+#'   ind <- which(result[, "wh"])
+#'   sh <- getShape(grid)
+#'   sh[names(sh) == "time"] <- nrow(result)
+#'   result <- NULL
+#'   arr <- array(data = NA, dim = sh)
+#'   arr[,,, ind ,,] <- grid[["Data"]] 
+#'   grid[["Data"]] <- arr
+#'   arr <- NULL
+#'   attr(grid[["Data"]], "dimensions") <- names(sh)
+#'   grid[["Dates"]][["start"]] <- strftime(xs, format = formato, tz = tz, usetz = usetz)
+#'   grid[["Dates"]][["end"]] <- strftime(xe, format = formato, tz = tz, usetz = usetz)
+#'   xs <- xe <- NULL
+#'   grid <- redim(grid, drop = TRUE, loc = station)
+#'   return(grid)
+#' }
+#' 
+#' # end
+
+  
\ No newline at end of file
diff --git a/R/fillGridDates.R b/R/fillGridDates.R
index 9561c22..15835ed 100644
--- a/R/fillGridDates.R
+++ b/R/fillGridDates.R
@@ -1,35 +1,56 @@
+#     fillGridDates.R Continuous time series filling the missing data with NA when necessary
+#
+#     Copyright (C) 2021 Santander Meteorology Group (http://www.meteo.unican.es)
+#
+#     This program is free software: you can redistribute it and/or modify
+#     it under the terms of the GNU General Public License as published by
+#     the Free Software Foundation, either version 3 of the License, or
+#     (at your option) any later version.
+# 
+#     This program is distributed in the hope that it will be useful,
+#     but WITHOUT ANY WARRANTY; without even the implied warranty of
+#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#     GNU General Public License for more details.
+# 
+#     You should have received a copy of the GNU General Public License
+#     along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
 #' @title Fill missing dates
 #' @description fill with NA missing dates in grids and station datasets
 #' @param grid grid or station data
-#' @param tz Optional. Time zone (by default CET is used, see e.g. \code{\{as.POSIXlt}).
+#' @param tz Optional. Time zone. See Details.
+#' @details The function attempts to recover the time zone of the input grid when this is correctly defined.
+#' Otherwise, the function will leave it as unknown. See \code{\link{timezones}} for more details.
 #' @return A grid filled with NAs in the previously missing date positions
 #' @author M Iturbide
+#' @family dateutils
 #' @export
 
 
-fillGridDates <- function(grid, tz = ""){
-  station <- ("loc" %in% getDim(grid)) 
+fillGridDates <- function(grid, tz = "") {
+  station <- ("loc" %in% getDim(grid))
+  grid <- setGridDates.asPOSIXlt(grid)
   grid <- redim(grid, runtime = TRUE, var = TRUE)
   start <- getRefDates(grid)
   end <- getRefDates(grid, which = "end")
-  day.step <- as.numeric(names(which.max(table(difftime(c(start, NA), c(NA, start), units = "days")))))
-  message("Time difference of ", day.step, " days")
-  formato <- "%Y-%m-%d %H:%M:%S"
-  if (day.step >= 1) formato <- "%Y-%m-%d"
-  start <- as.POSIXlt(start, format = formato, tz = tz)
-  end <- as.POSIXlt(end, format = formato, tz = tz)
-  xs <- as.POSIXlt(as.character(seq.POSIXt(start[1], start[length(start)],
-                                           by = day.step*24*60*60)),
-                   format = formato, tz = tz)
-  xe <- as.POSIXlt(as.character(seq.POSIXt(end[1], end[length(end)],
-                                           by = day.step*24*60*60)),
-                   format = formato, tz = tz)
+  timeres <- getTimeResolution(grid)
+  if (timeres == "unknown") stop("Unknown grid temporal resolution")
+  by <- switch(timeres,
+               "1h" = "hour",
+               "3h" = 3600*3,
+               "6h" = 3600*6,
+               "12h" = 3600*12,
+               "DD" = "day",
+               "MM" = "month",
+               "YY" = "year")
+  xs <- seq.POSIXt(from = start[1], to = start[length(start)], by = by)
+  xe <- seq.POSIXt(from = end[1], to = end[length(end)], by = by)
   end <- NULL
   test <- data.frame("date" = start, "wh" = TRUE)
   start <- NULL
-  result <- merge(data.frame("date" = xs), test,
-                  by.y = "date", by.x = "date", all.x = TRUE)
-  ind <- which(result[, "wh"])
+  result <- merge(data.frame("date" = xs), test, by.y = "date", by.x = "date", all.x = TRUE)
+  ind <- which(result[ , "wh"])
   sh <- getShape(grid)
   sh[names(sh) == "time"] <- nrow(result)
   result <- NULL
@@ -38,12 +59,144 @@ fillGridDates <- function(grid, tz = ""){
   grid[["Data"]] <- arr
   arr <- NULL
   attr(grid[["Data"]], "dimensions") <- names(sh)
-  grid[["Dates"]][["start"]] <- strftime(xs, format = formato, tz = tz, usetz = TRUE)
-  grid[["Dates"]][["end"]] <- strftime(xe, format = formato, tz = tz, usetz = TRUE)
+  grid[["Dates"]][["start"]] <- xs
+  grid[["Dates"]][["end"]] <- xe
   xs <- xe <- NULL
   grid <- redim(grid, drop = TRUE, loc = station)
   return(grid)
 }
-
 # end
-  
\ No newline at end of file
+
+
+#' @title Set Grid dates as POSIXlt (possibly with user-defined time zone)
+#' @description Internal utility for date format conversion to POSIX.lt
+#' @details The function attempts to identify the format of the grid dates and 
+#' to perform an adequate conversion to POSIXlt-class, including the time zone identification.
+#' 
+#' If no tz is specified, the function tries to "guess" it from the input. 
+#' Otherwise it is set to the current time zone, with a warning.
+#' 
+#' If no hour-minute-sec information is stored, the output will read 00:00:00
+#' 
+#' @param grid Input C4R grid
+#' @param tz Time zone to be coerced to. Default to \code{NULL}, 
+#' that will preserve the actual tz if defined, or leave it as undefined otherwise.
+#' @return The C4R grid with transformed dates
+#' @keywords internal
+#' @author J Bedia
+#' @family dateutils
+#' @examples \dontrun{
+#' library(climate4R.datasets)
+#' data("CMIP5_Iberia_hus850.rcp85")
+#' # Dates are represented as character string
+#' class(getRefDates(CMIP5_Iberia_hus850.rcp85))
+#' # Conversion to POSIXlt
+#' a <- setGridDates.asPOSIXlt(CMIP5_Iberia_hus850.rcp85)
+#' range(getRefDates(a))
+#' class(getRefDates(a))
+#' ## WARNING: Forcing a time zone will change the tz attribute, but not the time
+#' # (i.e.: it will alter the actual times):
+#' b <- setGridDates.asPOSIXlt(CMIP5_Iberia_hus850.rcp85, tz = "Europe/Madrid")
+#' range(getRefDates(b))
+#' # Therefore, tz should be used ONLY when the time zone is missing or wrongly specified in the 
+#' # input grid and we are sure about the actual TZ of our data
+#' # It may also be used to change the representation of the time zone, e.g., changing 
+#' # the denomination GMT to the equivalent UTC:
+#' c <- setGridDates.asPOSIXlt(CMIP5_Iberia_hus850.rcp85, tz = "UTC")
+#' range(getRefDates(c))
+#' }
+
+
+setGridDates.asPOSIXlt <- function(grid, tz = "") {
+  ds <- getRefDates(grid)
+  de <- getRefDates(grid)
+  dateclass <- class(ds)
+  ref <- ds[1]
+  format <- "%Y-%m-%d %H:%M:%S"
+  
+  if(any(grepl("POSIXlt", dateclass))) {
+    ## Retrieve tz from original data
+    if (tz == "") tz <- ref$zone
+    grid$Dates$start <- as.POSIXlt(ds, tz = tz, format = format)
+    grid$Dates$end <- as.POSIXlt(de, tz = tz, format = format)
+    
+  } else if (any(grepl("POSIXct", dateclass))) {
+    ## Retrieve tz from original data
+    if (tz == "") tz <- attr(ref, "tzone")
+    grid$Dates$start <- as.POSIXlt.POSIXct(ds, tz = tz, format = format)
+    grid$Dates$end <- as.POSIXlt.POSIXct(de, tz = tz, format = format)
+    
+  } else if (dateclass == "character" | dateclass == "array") {
+    if (tz == "") {# Try to guess
+      message("[", Sys.time(), "] Trying to determine the time zone...")
+      # If dates are defined as a character string, somehow we need to "guess" the format
+      # This is done by counting the number of empty spaces between character strings assuming 
+      # Year:Month:day Time TimeZone
+      split.string <- strsplit(ref, split = "\\s")[[1]]
+      # First character of the first part of the date string (should be ALWAYS a number)
+      if (isAlphaCharacter(getFirstChar(split.string[1]))) {
+        stop("Unrecognized Date Format")
+      }
+      # If the second part of the string starts with a number, it is a time definition:
+      if (is.na(split.string[2])) {# There is only a date with no more info attached to it
+        ## TZ is unknown
+        tz <- ""
+      } else if (isAlphaCharacter(getFirstChar(split.string[2]))) {# It is a time def
+        tz <- ""    
+      } else {# Should be a time zone def
+        # In  this case, the time zone should not be coercible to character.
+        # Time zones should never begin with a number, or at least this is the assumption here (see 'OlsonNames()')
+        tz <- split.string[2]
+      }
+      # The third part of the string must be a time zone definition
+      if (!is.na(split.string[3])) {
+        tz <- split.string[3]
+      } else {
+        tz <- ""
+      } 
+      if (tz == "") {
+        warning("[", Sys.time(), "] Time zone unknown. It was set to the current auto-detected time zone (",
+                as.POSIXlt(Sys.time())$zone,")")
+      } else {
+        message("[", Sys.time(), "] Time zone identified and set to ", tz, "\nSee \'setGridDates.asPOSIXlt\' to change the time zone")
+      }
+    }
+    grid$Dates$start <- as.POSIXlt.character(ds, tz = tz, format = format)
+    grid$Dates$end <- as.POSIXlt.character(de, tz = tz, format = format)
+  }
+  return(grid)
+}
+
+
+#' @title Get first character string before punctuation or space
+#' @description Utility for parsing date elements as character strings
+#' @param x character string
+#' @return First part of a character string before any punctuation or space character
+#' @keywords internal
+#' @author J Bedia
+#' @family dateutils
+#' @examples \dontrun{ 
+#'  getFirstChar("2008 Feb 12")
+#'  getFirstChar("2008-02-12")
+#'  getFirstChar("12:00:00 GMT")
+#'  getFirstChar("Etc/UTC")
+#'  }
+
+getFirstChar <- function(x) {
+  strsplit(x, split = "[[:punct:]]|[[:space:]]")[[1]][1]
+}
+
+
+#' @title Is alphabetic character
+#' @description Is alphabetic character? Internal utility for parsing dates
+#' @param x Character string
+#' @return Logical
+#' @keywords internal
+#' @author J Bedia
+#' @family dateutils
+
+isAlphaCharacter <- function(x) {
+  is.na(suppressWarnings(as.numeric(x)))
+}
+
+
diff --git a/R/helpers.R b/R/helpers.R
index fed689b..997383b 100644
--- a/R/helpers.R
+++ b/R/helpers.R
@@ -76,8 +76,18 @@ getGrid <- function(gridData) {
             }
       }else{
             if ("lon" %in% names(gridData$xyCoords)) rot <- TRUE
-            grid.x <- c(gridData$xyCoords$x[1], tail(gridData$xyCoords$x, 1))
-            grid.y <- c(gridData$xyCoords$y[1], tail(gridData$xyCoords$y, 1))
+            if (length(gridData$xyCoords$x) > 1){
+                  grid.x <- c(gridData$xyCoords$x[1], tail(gridData$xyCoords$x, 1))
+            }else{
+                  grid.x <- gridData$xyCoords$x
+            }
+            if (length(gridData$xyCoords$y) > 1){
+                  grid.y <- c(gridData$xyCoords$y[1], tail(gridData$xyCoords$y, 1))
+            }else{
+                  grid.y <- gridData$xyCoords$y
+            }
+            ## grid.x <- c(gridData$xyCoords$x[1], tail(gridData$xyCoords$x, 1))
+            ## grid.y <- c(gridData$xyCoords$y[1], tail(gridData$xyCoords$y, 1))
             out <- list(x = grid.x, y = grid.y)
             if (rot) {
                   out$lon <- gridData$xyCoords$lon 
@@ -648,7 +658,7 @@ isRegular <- function(grid) {
 #' \code{"1h", "3h", "6h", "12h", "DD", "MM", "YY"}. If none of these matches, \code{"unknown"} is returned 
 #' @author J Bedia
 #' @keywords internal
-#' @family get.helpers
+#' @family get.helpers dateutils
 #' @export
 #' @examples \donttest{ 
 #' require(climate4R.datasets)
diff --git a/R/interpGrid.R b/R/interpGrid.R
index b86b013..a4dbd64 100644
--- a/R/interpGrid.R
+++ b/R/interpGrid.R
@@ -162,6 +162,14 @@ interpGrid <- function(grid,
               y = outer(coords$y, coords$x*0, FUN = "+"))$x
     y <- list(x = outer(coords$y*0, coords$x, FUN = "+"),
               y = outer(coords$y, coords$x*0, FUN = "+"))$y
+    if ((length(coords$x) == 1) | (length(coords$y) == 1)){
+      if ((length(x) != getShape(grid,"lon")) & (isRegular(getGrid(grid)))){
+        x <- coords$x
+      }    
+      if ((length(y) != getShape(grid,"lat")) & (isRegular(getGrid(grid)))){
+        y <- coords$y
+      }
+    }
   }
   # New coordinates
   if (is.null(new.coordinates)) {
@@ -291,7 +299,7 @@ interpGrid <- function(grid,
               ind.NN.y[k,l] <- aux.ind[1,1] 
             } else {
               ind.NN.x[k,l] <- 1
-              ind.NN.y[k,l] <- aux.ind
+              ind.NN.y[k,l] <- aux.ind[1]
             }
           } else {
             warning("There are not values to interpolate.")
diff --git a/R/lambWT.R b/R/lambWT.R
index f21ef17..0ea9921 100644
--- a/R/lambWT.R
+++ b/R/lambWT.R
@@ -1,6 +1,6 @@
 #     lambWT.R Calculation of the Weather types (WT) circulation indices from grid
 #
-#     Copyright (C) 2019 Santander Meteorology Group (http://www.meteo.unican.es)
+#     Copyright (C) 2021 Santander Meteorology Group (http://www.meteo.unican.es)
 #
 #     This program is free software: you can redistribute it and/or modify
 #     it under the terms of the GNU General Public License as published by
@@ -20,6 +20,7 @@
 #' @param grid A grid (gridded or station dataset), or multimember grid object of MSLP values.
 #' @param center.point A two value vector that must include lon and lat from a location that will work as center point for the Lamb WT.
 #' See details. 
+#' @param typeU Logical. Should "Unclassified" type be included in the output. Default to \code{FALSE}, so the closest is assigned.
 #' @details According to Trigo and daCamara (2000), Int J Climatol, Lamb WT is only applied on North Atlantic domain. 
 #' The input grid units must be Pa, not hPa/mbar. If it is not in Pa, the units must be converted.
 #' A center location point must be specified by the user. Then, the function calculates from left to right and from first to 16st 
diff --git a/R/matchStations.R b/R/matchStations.R
index 396c151..6f2c1b7 100644
--- a/R/matchStations.R
+++ b/R/matchStations.R
@@ -29,6 +29,6 @@ matchStations <- function(obj,ref) {
   ind <- match(ref$Metadata$station_id,obj$Metadata$station_id)
   nStations <- dim(obj$Data)[which(getDim(obj) == "loc")]
   out <- lapply(1:nStations, FUN = function(z) {
-    subsetDimension(obj,dimension="loc",indices=ind[z])
+    subsetGrid(obj,station.id = ind[z])
   }) %>% bindGrid(dimension = "loc") %>% redim(drop = TRUE)
 }
diff --git a/R/upscaleGrid.R b/R/upscaleGrid.R
index 447dac5..072c2e6 100644
--- a/R/upscaleGrid.R
+++ b/R/upscaleGrid.R
@@ -44,11 +44,13 @@
 upscaleGrid <- function(grid, times = 5,
                         aggr.fun = list(FUN = max, na.rm = TRUE)) {
       x <- grid$xyCoords$x
-      fac <- rep(1:floor(length(x)/times), each = times)
-      indfac <- length(x) - length(fac)
-      fac <- c(fac, rep(max(fac) + 1, indfac))
+      fac0 <- rep(1:floor(length(x)/times), each = times)
+      nfac <- length(x) - length(fac0)
+      indfac <- max(fac0) + 1
+      fac <- c(fac0, rep(indfac, nfac))
       coords <- lapply(split(x, fac), function(k) range(k))
       newcoords <- unlist(lapply(split(x, fac), function(k) mean(k)))
+      if (nfac != 0)  newcoords[indfac] <- newcoords[(indfac - 1)] + (newcoords[(indfac - 1)] - newcoords[(indfac - 2)])
       grid.list <- lapply(coords, function(k) subsetGrid(grid, lonLim = k))
       suppressMessages(suppressWarnings(
             grid.list.lon <- lapply(grid.list, function(k) aggregateGrid(k, aggr.lon = aggr.fun))
@@ -56,19 +58,21 @@ upscaleGrid <- function(grid, times = 5,
       grid <- bindGrid(grid.list.lon, dimension = "lon")
       grid$xyCoords$x <- unname(newcoords)
       y <- grid$xyCoords$y
-      fac <- rep(1:floor(length(y)/times), each = times)
-      indfac <- length(y) - length(fac)
-      fac <- c(fac, rep(max(fac) + 1, indfac))
+      fac0 <- rep(1:floor(length(y)/times), each = times)
+      nfac <- length(y) - length(fac0)
+      indfac <- max(fac0) + 1
+      fac <- c(fac0, rep(indfac, nfac))
       coords <- lapply(split(y, fac), function(k) range(k))
       newcoords <- unlist(lapply(split(y, fac), function(k) mean(k)))
+      if (nfac != 0)  newcoords[indfac] <- newcoords[(indfac - 1)] + (newcoords[(indfac - 1)] - newcoords[(indfac - 2)])
       grid.list <- lapply(coords, function(k) subsetGrid(grid, latLim = k))
       suppressMessages(suppressWarnings(
             grid.list.lat <- lapply(grid.list, function(k) aggregateGrid(k, aggr.lat = aggr.fun, weight.by.lat = FALSE))
       ))
       grid <- bindGrid(grid.list.lat, dimension = "lat")
       grid$xyCoords$y <- unname(newcoords)
-      if(!is.null(attr(grid$xyCoords, "resX"))) attr(grid$xyCoords, "resX") <- attr(grid$xyCoords, "resX") * times
-      if(!is.null(attr(grid$xyCoords, "resY"))) attr(grid$xyCoords, "resY") <- attr(grid$xyCoords, "resY") * times
+      if (!is.null(attr(grid$xyCoords, "resX"))) attr(grid$xyCoords, "resX") <- attr(grid$xyCoords, "resX") * times
+      if (!is.null(attr(grid$xyCoords, "resY"))) attr(grid$xyCoords, "resY") <- attr(grid$xyCoords, "resY") * times
       return(grid)
 }
 
diff --git a/man/EOF2clim.Rd b/man/EOF2clim.Rd
index 9962608..ed7d558 100644
--- a/man/EOF2clim.Rd
+++ b/man/EOF2clim.Rd
@@ -4,16 +4,19 @@
 \alias{EOF2clim}
 \title{PCA object EOFs to climatology grid}
 \usage{
-EOF2clim(prinCompObj, ind.var, member, n.eofs)
+EOF2clim(prinCompObj, ind.var, member = 1L, n.eofs = 1L, rescale = TRUE)
 }
 \arguments{
 \item{prinCompObj}{PCA object}
 
 \item{ind.var}{index position of the target variable}
 
-\item{member}{index position of the member whose EOFs are to be displayed}
+\item{member}{Integer. index position of the member whose EOFs are to be displayed.
+Default to 1 (ignored if no members available)}
 
-\item{n.eofs}{Number of EOFs to be displayed (from 1 to n.eofs)}
+\item{n.eofs}{Integer. Number of EOFs to be displayed (from 1 to \code{n.eofs}). Default to one -first- EOF.}
+
+\item{rescale}{Logical flag. Default to \code{TRUE}, which rescales the EOF to the original input variable units.}
 }
 \value{
 A climatological grid. Note that EOFs are treated as members
@@ -21,6 +24,20 @@ A climatological grid. Note that EOFs are treated as members
 \description{
 Converts a \code{prinComp} output to a climatology of EOFs.
 }
+\examples{
+\dontrun{
+require(climate4R.datasets)
+data("NCEP_Iberia_hus850", "NCEP_Iberia_psl", "NCEP_Iberia_ta850")
+multigrid <- makeMultiGrid(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
+# In this example, we retain the PCs explaining the 99\\% of the variance
+pca <- prinComp(multigrid, v.exp = c(.95,0.90,.90), keep.orig = FALSE)
+require(visualizeR)
+# Original EOF
+EOF2clim(pca, ind.var = "psl") \%>\% spatialPlot(backdrop.theme = "coastline")
+# Rescaled EOF
+EOF2clim(pca, ind.var = "psl", rescale = FALSE) \%>\% spatialPlot(backdrop.theme = "coastline")
+}
+}
 \author{
 J Bedia
 }
diff --git a/man/aggregateGrid.Rd b/man/aggregateGrid.Rd
index 2003f72..a1b4fac 100644
--- a/man/aggregateGrid.Rd
+++ b/man/aggregateGrid.Rd
@@ -11,6 +11,7 @@ aggregateGrid(
   aggr.m = list(FUN = NULL),
   aggr.y = list(FUN = NULL),
   aggr.s = list(FUN = NULL, season = NULL),
+  aggr.spatial = list(FUN = NULL),
   aggr.lat = list(FUN = NULL),
   weight.by.lat = TRUE,
   aggr.lon = list(FUN = NULL),
@@ -26,7 +27,8 @@ aggregateGrid(
 \item{aggr.mem}{Same as \code{aggr.d}, but indicating the function for computing the member aggregation.}
 
 \item{aggr.d}{Daily aggregation function (for sub-daily data only). A list indicating the name of the
-aggregation function in first place, and other optional arguments to be passed to the aggregation function. See the examples.}
+aggregation function in first place, and other optional arguments to be passed to the aggregation function.
+To be on the safe side, the function in \code{FUN} should be always indicated as a character string. See the examples.}
 
 \item{aggr.m}{Same as \code{aggr.d}, but indicating the monthly aggregation function.}
 
@@ -35,11 +37,14 @@ aggregation function in first place, and other optional arguments to be passed t
 \item{aggr.s}{Same as \code{aggr.d}, but indicating the seasonal aggregation function. The season can be indicated
 as shown in this example:  aggr.s = list(FUN = list("mean", na.rm = TRUE), season = c(12,1,2))}
 
-\item{aggr.lat}{Same as \code{aggr.d}, indicating the aggregation function to be applied along latitude.}
+\item{aggr.spatial}{Same as \code{aggr.d}, but indicating the aggregation function in case of rectangular domains to be aggregated
+as a unique time series grid (or multimember time series grid.)}
+
+\item{aggr.lat}{Same as \code{aggr.d}, indicating the aggregation function to be applied along latitude only.}
 
 \item{weight.by.lat}{Logical. Should latitudinal averages be weighted by the cosine of latitude?.
-Default to \code{TRUE}. Ignored if no \code{aggr.lat} function is indicated, or a function different from \code{"mean"}
-is applied.}
+Default to \code{TRUE}. Ignored if no \code{aggr.lat} or \code{aggr.spatial} function is indicated,
+ or a function different from \code{"mean"} is applied.}
 
 \item{aggr.lon}{Same as \code{aggr.lat}, but for longitude.}
 
diff --git a/man/fillGrid.Rd b/man/fillGrid.Rd
new file mode 100644
index 0000000..1b6bd90
--- /dev/null
+++ b/man/fillGrid.Rd
@@ -0,0 +1,36 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fillGrid.R
+\name{fillGrid}
+\alias{fillGrid}
+\title{Fill missing dates and extend the latitude-longitude domain of a given grid.}
+\usage{
+fillGrid(grid, tz = "", lonLim = c(-180, 180), latLim = c(-90, 90))
+}
+\arguments{
+\item{grid}{grid or station data. The lonLim and latLim arguments are only valid for gridded data.}
+
+\item{tz}{Optional. Time zone. See Details in \code{\link{fillGridDates}}.}
+
+\item{lonLim}{Optional. A vector with the minimum and maximum longitude boundaries to be filled with NAs. Default to lonLim = c(-180,180).}
+
+\item{latLim}{Optional. A vector with the minimum and maximum latitude boundaries to be filled with NAs. Default to latLim = c(-90,90).}
+}
+\value{
+A grid filled with NAs in the previously missing date positions and/or in the latitude-longitude domain indicated.
+}
+\description{
+fill with NA missing dates in grids and station datasets
+}
+\examples{
+\donttest{
+require(climate4R.datasets) 
+require(visualizeR)
+grid <- get(data("NCEP_Iberia_psl"))
+spatialPlot(climatology(grid), backdrop.theme = "coastline")
+grid <- fillGrid(grid, tz = NULL, lonLim = c(-180,180), latLim = c(-90,90))
+spatialPlot(climatology(grid), backdrop.theme = "coastline")
+}
+}
+\author{
+J. Baño-Medina
+}
diff --git a/man/fillGridDates.Rd b/man/fillGridDates.Rd
index 96e32ed..c607a1e 100644
--- a/man/fillGridDates.Rd
+++ b/man/fillGridDates.Rd
@@ -9,7 +9,7 @@ fillGridDates(grid, tz = "")
 \arguments{
 \item{grid}{grid or station data}
 
-\item{tz}{Optional. Time zone (by default CET is used, see e.g. \code{\{as.POSIXlt}).}
+\item{tz}{Optional. Time zone. See Details.}
 }
 \value{
 A grid filled with NAs in the previously missing date positions
@@ -17,6 +17,17 @@ A grid filled with NAs in the previously missing date positions
 \description{
 fill with NA missing dates in grids and station datasets
 }
+\details{
+The function attempts to recover the time zone of the input grid when this is correctly defined.
+Otherwise, the function will leave it as unknown. See \code{\link{timezones}} for more details.
+}
+\seealso{
+Other dateutils: 
+\code{\link{getFirstChar}()},
+\code{\link{isAlphaCharacter}()},
+\code{\link{setGridDates.asPOSIXlt}()}
+}
 \author{
 M Iturbide
 }
+\concept{dateutils}
diff --git a/man/fillGridSpatial.Rd b/man/fillGridSpatial.Rd
new file mode 100644
index 0000000..9c5b8cc
--- /dev/null
+++ b/man/fillGridSpatial.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fillGrid.R
+\name{fillGridSpatial}
+\alias{fillGridSpatial}
+\title{Fill missing dates}
+\usage{
+fillGridSpatial(grid, lonLim = c(-180, 180), latLim = c(-90, 90))
+}
+\arguments{
+\item{grid}{grid or station data}
+
+\item{lonLim}{Optional. A vector with the minimum and maximum longitude boundaries to be filled with NAs. Default to lonLim = c(-180,180).}
+
+\item{latLim}{Optional. A vector with the minimum and maximum latitude boundaries to be filled with NAs. Default to latLim = c(-90,90).}
+}
+\value{
+A grid filled with NAs in the previously missing date positions
+}
+\description{
+fill with NA missing dates in grids and station datasets
+}
+\author{
+J. Baño-Medina
+}
diff --git a/man/get2DmatCoordinates.Rd b/man/get2DmatCoordinates.Rd
index b206584..1cb84c8 100644
--- a/man/get2DmatCoordinates.Rd
+++ b/man/get2DmatCoordinates.Rd
@@ -26,7 +26,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/getCoordinates.Rd b/man/getCoordinates.Rd
index 6827a01..d74e17c 100644
--- a/man/getCoordinates.Rd
+++ b/man/getCoordinates.Rd
@@ -41,7 +41,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/getDim.Rd b/man/getDim.Rd
index 3309eaf..129457b 100644
--- a/man/getDim.Rd
+++ b/man/getDim.Rd
@@ -37,7 +37,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/getFirstChar.Rd b/man/getFirstChar.Rd
new file mode 100644
index 0000000..7f90421
--- /dev/null
+++ b/man/getFirstChar.Rd
@@ -0,0 +1,36 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fillGridDates.R
+\name{getFirstChar}
+\alias{getFirstChar}
+\title{Get first character string before punctuation or space}
+\usage{
+getFirstChar(x)
+}
+\arguments{
+\item{x}{character string}
+}
+\value{
+First part of a character string before any punctuation or space character
+}
+\description{
+Utility for parsing date elements as character strings
+}
+\examples{
+\dontrun{ 
+ getFirstChar("2008 Feb 12")
+ getFirstChar("2008-02-12")
+ getFirstChar("12:00:00 GMT")
+ getFirstChar("Etc/UTC")
+ }
+}
+\seealso{
+Other dateutils: 
+\code{\link{fillGridDates}()},
+\code{\link{isAlphaCharacter}()},
+\code{\link{setGridDates.asPOSIXlt}()}
+}
+\author{
+J Bedia
+}
+\concept{dateutils}
+\keyword{internal}
diff --git a/man/getGrid.Rd b/man/getGrid.Rd
index e6fe0b3..8e0be49 100644
--- a/man/getGrid.Rd
+++ b/man/getGrid.Rd
@@ -63,7 +63,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/getGridProj.Rd b/man/getGridProj.Rd
index f9ef8e4..6f3b851 100644
--- a/man/getGridProj.Rd
+++ b/man/getGridProj.Rd
@@ -33,7 +33,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/getGridVerticalLevels.Rd b/man/getGridVerticalLevels.Rd
index 3213b61..345e700 100644
--- a/man/getGridVerticalLevels.Rd
+++ b/man/getGridVerticalLevels.Rd
@@ -51,7 +51,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/getRefDates.Rd b/man/getRefDates.Rd
index 3bd5bf6..770a782 100644
--- a/man/getRefDates.Rd
+++ b/man/getRefDates.Rd
@@ -47,7 +47,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/getSeason.Rd b/man/getSeason.Rd
index 7705c73..1918744 100644
--- a/man/getSeason.Rd
+++ b/man/getSeason.Rd
@@ -33,7 +33,6 @@ Other get.helpers:
 \code{\link{getRefDates}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/getShape.Rd b/man/getShape.Rd
index 855f6e2..80e56fe 100644
--- a/man/getShape.Rd
+++ b/man/getShape.Rd
@@ -37,7 +37,6 @@ Other get.helpers:
 \code{\link{getRefDates}()},
 \code{\link{getSeason}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/getStationID.Rd b/man/getStationID.Rd
index bdb5917..c2e4333 100644
--- a/man/getStationID.Rd
+++ b/man/getStationID.Rd
@@ -26,7 +26,6 @@ Other get.helpers:
 \code{\link{getRefDates}()},
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/getTimeResolution.Rd b/man/getTimeResolution.Rd
index dd22106..2234c7e 100644
--- a/man/getTimeResolution.Rd
+++ b/man/getTimeResolution.Rd
@@ -27,27 +27,8 @@ annual.grid <- aggregateGrid(monthly.grid, aggr.y = list(FUN = "mean", na.rm = T
 stopifnot(identical(getTimeResolution(annual.grid), "YY"))
 }
 }
-\seealso{
-Other get.helpers: 
-\code{\link{get2DmatCoordinates}()},
-\code{\link{getCoordinates}()},
-\code{\link{getDim}()},
-\code{\link{getGridProj}()},
-\code{\link{getGridVerticalLevels}()},
-\code{\link{getGrid}()},
-\code{\link{getRefDates}()},
-\code{\link{getSeason}()},
-\code{\link{getShape}()},
-\code{\link{getStationID}()},
-\code{\link{getVarNames}()},
-\code{\link{getWT}()},
-\code{\link{getYearsAsINDEX}()},
-\code{\link{setGridProj}()},
-\code{\link{typeofGrid}()},
-\code{\link{which.leap}()}
-}
 \author{
 J Bedia
 }
-\concept{get.helpers}
+\concept{get.helpers dateutils}
 \keyword{internal}
diff --git a/man/getVarNames.Rd b/man/getVarNames.Rd
index 9194c47..225e8f1 100644
--- a/man/getVarNames.Rd
+++ b/man/getVarNames.Rd
@@ -70,7 +70,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
 \code{\link{setGridProj}()},
diff --git a/man/getWT.Rd b/man/getWT.Rd
index 6ba2f2f..3246b08 100644
--- a/man/getWT.Rd
+++ b/man/getWT.Rd
@@ -27,7 +27,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getYearsAsINDEX}()},
 \code{\link{setGridProj}()},
diff --git a/man/getYearsAsINDEX.Rd b/man/getYearsAsINDEX.Rd
index f41aa66..eb29c47 100644
--- a/man/getYearsAsINDEX.Rd
+++ b/man/getYearsAsINDEX.Rd
@@ -70,7 +70,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{setGridProj}()},
diff --git a/man/isAlphaCharacter.Rd b/man/isAlphaCharacter.Rd
new file mode 100644
index 0000000..a037626
--- /dev/null
+++ b/man/isAlphaCharacter.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fillGridDates.R
+\name{isAlphaCharacter}
+\alias{isAlphaCharacter}
+\title{Is alphabetic character}
+\usage{
+isAlphaCharacter(x)
+}
+\arguments{
+\item{x}{Character string}
+}
+\value{
+Logical
+}
+\description{
+Is alphabetic character? Internal utility for parsing dates
+}
+\seealso{
+Other dateutils: 
+\code{\link{fillGridDates}()},
+\code{\link{getFirstChar}()},
+\code{\link{setGridDates.asPOSIXlt}()}
+}
+\author{
+J Bedia
+}
+\concept{dateutils}
+\keyword{internal}
diff --git a/man/lambWT.Rd b/man/lambWT.Rd
index fcdf945..b6c7b73 100644
--- a/man/lambWT.Rd
+++ b/man/lambWT.Rd
@@ -4,13 +4,15 @@
 \alias{lambWT}
 \title{Calculation of Lamb Weather types (WT).}
 \usage{
-lambWT(grid, center.point = c(-5, 55))
+lambWT(grid, center.point = c(-5, 55), typeU = FALSE)
 }
 \arguments{
 \item{grid}{A grid (gridded or station dataset), or multimember grid object of MSLP values.}
 
 \item{center.point}{A two value vector that must include lon and lat from a location that will work as center point for the Lamb WT.
 See details.}
+
+\item{typeU}{Logical. Should "Unclassified" type be included in the output. Default to \code{FALSE}, so the closest is assigned.}
 }
 \value{
 The Lamb WT circulation index (and members, if applicable) with:
diff --git a/man/setGridDates.asPOSIXlt.Rd b/man/setGridDates.asPOSIXlt.Rd
new file mode 100644
index 0000000..c8ddbdb
--- /dev/null
+++ b/man/setGridDates.asPOSIXlt.Rd
@@ -0,0 +1,62 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fillGridDates.R
+\name{setGridDates.asPOSIXlt}
+\alias{setGridDates.asPOSIXlt}
+\title{Set Grid dates as POSIXlt (possibly with user-defined time zone)}
+\usage{
+setGridDates.asPOSIXlt(grid, tz = "")
+}
+\arguments{
+\item{grid}{Input C4R grid}
+
+\item{tz}{Time zone to be coerced to. Default to \code{NULL}, 
+that will preserve the actual tz if defined, or leave it as undefined otherwise.}
+}
+\value{
+The C4R grid with transformed dates
+}
+\description{
+Internal utility for date format conversion to POSIX.lt
+}
+\details{
+The function attempts to identify the format of the grid dates and 
+to perform an adequate conversion to POSIXlt-class, including the time zone identification.
+
+If no tz is specified, the function tries to "guess" it from the input. 
+Otherwise it is set to the current time zone, with a warning.
+
+If no hour-minute-sec information is stored, the output will read 00:00:00
+}
+\examples{
+\dontrun{
+library(climate4R.datasets)
+data("CMIP5_Iberia_hus850.rcp85")
+# Dates are represented as character string
+class(getRefDates(CMIP5_Iberia_hus850.rcp85))
+# Conversion to POSIXlt
+a <- setGridDates.asPOSIXlt(CMIP5_Iberia_hus850.rcp85)
+range(getRefDates(a))
+class(getRefDates(a))
+## WARNING: Forcing a time zone will change the tz attribute, but not the time
+# (i.e.: it will alter the actual times):
+b <- setGridDates.asPOSIXlt(CMIP5_Iberia_hus850.rcp85, tz = "Europe/Madrid")
+range(getRefDates(b))
+# Therefore, tz should be used ONLY when the time zone is missing or wrongly specified in the 
+# input grid and we are sure about the actual TZ of our data
+# It may also be used to change the representation of the time zone, e.g., changing 
+# the denomination GMT to the equivalent UTC:
+c <- setGridDates.asPOSIXlt(CMIP5_Iberia_hus850.rcp85, tz = "UTC")
+range(getRefDates(c))
+}
+}
+\seealso{
+Other dateutils: 
+\code{\link{fillGridDates}()},
+\code{\link{getFirstChar}()},
+\code{\link{isAlphaCharacter}()}
+}
+\author{
+J Bedia
+}
+\concept{dateutils}
+\keyword{internal}
diff --git a/man/setGridProj.Rd b/man/setGridProj.Rd
index 697258b..0f6d79e 100644
--- a/man/setGridProj.Rd
+++ b/man/setGridProj.Rd
@@ -44,7 +44,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/spatialAggregation.Rd b/man/spatialAggregation.Rd
new file mode 100644
index 0000000..8bf116e
--- /dev/null
+++ b/man/spatialAggregation.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aggregateGrid.R
+\name{spatialAggregation}
+\alias{spatialAggregation}
+\title{Spatial Aggregation}
+\usage{
+spatialAggregation(grid, aggr.fun, weight.by.lat, parallel, max.ncores, ncores)
+}
+\arguments{
+\item{grid}{Input grid.}
+
+\item{aggr.fun}{Aggregation function}
+
+\item{weight.by.lat}{Logical flag}
+}
+\value{
+A spatially averaged time series grid (possibly multimember)
+}
+\description{
+Spatial aggregation for rectangular domains.
+}
+\author{
+J Bedia
+}
+\keyword{internal}
diff --git a/man/typeofGrid.Rd b/man/typeofGrid.Rd
index 0b79a8f..8f615a3 100644
--- a/man/typeofGrid.Rd
+++ b/man/typeofGrid.Rd
@@ -27,7 +27,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},
diff --git a/man/which.leap.Rd b/man/which.leap.Rd
index 56a2c38..b09fa56 100644
--- a/man/which.leap.Rd
+++ b/man/which.leap.Rd
@@ -34,7 +34,6 @@ Other get.helpers:
 \code{\link{getSeason}()},
 \code{\link{getShape}()},
 \code{\link{getStationID}()},
-\code{\link{getTimeResolution}()},
 \code{\link{getVarNames}()},
 \code{\link{getWT}()},
 \code{\link{getYearsAsINDEX}()},