diff --git a/.Rbuildignore b/.Rbuildignore
index 19a1af65ac5ae3e96e68397e9d01cc76b0c9e8ac..2814b11b1187b7963ea8c1cf6a04e019ebe4392b 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -9,4 +9,8 @@ README\.md$
 \..*\.RData$
 vignettes
 .gitlab-ci.yml
+# unit tests should be ignored when building the package for CRAN
 ^tests$
+# CDO is not in windbuilder, so we can test the unit tests by winbuilder
+# but test-CDORemap.R needs to be hidden
+#tests/testthat/test-CDORemap.R
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index cbe09a22ac03d801b1ae1a2c0f07cc6c4b59f8ef..3eb90d50935404fa29137db76f00eae65305193a 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -4,7 +4,7 @@ build:
   stage: build
   script:
     - module load R/3.6.1-foss-2015a-bare
-#    - module load CDO
+    - module load CDO/1.9.8-foss-2015a
     - R CMD build --resave-data .
     - R CMD check --as-cran --no-manual --run-donttest s2dv_*.tar.gz
-#    - R -e 'covr::package_coverage()'
+    - R -e 'covr::package_coverage()'
diff --git a/DESCRIPTION b/DESCRIPTION
index c2dbf9d5120f8e134ee8b7d3dd91872960c618ff..f2577bfeebb720f2d76e89bde01f267e54ade2ae 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: s2dv
 Title: A Set of Common Tools for Seasonal to Decadal Verification
-Version: 1.0.0
+Version: 1.1.0
 Authors@R: c(
     person("BSC-CNS", role = c("aut", "cph")),
     person("An-Chi", "Ho", , "an.ho@bsc.es", role = c("aut", "cre")),
@@ -26,8 +26,6 @@ Depends:
 Imports:
     abind,
     bigmemory,
-    GEOmap,
-    geomapdata,
     graphics,
     grDevices,
     mapproj,
diff --git a/NAMESPACE b/NAMESPACE
index c02f502122458ebeec850a5c93b498ab9de625f7..ca5d500368ca97734ff3b1c4bbc51b3b547e948c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -55,12 +55,14 @@ export(REOF)
 export(RMS)
 export(RMSSS)
 export(RandomWalkTest)
+export(RatioPredictableComponents)
 export(RatioRMS)
 export(RatioSDRMS)
 export(Regression)
 export(Reorder)
 export(SPOD)
 export(Season)
+export(SignalNoiseRatio)
 export(Smoothing)
 export(Spectrum)
 export(Spread)
@@ -71,11 +73,9 @@ export(Trend)
 export(UltimateBrier)
 export(clim.colors)
 export(clim.palette)
-import(GEOmap)
 import(NbClust)
 import(SpecsVerification)
 import(bigmemory)
-import(geomapdata)
 import(graphics)
 import(mapproj)
 import(maps)
@@ -84,6 +84,7 @@ import(multiApply)
 import(ncdf4)
 import(parallel)
 import(plyr)
+import(stats)
 importFrom(ClimProjDiags,CombineIndices)
 importFrom(ClimProjDiags,Subset)
 importFrom(ClimProjDiags,WeightedMean)
@@ -130,5 +131,6 @@ importFrom(stats,sd)
 importFrom(stats,setNames)
 importFrom(stats,spectrum)
 importFrom(stats,ts)
+importFrom(stats,var)
 importFrom(stats,varimax)
 importFrom(stats,window)
diff --git a/NEWS.md b/NEWS.md
index e537e95964bc8453030f091399a1a66bc81cd3f0..71993f0e554bb260e19b2678f9c455d5e0923f2f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,13 @@
+# s2dv 1.1.0 (Release date: 2021-12-14)
+- New functions: RatioPredictableComponents, SignalNoiseRatio  
+- CDORemap(): Able to interpolate irregular grid to regular grid; include new cdo methods 'con2', 'laf' and 'nn'
+- PlotEquiMap(): Discard the dependency on 'GEOmap' and 'geomapdata' and only use package 'map' now;
+new parameters 'country.borders', 'shapefile', 'shapefile_color', and 'shapefile_lwd' for plotting the national borders and shapefile  
+- PlotLayout(): new parameter 'layout_by_rows' for changing the layout order  
+- MeanDims(): new parameter 'drop' to choose whether to drop the averaged dimension or not;
+Bugfix for making the result as array even if the result is only a number  
+- Season(): Add dimension name even if the result is only a number
+
 # s2dv 1.0.0 (Release date: 2021-06-16)
 - New functions:
 ACC, Ano_CrossValid, BrierScore, CDORemap, Cluster, Consistent_Trend, EOF, EuroAtlanticTC, Filter, Histo2Hindcast, 
diff --git a/R/Ano.R b/R/Ano.R
index e0a69db232230da28d8dde61fc1d92d07f123de5..c4c70a330bf00074632ec18f1b6ea243e621ced8 100644
--- a/R/Ano.R
+++ b/R/Ano.R
@@ -21,10 +21,10 @@
 #'clim <- Clim(sampleData$mod, sampleData$obs)
 #'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
 #'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
-#'\donttest{
+#'\dontrun{
 #'PlotAno(ano_exp, ano_obs, startDates, 
 #'        toptitle = 'Anomaly', ytitle = c('K', 'K', 'K'), 
-#'        legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano.png')
+#'        legends = 'ERSST', biglab = FALSE)
 #'}
 #'@import multiApply
 #'@export
diff --git a/R/Ano_CrossValid.R b/R/Ano_CrossValid.R
index 22e710a1a3ded8b428eb7493a1dd0c0aec489d18..99205020b7412afdf90513ad18dedf0b9302320d 100644
--- a/R/Ano_CrossValid.R
+++ b/R/Ano_CrossValid.R
@@ -47,7 +47,7 @@
 #'\dontrun{
 #'PlotAno(anomalies$exp, anomalies$obs, startDates, 
 #'        toptitle = paste('anomalies'), ytitle = c('K', 'K', 'K'), 
-#'        legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano_crossvalid.eps')
+#'        legends = 'ERSST', biglab = FALSE)
 #'}
 #'@import multiApply
 #'@importFrom ClimProjDiags Subset
diff --git a/R/BrierScore.R b/R/BrierScore.R
index 1363f6155ff53a7c73cb8afd3d630f05fceb343a..c3673ad4273089059f13debc91d244554c2f90cd 100644
--- a/R/BrierScore.R
+++ b/R/BrierScore.R
@@ -174,7 +174,7 @@ BrierScore <- function(exp, obs, thresholds = seq(0.1, 0.9, 0.1), time_dim = 'sd
     name_exp <- name_exp[-which(name_exp == dat_dim)]
     name_obs <- name_obs[-which(name_obs == dat_dim)]
   }
-  if (any(name_exp != name_obs)) {
+  if (any(!name_exp %in% name_obs) | any(!name_obs %in% name_exp)) {
     stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
                 "of all the dimensions expect 'dat_dim' and 'memb_dim'."))
   }
diff --git a/R/CDORemap.R b/R/CDORemap.R
index fc25b527de4036257c462b28a0d5eb27eac89154..00824de207b5649e541451b009da43eafa4ab321 100644
--- a/R/CDORemap.R
+++ b/R/CDORemap.R
@@ -20,8 +20,7 @@
 #'  provided, it must have at least a longitude and a latitude dimensions, 
 #'  identified by the array dimension names. The names for these dimensions 
 #'  must be one of the recognized by s2dverification (can be checked with 
-#'  \code{s2dverification:::.KnownLonNames()} and 
-#'  \code{s2dverification:::.KnownLatNames()}).
+#'  \code{s2dv:::.KnownLonNames()} and \code{s2dv:::.KnownLatNames()}).
 #'@param lons Numeric vector or array of longitudes of the centers of the grid 
 #'  cells. Its size must match the size of the longitude/latitude dimensions 
 #'  of the input array.
@@ -33,9 +32,9 @@
 #'  NetCDF file which to read the target grid from (a single grid must be 
 #'  defined in such file).
 #'@param method Character string specifying an interpolation method 
-#'  (recognized by CDO; e.g.: 'con', 'bil', 'bic', 'dis'). The following 
-#'  long names are also supported: 'conservative', 'bilinear', 'bicubic' and 
-#'  'distance-weighted'.
+#'  (recognized by CDO; e.g.: 'con', 'bil', 'bic', 'dis', 'con2', 'laf', 'nn').
+#'  The following long names are also supported: 'conservative', 'bilinear', 
+#'  'bicubic' and 'distance-weighted'.
 #'@param avoid_writes The step of permutation is needed when the input array 
 #'  has more than 3 dimensions and none of the longitude or latitude dimensions
 #'   in the right-most position (CDO would not accept it without permuting 
@@ -251,22 +250,27 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
       array_dims <- c(length(lats), length(lons))
       new_lon_dim_name <- 'lon'
       new_lat_dim_name <- 'lat'
-    } else {
-      array_dims <- dim(lons)
-      new_lon_dim_name <- 'i'
-      new_lat_dim_name <- 'j'
-    }
-    if (!is.null(names(dim(lons)))) {
-      if (any(known_lon_names %in% names(dim(lons)))) {
-        new_lon_dim_name <- known_lon_names[which(known_lon_names %in% names(dim(lons)))[1]]
+
+      if (!is.null(names(dim(lons)))) {
+        if (any(known_lon_names %in% names(dim(lons)))) {
+          new_lon_dim_name <- known_lon_names[which(known_lon_names %in% names(dim(lons)))[1]]
+        }
       }
-    }
-    if (!is.null(names(dim(lats)))) {
-      if (any(known_lat_names %in% names(dim(lats)))) {
-        new_lat_dim_name <- known_lat_names[which(known_lat_names %in% names(dim(lats)))[1]]
+      if (!is.null(names(dim(lats)))) {
+        if (any(known_lat_names %in% names(dim(lats)))) {
+          new_lat_dim_name <- known_lat_names[which(known_lat_names %in% names(dim(lats)))[1]]
+        }
+      }
+      names(array_dims) <- c(new_lat_dim_name, new_lon_dim_name)
+
+    } else {  # irregular
+      array_dims <- dim(lons)
+      if (is.null(names(array_dims))) {
+        new_lon_dim_name <- 'i'
+        new_lat_dim_name <- 'j'
       }
     }
-    names(array_dims) <- c(new_lat_dim_name, new_lon_dim_name)
+
     data_array <- array(as.numeric(NA), array_dims)
   }
   if (!(is.logical(data_array) || is.numeric(data_array)) || !is.array(data_array)) {
@@ -385,8 +389,14 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
     method <- 'con'
   } else if (method %in% c('dis', 'distance-weighted')) {
     method <- 'dis'
+  } else if (method %in% 'nn') {
+    method <- 'nn'
+  } else if (method %in% 'laf') {
+    method <- 'laf'
+  } else if (method %in% 'con2') {
+    method <- 'con2'
   } else {
-    stop("Unsupported CDO remap method. 'bilinear', 'bicubic', 'conservative' or 'distance-weighted' supported only.")
+    stop("Unsupported CDO remap method. Only 'bilinear', 'bicubic', 'conservative', 'distance-weighted', 'nn', 'laf', and 'con2' are supported.")
   }
   # Check avoid_writes
   if (!is.logical(avoid_writes)) {
@@ -890,7 +900,16 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
       } else {
         new_dims <- dim(data_array)
         new_dims[c(lon_dim, lat_dim)] <- c(found_lon_dim_size, found_lat_dim_size)
-
+        if (is_irregular) {
+          lon_pos <- which(names(new_dims) == lon_dim)
+          lat_pos <- which(names(new_dims) == lat_dim)
+          if (lon_pos > lat_pos) {
+            new_pos <- 1:length(new_dims)
+            new_pos[lon_pos] <- lat_pos
+            new_pos[lat_pos] <- lon_pos
+            new_dims <- new_dims[new_pos]
+          }
+        }
         result_array <- ncvar_get(ncdf_remapped, 'var', collapse_degen = FALSE)
         dim(result_array) <- new_dims
       }
diff --git a/R/Clim.R b/R/Clim.R
index 21f97b67ebcd586efb284f6654523b51f53eadd7..d01cf5d91ec219dd7b9124abccabb392d1f35215 100644
--- a/R/Clim.R
+++ b/R/Clim.R
@@ -56,11 +56,11 @@
 #'example(Load)
 #'clim <- Clim(sampleData$mod, sampleData$obs)
 #'clim2 <- Clim(sampleData$mod, sampleData$obs, method = 'kharin', memb = FALSE)
-#'\donttest{
+#'\dontrun{
 #'PlotClim(clim$clim_exp, clim$clim_obs, 
 #'         toptitle = paste('sea surface temperature climatologies'), 
 #'         ytitle = 'K', monini = 11, listexp = c('CMIP5 IC3'), 
-#'         listobs = c('ERSST'), biglab = FALSE, fileout = 'tos_clim.eps')
+#'         listobs = c('ERSST'), biglab = FALSE)
 #'}
 #'@importFrom abind adrop
 #'@importFrom ClimProjDiags Subset
diff --git a/R/Cluster.R b/R/Cluster.R
index 33527aea9b0e8f43bcce78be874cd07d5a02629a..8f401182b93a1d1b62a1041de9779910cbcc8d2d 100644
--- a/R/Cluster.R
+++ b/R/Cluster.R
@@ -181,7 +181,7 @@ Cluster <- function(data, weights = NULL, time_dim = 'sdate', space_dim = NULL,
       if (!(length(space_dim) == length(dim(weights)) & all(space_dim %in% names(dim(weights))))) {
         stop("Parameter 'weights' must have dimension names the same as 'space_dim'.")
       }
-      if (space_dim != names(dim(weights))) {
+      if (any(space_dim != names(dim(weights)))) {
         space_dim <- names(dim(weights))
       }
     }
diff --git a/R/Histo2Hindcast.R b/R/Histo2Hindcast.R
index 860f56b96bae59a9b06e0eff937a24d5df79668d..5657457d7a47933153496ada5062cb8a94bc95fa 100644
--- a/R/Histo2Hindcast.R
+++ b/R/Histo2Hindcast.R
@@ -66,7 +66,7 @@ Histo2Hindcast <- function(data, sdatesin, sdatesout, nleadtimesout,
   if (is.null(sdatesin)) {
     stop("Parameter 'sdatesin' cannot be NULL.")
   }
-  if (!is.character(sdatesin) | length(sdatesin) > 1) {
+  if (!is.character(sdatesin) || length(sdatesin) > 1) {
     stop(paste0("Parameter 'sdatesin' must be a character string in the format",
                 " 'YYYYMMDD' or 'YYYYMM'."))
   } else if (!nchar(sdatesin) %in% c(6, 8) | is.na(as.numeric(sdatesin))) {
@@ -88,12 +88,12 @@ Histo2Hindcast <- function(data, sdatesin, sdatesout, nleadtimesout,
   if (is.null(nleadtimesout)) {
     stop("Parameter 'nleadtimesout' cannot be NULL.")
   }
-  if (!is.numeric(nleadtimesout) | nleadtimesout %% 1 != 0 | 
-      nleadtimesout < 0 | length(nleadtimesout) > 1) {
+  if (!is.numeric(nleadtimesout) | any(nleadtimesout %% 1 != 0) | 
+      any(nleadtimesout < 0) | length(nleadtimesout) > 1) {
     stop("Parameter 'nleadtimesout' must be a positive integer.")
   }
   # sdate_dim
-  if (!is.character(sdate_dim) | length(sdate_dim) > 1) {
+  if (!is.character(sdate_dim) || length(sdate_dim) > 1) {
     stop("Parameter 'sdate_dim' must be a character string.")
   }
   if (!sdate_dim %in% names(dim(data))) {
@@ -103,7 +103,7 @@ Histo2Hindcast <- function(data, sdatesin, sdatesout, nleadtimesout,
     stop("The dimension length of sdate_dim of 'data' must be 1.")
   }
   # ftime_dim
-  if (!is.character(ftime_dim) | length(ftime_dim) > 1) {
+  if (!is.character(ftime_dim) || length(ftime_dim) > 1) {
     stop("Parameter 'ftime_dim' must be a character string.")
   }
   if (!ftime_dim %in% names(dim(data))) {
@@ -111,7 +111,7 @@ Histo2Hindcast <- function(data, sdatesin, sdatesout, nleadtimesout,
   }
   # ncores
   if (!is.null(ncores)) {
-    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+    if (!is.numeric(ncores) || ncores %% 1 != 0 || ncores <= 0 |
       length(ncores) > 1) {
       stop("Parameter 'ncores' must be a positive integer.")
     }
diff --git a/R/MeanDims.R b/R/MeanDims.R
index 7d3cb445487f60aefb8020f2cb598c28d2777fd2..9e5dd49fe47ed21008ed0c95da4a2f3163d98acd 100644
--- a/R/MeanDims.R
+++ b/R/MeanDims.R
@@ -8,17 +8,21 @@
 #'  dimensions to average.
 #'@param na.rm A logical value indicating whether to ignore NA values (TRUE) or 
 #'  not (FALSE).
-#'@return An array with the same dimension as parameter 'data' except the 'dims' 
-#'  dimensions. 
-#'  removed.
+#'@param drop A logical value indicating whether to keep the averaged 
+#'  dimension (FALSE) or drop it (TRUE). The default value is TRUE.
+#'@return An array with the same dimension as parameter 'data' except the 
+#'  'dims' dimensions. If 'drop' is TRUE, 'dims' will be removed; if 'drop' is
+#'  FALSE, 'dims' will be preserved and the length will be 1.
 #'
 #'@examples
-#'a <- array(rnorm(24), dim = c(2, 3, 4))
-#'MeanDims(a, 2)
-#'MeanDims(a, c(2, 3))
+#'a <- array(rnorm(24), dim = c(dat = 2, member = 3, time = 4))
+#'ens_mean <- MeanDims(a, 'member')
+#'dim(ens_mean)
+#'ens_time_mean <- MeanDims(a, c(2, 3), drop = FALSE)
+#'dim(ens_time_mean)
 #'@import multiApply
 #'@export
-MeanDims <- function(data, dims, na.rm = FALSE) {
+MeanDims <- function(data, dims, na.rm = FALSE, drop = TRUE) {
 
   # Check inputs 
   ## data
@@ -54,18 +58,42 @@ MeanDims <- function(data, dims, na.rm = FALSE) {
   if (!is.logical(na.rm) | length(na.rm) > 1) {
     stop("Parameter 'na.rm' must be one logical value.")
   }
+  ## drop
+  if (!is.logical(drop) | length(drop) > 1) {
+    stop("Parameter 'drop' must be one logical value.")
+  }
 
   ###############################
   # Calculate MeanDims
-  if (length(dims) == length(dim(data))) {
-    data <- mean(data, na.rm = na.rm)
+  dim_data <- dim(data)
+
+  if (length(dims) == length(dim_data)) {
+    if (drop) {
+      data <- as.array(mean(data, na.rm = na.rm))
+    } else {
+      data <- array(mean(data, na.rm = na.rm), 
+                    dim = rep(1, length(dim_data)))
+      names(dim(data)) <- names(dim_data)
+    }
   } else {
     if (is.character(dims)) {
-      dims <- which(names(dim(data)) %in% dims)
+      dims <- which(names(dim_data) %in% dims)
     }
-    pos <- (1:length(dim(data)))[-dims]
+    pos <- (1:length(dim_data))[-dims]
     data <- apply(data, pos, mean, na.rm = na.rm)
+
+    # If data is vector
+    if (is.null(dim(data))) {
+      data <- array(data, dim = dim_data[-dims])
+    }
+    if (!drop) {
+      restore_dim <- as.array(rep(1, length(dims)))
+      names(restore_dim) <- names(dim_data)[dims]
+      data <- array(data, dim = c(dim(data), restore_dim))
+      data <- Reorder(data, names(dim_data))
+    }
   }
+
   return(data)
 }
 
diff --git a/R/PlotEquiMap.R b/R/PlotEquiMap.R
index f6a85380bedaa0649ae56026c2721a6e308f4f2f..3de30e984afe36fdd29b1ece13069d5f67abd1d1 100644
--- a/R/PlotEquiMap.R
+++ b/R/PlotEquiMap.R
@@ -62,14 +62,23 @@
 #'  Takes the value gray(0.5) by default or, if 'square = FALSE', takes the 
 #'  value FALSE. If set to FALSE, continents are not filled in.
 #'@param filled.oceans A logical value or the color name to fill in drawn 
-#'  projected oceans. The default value is FALSE. If it is TRUE, the default 
+#'  projected oceans. The default value is FALSE. If it is TRUE, the default
 #'  colour is "light blue".
+#'@param country.borders A logical value indicating if the country borders 
+#'  should be plotted (TRUE) or not (FALSE). It only works when 
+#'  'filled.continents' is FALSE. The default value is FALSE.
 #'@param coast_color Colour of the coast line of the drawn projected continents.
 #'   Takes the value gray(0.5) by default.
 #'@param coast_width Line width of the coast line of the drawn projected 
 #'  continents. Takes the value 1 by default.
 #'@param lake_color Colour of the lake or other water body inside continents.
-#'  The default value is NULL. 
+#'  The default value is NULL.
+#'@param shapefile A character string of the path to a .rds file or a list 
+#'  object containinig shape file data. If it is a .rds file, it should contain
+#'  a list. The list should contains 'x' and 'y' at least, which indicate the 
+#'  location of the shape. The default value is NULL.
+#'@param shapefile_color Line color of the shapefile.
+#'@param shapefile_lwd Line width of the shapefile. The default value is 1. 
 #'@param contours Array of same dimensions as 'var' to be added to the plot 
 #'  and displayed with contours. Parameter 'brks2' is required to define the 
 #'  magnitude breaks for each contour curve. Disregarded if 'square = FALSE'.
@@ -218,7 +227,7 @@
 #'PlotEquiMap(sampleData$mod[1, 1, 1, 1, , ], sampleData$lon, sampleData$lat, 
 #'            toptitle = 'Predicted sea surface temperature for Nov 1960 from 1st Nov',
 #'            sizetit = 0.5)
-#'@import graphics GEOmap geomapdata maps
+#'@import graphics maps
 #'@importFrom grDevices dev.cur dev.new dev.off gray
 #'@importFrom stats cor
 #'@export
@@ -228,8 +237,9 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
                         triangle_ends = NULL, col_inf = NULL, col_sup = NULL, 
                         colNA = NULL, color_fun = clim.palette(),
                         square = TRUE, filled.continents = NULL,
-                        filled.oceans = FALSE,
+                        filled.oceans = FALSE, country.borders = FALSE,
                         coast_color = NULL, coast_width = 1, lake_color = NULL,
+                        shapefile = NULL, shapefile_color = NULL, shapefile_lwd = 1,
                         contours = NULL, brks2 = NULL, contour_lwd = 0.5,
                         contour_color = 'black', contour_lty = 1,
                         contour_draw_label = TRUE, contour_label_scale = 1,
@@ -434,6 +444,11 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
     ocean_color <- "light blue"
   }
 
+  # Check country.borders
+  if (!is.logical(country.borders)) {
+    stop("Parameter 'country.borders' must be logical.")
+  }
+
   # Check coast_color
   if (is.null(coast_color)) {
     if (filled.continents) {
@@ -453,15 +468,57 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
 
   # Check lake_color
   if (!is.null(lake_color)) {
-#    if (filled.continents) {
-#      lake_color <- 'white'
-#    }
-#  } else {
     if (!.IsColor(lake_color)) {
       stop("Parameter 'lake_color' must be a valid colour identifier.")
     }
   }
 
+  # Check shapefile
+  if (!is.null(shapefile)) {
+    if (is.list(shapefile)) {
+      shape <- shapefile
+      if (any(!c('x', 'y') %in% names(shape))) {
+        stop("The list names of the object in 'shapefile' .rds file should ", 
+             "have at least 'x' and 'y'.")
+      }
+      if (length(shape$x) != length(shape$y)) {
+        stop("The length of x and y in 'shapefile' list should be equal.")
+      }
+    } else if (!is.character(shapefile)) {
+      stop("Parameter 'shapefile' must be a .rds file or a list.")
+    } else {  # .rds file
+      if (!file.exists(shapefile)) {
+        stop("Parameter 'shapefile' is not a valid file.")
+      }
+      if (!grepl("\\.rds$", shapefile)) {
+        stop("Parameter 'shapefile' must be a .rds file or a list.")
+      }
+      shape <- readRDS(file = shapefile)
+      if (!is.list(shape)) {
+        stop("Parameter 'shapefile' should be a .rds file of a list object.")
+      }
+      if (any(!c('x', 'y') %in% names(shape))) {
+        stop("The list names of the object in 'shapefile' .rds file should ",
+             "have at least 'x' and 'y'.")
+      }
+      if (length(shape$x) != length(shape$y)) {
+        stop("The length of x and y in 'shapefile' list should be equal.")
+      }
+    }
+  }
+
+  # Check shapefile_col
+  if (is.null(shapefile_color)) {
+    if (filled.continents) {
+      shapefile_color <- continent_color
+    } else {
+      shapefile_color <- 'black'
+    }
+  }
+  if (!.IsColor(shapefile_color)) {
+    stop("Parameter 'shapefile_color' must be a valid colour identifier.")
+  }
+
   # Check contours
   if (!is.null(contours)) {
     if (dim(contours)[1] != dims[1] || dim(contours)[2] != dims[2]) {
@@ -656,10 +713,6 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
     }
   }
 
-  #library(GEOmap)
-  #library(geomapdata)
-  #library(maps)
-  utils::data(coastmap, package = 'GEOmap', envir = environment())
   #
   #  Input arguments 
   # ~~~~~~~~~~~~~~~~~
@@ -844,47 +897,21 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
   }
   old_lwd <- par('lwd')
   par(lwd = coast_width)
-  # If [0, 360], use GEOmap; if [-180, 180], use maps
-  if (min(lon) >= 0) {
-    ylat <- latmin:latmax
-    xlon <- lonmin:lonmax
-    proj <- GEOmap::setPROJ(1, LON0 = mean(xlon),
-                            LAT0 = mean(ylat), LATS = ylat, LONS = xlon)
-    lakes <- which(coastmap$STROKES$col == "blue")
-    par(new = TRUE)
-    if (filled.continents) {
-      coastmap$STROKES$col[which(coastmap$STROKES$col != "blue")] <- continent_color
-      if (is.null(lake_color)) {
-        coastmap$STROKES$col[lakes] <- continent_color
-      } else {
-        coastmap$STROKES$col[lakes] <- lake_color #"white"
-      }
-      GEOmap::plotGEOmap(coastmap, PROJ = proj, border = coast_color, 
-                         add = TRUE, lwd = coast_width)
-    } else {
-      coastmap$STROKES$col[which(coastmap$STROKES$col != "blue")] <- coast_color
-      if (is.null(lake_color)) {
-        coastmap$STROKES$col[lakes] <- coast_color
-      } else {
-        coastmap$STROKES$col[lakes] <- lake_color #"white"
-      }
-      GEOmap::plotGEOmap(coastmap, PROJ = proj, #MAPcol = coast_color,
-                         add = TRUE, lwd = coast_width, MAPstyle = 2)
-    }
-
-  } else {
-    # [-180, 180]
-    coast <- map(continents, interior = FALSE, wrap = TRUE,
+  # If [0, 360], use GEOmap; if [-180, 180], use maps::map
+  # UPDATE: Use maps::map for both cases. The difference between GEOmap and
+  #         maps is trivial. The only thing we can see for now is that 
+  #         GEOmap has better lakes.
+    coast <- maps::map(continents, interior = country.borders, wrap = TRUE,
                  xlim = xlim_conti, ylim = c(-89.99, 89.99),
                  fill = filled.continents, add = TRUE, plot = FALSE)
-    if (filled.continents) {
-      polygon(coast, col = continent_color, border = coast_color, lwd = coast_width)
-    } else {
-      lines(coast, col = coast_color, lwd = coast_width)
-    }
-    if (!is.null(lake_color)) {
-      map('lakes', add = TRUE, fill = filled.continents, col = lake_color) 
-    }
+
+  if (filled.continents) {
+    polygon(coast, col = continent_color, border = coast_color, lwd = coast_width)
+  } else {
+    lines(coast, col = coast_color, lwd = coast_width)
+  }
+  if (!is.null(lake_color)) {
+    maps::map('lakes', add = TRUE, fill = filled.continents, col = lake_color) 
   }
   par(lwd = old_lwd)
 
@@ -893,7 +920,7 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
       old_lwd <- par('lwd')
       par(lwd = coast_width)
 
-      outline <- map(continents, fill = T, plot = FALSE)  # must be fill = T
+      outline <- maps::map(continents, fill = T, plot = FALSE)  # must be fill = T
       xbox <- xlim_conti + c(-2, 2)
       ybox <- c(-92, 92) 
       outline$x <- c(outline$x, NA, c(xbox, rev(xbox), xbox[1]))
@@ -903,6 +930,14 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
       par(lwd = old_lwd)
   }
 
+  # Plot shapefile
+  if (!is.null(shapefile)) {
+    maps::map(shape, interior = country.borders, #wrap = TRUE,
+              xlim = xlim_conti, ylim = c(-89.99, 89.99),
+              fill = filled.continents, add = TRUE, plot = TRUE, 
+              lwd = shapefile_lwd, col = shapefile_color)
+  }
+
   box()
   # Draw rectangle on the map
   if (!is.null(boxlim)) {
diff --git a/R/PlotLayout.R b/R/PlotLayout.R
index 0d333cb6a5a67dd3479752c3b4beacaa1ecc948d..50a87f26a5c354ce0e1ed09782c615207ab9a366 100644
--- a/R/PlotLayout.R
+++ b/R/PlotLayout.R
@@ -122,6 +122,8 @@
 #'@param extra_margin Extra margins to be added around the layout, in the 
 #'  format c(y1, x1, y2, x2). The units are margin lines. Takes rep(0, 4) 
 #'  by default.
+#'@param layout_by_rows Logical indicating wether the panels should be filled 
+#'  by columns (FALSE) or by raws (TRUE, default).
 #'@param fileout File where to save the plot. If not specified (default) a 
 #'  graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, 
 #'  bmp and tiff.
@@ -196,7 +198,7 @@
 #'@importFrom grDevices dev.cur dev.new dev.off
 #'@export
 PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL, 
-                       nrow = NULL, ncol = NULL, toptitle = NULL, 
+                       nrow = NULL, ncol = NULL, toptitle = NULL,
                        row_titles = NULL, col_titles = NULL, bar_scale = 1, 
                        title_scale = 1, title_margin_scale = 1,
                        title_left_shift_scale = 1,
@@ -210,11 +212,11 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL,
                        units = NULL, units_scale = 1, bar_label_scale = 1, 
                        bar_tick_scale = 1, bar_extra_margin = rep(0, 4), 
                        bar_left_shift_scale = 1, bar_label_digits = 4, 
-                       extra_margin = rep(0, 4),
+                       extra_margin = rep(0, 4), layout_by_rows = TRUE,
                        fileout = NULL, width = NULL, height = NULL, 
                        size_units = 'in', res = 100, close_device = TRUE) {
   # If there is any filenames to store the graphics, process them
-  # to select the right device 
+  # to select the right device
   if (!is.null(fileout)) {
     deviceInfo <- .SelectDevice(fileout = fileout, width = width, height = height, units = size_units, res = res)
     saveToFile <- deviceInfo$fun
@@ -279,6 +281,10 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL,
     }
     ncol <- round(ncol)
   }
+  # Check layout_by_rows
+  if (!is.logical(layout_by_rows)) {
+     stop("Parameter 'layout_by_rows' must be logical.")
+  }
 
   # Check toptitle
   if (is.null(toptitle) || is.na(toptitle)) {
@@ -522,7 +528,7 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL,
   subtitle_cex <- 1.5 * subtitle_scale
   subtitle_margin <- 0.5 * sqrt(nrow * ncol) * subtitle_cex * subtitle_margin_scale
   mat_layout <- 1:(nrow * ncol) + ifelse(drawleg != FALSE, 1, 0)
-  mat_layout <- matrix(mat_layout, nrow, ncol, byrow = TRUE)
+  mat_layout <- matrix(mat_layout, nrow, ncol, byrow = layout_by_rows)
   fsu <- figure_size_units <- 10  # unitless
   widths <- rep(fsu, ncol)
   heights <- rep(fsu, nrow)
diff --git a/R/ProbBins.R b/R/ProbBins.R
index 327ceb34fafb48b6425373536e2b1bf7d70041dc..ef293d94710fa7d70bbb3a8984a25898ddc68f72 100644
--- a/R/ProbBins.R
+++ b/R/ProbBins.R
@@ -100,7 +100,7 @@ ProbBins <- function(data, thr, fcyr = 'all', time_dim = 'sdate', memb_dim = 'me
     stop("Parameter 'memb_dim' is not found in 'data' dimension.")
   }
   ## fcyr
-  if (fcyr != 'all') {
+  if (!identical(fcyr, 'all')) {
     if (!is.numeric(fcyr) | !is.vector(fcyr)) {
       stop("Parameter 'fcyr' must be a numeric vector or 'all'.")
     } else if (any(fcyr %% 1 != 0) | min(fcyr) < 1 | max(fcyr) > dim(data)[time_dim]) {
diff --git a/R/REOF.R b/R/REOF.R
index c9c82cf94e1091f88ea7a8c71dd2957e5079de80..2def222614b84ff7f15a772eb5092ff22691b6ce 100644
--- a/R/REOF.R
+++ b/R/REOF.R
@@ -146,9 +146,9 @@ REOF <- function(ano, lat, lon, ntrunc = 15, time_dim = 'sdate',
   # ntrunc is bounded
   if (ntrunc != min(dim(ano)[time_dim], prod(dim(ano)[space_dim]), ntrunc)) {
     ntrunc <- min(dim(ano)[time_dim], prod(dim(ano)[space_dim]), ntrunc)
-    .warning(paste0("Parameter 'ntrunc' is changed to ", ntrunc, ", the minimum among ",
-                    "the length of time_dim, the production of the length of space_dim, ",
-                    "and ntrunc."))
+    warning(paste0("Parameter 'ntrunc' is changed to ", ntrunc, ", the minimum among ",
+                   "the length of time_dim, the production of the length of space_dim, ",
+                   "and ntrunc."))
   }
 
   # Area weighting is needed to compute the fraction of variance explained by 
diff --git a/R/RatioPredictableComponents.R b/R/RatioPredictableComponents.R
new file mode 100644
index 0000000000000000000000000000000000000000..163d18fcdbd10524d24b6111a7b759527e5f611f
--- /dev/null
+++ b/R/RatioPredictableComponents.R
@@ -0,0 +1,105 @@
+#'Calculate ratio of predictable components (RPC)
+#'
+#'This function computes the ratio of predictable components (RPC; Eade et al., 2014).
+#'
+#'@param exp A numerical array with, at least, 'time_dim' and 'member_dim' 
+#'  dimensions.
+#'@param obs A numerical array with the same dimensions than 'exp' except the
+#'  'member_dim' dimension.
+#'@param time_dim A character string indicating the name of the time dimension.
+#'  The default value is 'year'.
+#'@param member_dim A character string indicating the name of the member 
+#'  dimension. The default value is 'member'.
+#'@param na.rm A logical value indicating whether to remove NA values during
+#'  the computation. The default value is FALSE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return An array of the ratio of the predictable components. it has the same
+#'  dimensions as 'exp' except 'time_dim' and 'member_dim' dimensions.
+#'
+#'@examples
+#'exp <- array(data = runif(600), dim = c(year = 15, member = 10, lat = 2, lon = 2))
+#'obs <- array(data = runif(60), dim = c(year = 15, lat = 2, lon = 2))
+#'RatioPredictableComponents(exp, obs)
+#'
+#'@import multiApply stats
+#'@export
+RatioPredictableComponents <- function(exp, obs, time_dim = 'year', member_dim = 'member', na.rm = FALSE, ncores = NULL) {
+
+  ## Checkings
+  if (is.null(exp)) {
+    stop("Parameter 'exp' cannot be NULL.")
+  }
+  if (!is.numeric(exp)) {
+    stop("Parameter 'exp' must be a numeric array.")
+  }
+  if (is.null(obs)) {
+    stop("Parameter 'obs' cannot be NULL.")
+  }
+  if (!is.numeric(obs)) {
+    stop("Parameter 'obs' must be a numeric array.")
+  }
+  if (!(is.character(time_dim) & length(time_dim) == 1)) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!(is.character(member_dim) & length(member_dim) == 1)) {
+    stop("Parameter 'member_dim' must be a character string.")
+  }
+  if (!time_dim %in% names(dim(exp))) {
+    stop("'exp' must have 'time_dim' dimension.")
+  }
+  if (!member_dim %in% names(dim(exp))) {
+    stop("'exp' must have 'member_dim' dimension.")
+  }
+  if (!time_dim %in% names(dim(obs))) {
+    stop("'obs' must have 'time_dim' dimension.")
+  }
+  if (!identical(dim(exp)[time_dim], dim(obs)[time_dim])) {
+    stop("'exp' and 'obs' must have the same length for 'time_dim'.")
+  }
+  if (!is.logical(na.rm)) {
+    stop("Parameter 'na.rm' must be TRUE or FALSE.")
+  }
+  if (!is.null(ncores)){
+    if (!is.numeric(ncores) | ncores %% 1 != 0 |  ncores <= 0 | length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
+  RPC <- multiApply::Apply(data = list(exp, obs), 
+                           target_dims = list(exp = c(time_dim, member_dim),
+                                              obs = time_dim), 
+                           output_dims = NULL,
+                           fun = .RatioPredictableComponents,
+                           na.rm = na.rm,
+                           ncores = ncores)$output1
+  return(RPC)
+}
+
+.RatioPredictableComponents <- function(exp, obs, na.rm = na.rm) {
+  # exp: [time, member]
+  # obs: [time]
+
+  ## Ensemble mean and spread
+  ens_mean <- apply(exp, 1, mean, na.rm = na.rm)
+  ens_spread <- apply(exp, 2, "-", ens_mean)
+
+  ## Ensemble mean variance -> signal
+  var_signal <- var(ens_mean, na.rm = na.rm)
+  ## Variance of ensemble members about ensemble mean (= spread) -> noise
+  # var_noise <- var(as.vector(ens_spread), na.rm = na.rm)
+
+  ## Total variance
+  # var_total <- var_signal + var_noise
+  var_total <- mean(apply(exp, 2, var, na.rm = na.rm), na.rm = na.rm)
+
+  ## Correlation between observations and the ensemble mean
+  r <- cor(ens_mean, obs, method = 'pearson')
+
+  ## Ratio of predictable components
+  RPC <- r / sqrt(var_signal / var_total)
+  
+  return(RPC)
+}
+
diff --git a/R/RatioSDRMS.R b/R/RatioSDRMS.R
index 4d833d7ff12e05085c9bea899e8950ada803b4fd..544ca6f8db3b6f85c73788fdecd5caa9c8557586 100644
--- a/R/RatioSDRMS.R
+++ b/R/RatioSDRMS.R
@@ -44,11 +44,10 @@
 #'rsdrms_plot <- array(dim = c(dim(rsdrms$ratio)[1:2], 4, dim(rsdrms$ratio)[3]))
 #'rsdrms_plot[, , 2, ] <- rsdrms$ratio
 #'rsdrms_plot[, , 4, ] <- rsdrms$p.val
-#'\donttest{
+#'\dontrun{
 #'PlotVsLTime(rsdrms_plot, toptitle = "Ratio ensemble spread / RMSE", ytitle = "", 
 #'            monini = 11, limits = c(-1, 1.3), listexp = c('CMIP5 IC3'), 
-#'            listobs = c('ERSST'), biglab = FALSE, siglev = TRUE, 
-#'            fileout = 'tos_rsdrms.eps')
+#'            listobs = c('ERSST'), biglab = FALSE, siglev = TRUE) 
 #'}
 #'
 #'@import multiApply
@@ -183,7 +182,7 @@ RatioSDRMS <- function(exp, obs, dat_dim = 'dataset', memb_dim = 'member',
 
       if (pval) {
         F <- (enosd[jexp] * std[jexp]^2 / (enosd[jexp] - 1)) / (enorms * rms^2 / (enorms - 1))
-        if (!is.na(F) & !is.na(enosd) & !is.na(enorms) & enosd > 2 && enorms > 2) {
+        if (!is.na(F) & !is.na(enosd[jexp]) & !is.na(enorms) & any(enosd > 2) & enorms > 2) {
           p.val[jexp, jobs] <- 1 - pf(F, enosd[jexp] - 1, enorms - 1)
         } else {
           ratiosdrms[jexp, jobs] <- NA
diff --git a/R/Season.R b/R/Season.R
index d56aa1637fd77204db76061492a7d14fea9e8a87..086dc529471fb69a385f9c9b21eefc233c29b3c1 100644
--- a/R/Season.R
+++ b/R/Season.R
@@ -136,6 +136,9 @@ Season <- function(data, time_dim = 'ftime', monini, moninf, monsup,
         res <- apply(data, c(1:length(dim(data)))[-time_dim_ind], .Season,
                      monini = monini, moninf = moninf, monsup = monsup,
                      method = method, na.rm = na.rm)
+        if (is.null(dim(res))) {
+          res <- array(res, dim = dim(data)[-time_dim_ind])
+        }
         if (length(dim(res)) < length(dim(data))) {
           res <- InsertDim(res, posdim = 1, lendim = 1, name = time_dim)
         } else {
diff --git a/R/SignalNoiseRatio.R b/R/SignalNoiseRatio.R
new file mode 100644
index 0000000000000000000000000000000000000000..18dccf407a52c7eb3c7d6e193edf0fe2d1e470b5
--- /dev/null
+++ b/R/SignalNoiseRatio.R
@@ -0,0 +1,83 @@
+#'Calculate Signal-to-noise ratio
+#'
+#'This function computes the signal-to-noise ratio, where the signal is the 
+#'ensemble mean variance and the noise is the variance of the ensemble members 
+#'about the ensemble mean (Eade et al., 2014; Scaife and Smith, 2018).
+#'
+#'@param data A numerical array with, at least, 'time_dim' and 'member_dim' 
+#'  dimensions.
+#'@param time_dim A character string indicating the name of the time dimension
+#'  in 'data'. The default value is 'year'.
+#'@param member_dim A character string indicating the name of the member 
+#'  dimension in 'data'. The default value is 'member'.
+#'@param na.rm A logical value indicating whether to remove NA values during the
+#'  computation. The default value is FALSE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return An array with of the signal-to-noise ratio. It has the same dimensions
+#'  as 'data' except 'time_dim' and 'member_dim' dimensions.
+#'
+#'@examples
+#' exp <- array(data = runif(600), dim = c(year = 15, member = 10, lat = 2, lon = 2))
+#' SignalNoiseRatio(exp)
+#'
+#'@import multiApply
+#'@importFrom stats var
+#'@export
+SignalNoiseRatio <- function(data, time_dim = 'year', member_dim = 'member', na.rm = FALSE, ncores = NULL) {
+  
+  ## Input Check
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
+  }
+  if (!is.numeric(data)) {
+    stop("Parameter 'data' must be a numeric array.")
+  }
+  if (!(is.character(time_dim) & length(time_dim) == 1)) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!(is.character(member_dim) & length(member_dim) == 1)) {
+    stop("Parameter 'member_dim' must be a character string.")
+  }
+  if (!time_dim %in% names(dim(data))) {
+    stop("'data' must have 'time_dim' dimension.")
+  }
+  if (!member_dim %in% names(dim(data))) {
+    stop("'data' must have 'member_dim' dimension.")
+  }
+  if (!is.logical(na.rm)) {
+    stop("Parameter 'na.rm' must be TRUE or FALSE.")
+  }
+  if (!is.null(ncores)){
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | length(ncores) > 1){
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
+  SNR <- multiApply::Apply(data = data, 
+                           target_dims = c(time_dim, member_dim), 
+                           output_dims = NULL,
+                           fun = .SignalNoiseRatio,
+                           na.rm = na.rm,
+                           ncores = ncores)$output1
+  return(SNR)
+}
+
+.SignalNoiseRatio <- function(data, na.rm = na.rm) {
+  # data: [time, member]
+
+  ## Ensemble mean and spread
+  ens_mean <- apply(data, 1, mean, na.rm = na.rm)
+  ens_spread <- apply(data, 2, "-", ens_mean)
+
+  ## Ensemble mean variance -> signal
+  var_signal <- stats::var(ens_mean, na.rm = na.rm)
+  ## Variance of ensemble members about ensemble mean (= spread) -> noise
+  var_noise <- stats::var(as.vector(ens_spread), na.rm = na.rm)
+
+  ## Signal to noise ratio
+  SNR <- var_signal / var_noise
+
+  return(SNR)
+}
diff --git a/R/Smoothing.R b/R/Smoothing.R
index d5fd2a5eac95d11d7217b05d2c5b795f8792d71d..1b31e6594858262fbea5353392e04a7435e2aec8 100644
--- a/R/Smoothing.R
+++ b/R/Smoothing.R
@@ -25,10 +25,9 @@
 #'smooth_ano_obs <- Smoothing(ano_obs, time_dim = 'ftime', runmeanlen = 12)
 #'smooth_ano_exp <- Reorder(smooth_ano_exp, c(2, 3, 4, 1))
 #'smooth_ano_obs <- Reorder(smooth_ano_obs, c(2, 3, 4, 1))
-#'  \donttest{
+#'  \dontrun{
 #'PlotAno(smooth_ano_exp, smooth_ano_obs, startDates, 
-#'        toptitle = "Smoothed Mediterranean mean SST", ytitle = "K",
-#'        fileout = "tos_smoothed_ano.png")
+#'        toptitle = "Smoothed Mediterranean mean SST", ytitle = "K")
 #'  }
 #'@import plyr multiApply
 #'@export
diff --git a/R/Spread.R b/R/Spread.R
index 4b3bc6b1167534dc9004a91d99da686851ed4611..d1d8f6d159bceab7781f9d131ad47589381ed708 100644
--- a/R/Spread.R
+++ b/R/Spread.R
@@ -54,27 +54,27 @@
 #'                                                   name = 'member')
 #'spread <- Spread(smooth_ano_exp_m_sub, compute_dim = c('member', 'sdate'))
 #'
-#'\donttest{
+#'\dontrun{
 #'PlotVsLTime(Reorder(spread$iqr, c('dataset', 'stats', 'ftime')), 
 #'            toptitle = "Inter-Quartile Range between ensemble members",
 #'            ytitle = "K", monini = 11, limits = NULL, 
 #'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-#'            hlines = c(0), fileout = 'tos_iqr.png')
+#'            hlines = c(0))
 #'PlotVsLTime(Reorder(spread$maxmin, c('dataset', 'stats', 'ftime')), 
 #'            toptitle = "Maximum minus minimum of the members", 
 #'            ytitle = "K", monini = 11, limits = NULL, 
 #'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-#'            hlines = c(0), fileout = 'tos_maxmin.png')
+#'            hlines = c(0))
 #'PlotVsLTime(Reorder(spread$sd, c('dataset', 'stats', 'ftime')), 
 #'            toptitle = "Standard deviation of the members", 
 #'            ytitle = "K", monini = 11, limits = NULL, 
 #'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-#'            hlines = c(0), fileout = 'tos_sd.png')
+#'            hlines = c(0))
 #'PlotVsLTime(Reorder(spread$mad, c('dataset', 'stats', 'ftime')), 
 #'            toptitle = "Median Absolute Deviation of the members",
 #'            ytitle = "K", monini = 11, limits = NULL, 
 #'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-#'            hlines = c(0), fileout = 'tos_mad.png')
+#'            hlines = c(0))
 #'}
 #'
 #'@import multiApply
@@ -114,12 +114,12 @@ Spread <- function(data, compute_dim = 'member', na.rm = TRUE,
     stop("Parameter 'conf' must be one logical value.")
   }
   ## conf.lev
-  if (!is.numeric(conf.lev) | conf.lev < 0 | conf.lev > 1 | length(conf.lev) > 1) {
+  if (!is.numeric(conf.lev) | any(conf.lev < 0) | any(conf.lev > 1) | length(conf.lev) > 1) {
     stop("Parameter 'conf.lev' must be a numeric number between 0 and 1.")
   }
   ## ncores
   if (!is.null(ncores)) {
-    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores <= 0) |
       length(ncores) > 1) {
       stop("Parameter 'ncores' must be a positive integer.")
     }
diff --git a/R/Utils.R b/R/Utils.R
index e0e1f91fddd2bcd61460149d49060ec65ab4d018..76ef63a63d1510cadf6974f44185a74823e69f7e 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -1671,9 +1671,9 @@
   if (type == 'dcpp') {
     
     fyear_dim <- 'fyear'
-    data <- s2dv::Season(data = data, time_dim = fmonth_dim,
-                         monini = monini, moninf = 1, monsup = 12,
-                         method = mean, na.rm = na.rm)
+    data <- Season(data = data, time_dim = fmonth_dim,
+                   monini = monini, moninf = 1, monsup = 12,
+                   method = mean, na.rm = na.rm)
     names(dim(data))[which(names(dim(data))==fmonth_dim)] <- fyear_dim
     
     if (identical(indices_for_clim, FALSE)) { ## data is already anomalies
diff --git a/inst/doc/FAQ.md b/inst/doc/FAQ.md
index a8e49db9024f6e063a2fbe59db9274a4efd20b18..fe37eccc77209aeb192b53d6240ff46c4d6c9275 100644
--- a/inst/doc/FAQ.md
+++ b/inst/doc/FAQ.md
@@ -5,6 +5,7 @@ This document intends to be the first reference for any doubts that you may have
 ## Index
 1. **How to**
    1. [Global Map with non-standard longitudinal boundaries](#1-global-map-with-non-standard-longitudinal-boundaries)
+   2. [Plot a region boundary from shape file with PlotEquiMap](#2-plot-a-region-boundary-from-shape-file-with-plotequimap)
 
 2. **Something goes wrong...**  
    1. [CDORemap() returns errors or warnings with specific module versions](#1-cdoremap-returns-errors-or-warnings-with-specific-module-versions)
@@ -45,6 +46,49 @@ Note: You can adjust many parameters to visualize the plot, here we are just sho
 If you want to add other information to the plot (e.g.: hatching, points, countours, ...), you can add it just before ColorBar() function.
 
 
+### 2. Plot a region boundary from shape file with PlotEquiMap()
+
+PlotEquiMap() provides the option to plot the country boundaries by argument
+`country.borders = TRUE`. If you want to plot specific area, you can save 
+the shapefile data as an RDS file or read the shapefile as a list, and use 
+argument `shapefile` to plot it. 
+
+If you have more than one shapefile to plot, you have to plot it on top of the 
+figure after PlotEquiMap(). 
+The following script shows an example of plotting the Douro region on top of
+default country boundaries. The data in this example is synthetic.
+To use "rgdal" package, `module load GDAL` first (see wiki [R tips 3. How to load dependencies of R package rgdal](https://earth.bsc.es/wiki/doku.php?id=tools:Rtools&s%5B%5D=Rtools#r_tips) for details.)
+
+```r
+library(rgdal)
+library(sp)
+library(maps)
+library(s2dv)
+
+# Read the Douro shapefile and transform:
+shp1 <- readOGR(dsn = '/esarchive/scratch/cchou/MEDGOLD/others/Shapefiles/', layer = 'RDD_ETRS89')
+s1 <- spTransform(shp1, CRS("+proj=longlat"))
+douro <- SpatialPolygons2map(s1)
+str(douro)
+#List of 4
+# $ x    : num [1:37826] -7.77 -7.77 -7.77 -7.77 -7.77 ...
+# $ y    : num [1:37826] 41.1 41.1 41.1 41.1 41.1 ...
+# $ names: chr [1:10] "0:1" "0:2" "0:3" "0:4" ...
+# $ range: num [1:4] -7.92 -6.75 40.92 41.53
+# - attr(*, "class")= chr "map"
+
+# Synthetic data
+dat <- 1 : (21 * 10)
+dim(dat) <- c(lon = 21, lat = 10)
+
+# Plot Douro region using argument 'shapefile'
+PlotEquiMap(dat, lon = c(1:10, 350:360), lat = 36 : 45, drawleg = FALSE, 
+            filled.continents = FALSE, country.border = TRUE, 
+            shapefile = douro, shapefile_col = 'blue')
+
+# Enhance the border line and change the color
+map(douro, interior = FALSE, add = TRUE, lwd = 2, col = "red")
+```
 
 ## 2. Something goes wrong...  
 
@@ -65,6 +109,3 @@ R/3.6.1-foss-2015a-bare
 HDF5/1.10.5-foss-2015a  
 
 
-
-
-
diff --git a/man/Ano.Rd b/man/Ano.Rd
index d15ffd14bdbfbb52a3975d8146267edb20eca0b9..e30467e0f3070aacae112aa9086322f3320f2903 100644
--- a/man/Ano.Rd
+++ b/man/Ano.Rd
@@ -32,9 +32,9 @@ example(Load)
 clim <- Clim(sampleData$mod, sampleData$obs)
 ano_exp <- Ano(sampleData$mod, clim$clim_exp)
 ano_obs <- Ano(sampleData$obs, clim$clim_obs)
-\donttest{
+\dontrun{
 PlotAno(ano_exp, ano_obs, startDates, 
        toptitle = 'Anomaly', ytitle = c('K', 'K', 'K'), 
-       legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano.png')
+       legends = 'ERSST', biglab = FALSE)
 }
 }
diff --git a/man/Ano_CrossValid.Rd b/man/Ano_CrossValid.Rd
index 1e91528335c2e68b812594311f1446a8d7eabeb3..d2234a1623901efa767aba83bbd4c5eb532940a7 100644
--- a/man/Ano_CrossValid.Rd
+++ b/man/Ano_CrossValid.Rd
@@ -69,6 +69,6 @@ anomalies <- Ano_CrossValid(sampleData$mod, sampleData$obs)
 \dontrun{
 PlotAno(anomalies$exp, anomalies$obs, startDates, 
        toptitle = paste('anomalies'), ytitle = c('K', 'K', 'K'), 
-       legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano_crossvalid.eps')
+       legends = 'ERSST', biglab = FALSE)
 }
 }
diff --git a/man/CDORemap.Rd b/man/CDORemap.Rd
index c7a00a8f7789ec5c8162a2a3b03d0a8a0b5eafd3..b49f94493442665c162a9c10ae595738dbc2baa1 100644
--- a/man/CDORemap.Rd
+++ b/man/CDORemap.Rd
@@ -21,8 +21,7 @@ CDORemap(
 provided, it must have at least a longitude and a latitude dimensions, 
 identified by the array dimension names. The names for these dimensions 
 must be one of the recognized by s2dverification (can be checked with 
-\code{s2dverification:::.KnownLonNames()} and 
-\code{s2dverification:::.KnownLatNames()}).}
+\code{s2dv:::.KnownLonNames()} and \code{s2dv:::.KnownLatNames()}).}
 
 \item{lons}{Numeric vector or array of longitudes of the centers of the grid 
 cells. Its size must match the size of the longitude/latitude dimensions 
@@ -38,9 +37,9 @@ NetCDF file which to read the target grid from (a single grid must be
 defined in such file).}
 
 \item{method}{Character string specifying an interpolation method 
-(recognized by CDO; e.g.: 'con', 'bil', 'bic', 'dis'). The following 
-long names are also supported: 'conservative', 'bilinear', 'bicubic' and 
-'distance-weighted'.}
+(recognized by CDO; e.g.: 'con', 'bil', 'bic', 'dis', 'con2', 'laf', 'nn').
+The following long names are also supported: 'conservative', 'bilinear', 
+'bicubic' and 'distance-weighted'.}
 
 \item{avoid_writes}{The step of permutation is needed when the input array 
 has more than 3 dimensions and none of the longitude or latitude dimensions
diff --git a/man/Clim.Rd b/man/Clim.Rd
index 78559bdbefc9c5253d510f2c7eee1fbc743324aa..50ec0d9c90d7e2dc1b0486de820e657ee84aa325 100644
--- a/man/Clim.Rd
+++ b/man/Clim.Rd
@@ -84,10 +84,10 @@ the 'exp' and 'obs' are excluded when computing the climatologies.
 example(Load)
 clim <- Clim(sampleData$mod, sampleData$obs)
 clim2 <- Clim(sampleData$mod, sampleData$obs, method = 'kharin', memb = FALSE)
-\donttest{
+\dontrun{
 PlotClim(clim$clim_exp, clim$clim_obs, 
         toptitle = paste('sea surface temperature climatologies'), 
         ytitle = 'K', monini = 11, listexp = c('CMIP5 IC3'), 
-        listobs = c('ERSST'), biglab = FALSE, fileout = 'tos_clim.eps')
+        listobs = c('ERSST'), biglab = FALSE)
 }
 }
diff --git a/man/MeanDims.Rd b/man/MeanDims.Rd
index f70b78b391f1205a30614735a8af875a64c9b608..721866b9fbc58225741f89a5866445d600ad1dfb 100644
--- a/man/MeanDims.Rd
+++ b/man/MeanDims.Rd
@@ -4,7 +4,7 @@
 \alias{MeanDims}
 \title{Average an array along multiple dimensions}
 \usage{
-MeanDims(data, dims, na.rm = FALSE)
+MeanDims(data, dims, na.rm = FALSE, drop = TRUE)
 }
 \arguments{
 \item{data}{An array to be averaged.}
@@ -14,18 +14,23 @@ dimensions to average.}
 
 \item{na.rm}{A logical value indicating whether to ignore NA values (TRUE) or 
 not (FALSE).}
+
+\item{drop}{A logical value indicating whether to keep the averaged 
+dimension (FALSE) or drop it (TRUE). The default value is TRUE.}
 }
 \value{
-An array with the same dimension as parameter 'data' except the 'dims' 
- dimensions. 
- removed.
+An array with the same dimension as parameter 'data' except the 
+ 'dims' dimensions. If 'drop' is TRUE, 'dims' will be removed; if 'drop' is
+ FALSE, 'dims' will be preserved and the length will be 1.
 }
 \description{
 This function returns the mean of an array along a set of dimensions and 
 preserves the dimension names if it has.
 }
 \examples{
-a <- array(rnorm(24), dim = c(2, 3, 4))
-MeanDims(a, 2)
-MeanDims(a, c(2, 3))
+a <- array(rnorm(24), dim = c(dat = 2, member = 3, time = 4))
+ens_mean <- MeanDims(a, 'member')
+dim(ens_mean)
+ens_time_mean <- MeanDims(a, c(2, 3), drop = FALSE)
+dim(ens_time_mean)
 }
diff --git a/man/PlotEquiMap.Rd b/man/PlotEquiMap.Rd
index b574904cedbc13e914b8059e21ef3e8e5119b7a1..31fe4d82eeac0bcffe71d6485965632da750be4f 100644
--- a/man/PlotEquiMap.Rd
+++ b/man/PlotEquiMap.Rd
@@ -24,9 +24,13 @@ PlotEquiMap(
   square = TRUE,
   filled.continents = NULL,
   filled.oceans = FALSE,
+  country.borders = FALSE,
   coast_color = NULL,
   coast_width = 1,
   lake_color = NULL,
+  shapefile = NULL,
+  shapefile_color = NULL,
+  shapefile_lwd = 1,
   contours = NULL,
   brks2 = NULL,
   contour_lwd = 0.5,
@@ -143,9 +147,13 @@ Takes the value gray(0.5) by default or, if 'square = FALSE', takes the
 value FALSE. If set to FALSE, continents are not filled in.}
 
 \item{filled.oceans}{A logical value or the color name to fill in drawn 
-projected oceans. The default value is FALSE. If it is TRUE, the default 
+projected oceans. The default value is FALSE. If it is TRUE, the default
 colour is "light blue".}
 
+\item{country.borders}{A logical value indicating if the country borders 
+should be plotted (TRUE) or not (FALSE). It only works when 
+'filled.continents' is FALSE. The default value is FALSE.}
+
 \item{coast_color}{Colour of the coast line of the drawn projected continents.
 Takes the value gray(0.5) by default.}
 
@@ -155,6 +163,15 @@ continents. Takes the value 1 by default.}
 \item{lake_color}{Colour of the lake or other water body inside continents.
 The default value is NULL.}
 
+\item{shapefile}{A character string of the path to a .rds file or a list 
+object containinig shape file data. If it is a .rds file, it should contain
+a list. The list should contains 'x' and 'y' at least, which indicate the 
+location of the shape. The default value is NULL.}
+
+\item{shapefile_color}{Line color of the shapefile.}
+
+\item{shapefile_lwd}{Line width of the shapefile. The default value is 1.}
+
 \item{contours}{Array of same dimensions as 'var' to be added to the plot 
 and displayed with contours. Parameter 'brks2' is required to define the 
 magnitude breaks for each contour curve. Disregarded if 'square = FALSE'.}
diff --git a/man/PlotLayout.Rd b/man/PlotLayout.Rd
index 453cf2e924074f449efd4af697d21a60e3d5afb6..10bb39e94e733ed78e12d10de324f69ddc06530a 100644
--- a/man/PlotLayout.Rd
+++ b/man/PlotLayout.Rd
@@ -43,6 +43,7 @@ PlotLayout(
   bar_left_shift_scale = 1,
   bar_label_digits = 4,
   extra_margin = rep(0, 4),
+  layout_by_rows = TRUE,
   fileout = NULL,
   width = NULL,
   height = NULL,
@@ -180,6 +181,9 @@ disregarded if no 'row_titles' are provided.}
 format c(y1, x1, y2, x2). The units are margin lines. Takes rep(0, 4) 
 by default.}
 
+\item{layout_by_rows}{Logical indicating wether the panels should be filled 
+by columns (FALSE) or by raws (TRUE, default).}
+
 \item{fileout}{File where to save the plot. If not specified (default) a 
 graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, 
 bmp and tiff.}
diff --git a/man/RatioPredictableComponents.Rd b/man/RatioPredictableComponents.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..3e7fbad6c000b42fc6fb1bd689db8007d250b853
--- /dev/null
+++ b/man/RatioPredictableComponents.Rd
@@ -0,0 +1,47 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RatioPredictableComponents.R
+\name{RatioPredictableComponents}
+\alias{RatioPredictableComponents}
+\title{Calculate ratio of predictable components (RPC)}
+\usage{
+RatioPredictableComponents(
+  exp,
+  obs,
+  time_dim = "year",
+  member_dim = "member",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A numerical array with, at least, 'time_dim' and 'member_dim' 
+dimensions.}
+
+\item{obs}{A numerical array with the same dimensions than 'exp' except the
+'member_dim' dimension.}
+
+\item{time_dim}{A character string indicating the name of the time dimension.
+The default value is 'year'.}
+
+\item{member_dim}{A character string indicating the name of the member 
+dimension. The default value is 'member'.}
+
+\item{na.rm}{A logical value indicating whether to remove NA values during
+the computation. The default value is FALSE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+An array of the ratio of the predictable components. it has the same
+ dimensions as 'exp' except 'time_dim' and 'member_dim' dimensions.
+}
+\description{
+This function computes the ratio of predictable components (RPC; Eade et al., 2014).
+}
+\examples{
+exp <- array(data = runif(600), dim = c(year = 15, member = 10, lat = 2, lon = 2))
+obs <- array(data = runif(60), dim = c(year = 15, lat = 2, lon = 2))
+RatioPredictableComponents(exp, obs)
+
+}
diff --git a/man/RatioSDRMS.Rd b/man/RatioSDRMS.Rd
index 7dbd68283599edcfa599768a1946b00ce89fdea1..f1f6f3ddea92b777910d2fdc0176e9fd02a5f67a 100644
--- a/man/RatioSDRMS.Rd
+++ b/man/RatioSDRMS.Rd
@@ -67,11 +67,10 @@ rsdrms <- RatioSDRMS(sampleData$mod, sampleData$obs)
 rsdrms_plot <- array(dim = c(dim(rsdrms$ratio)[1:2], 4, dim(rsdrms$ratio)[3]))
 rsdrms_plot[, , 2, ] <- rsdrms$ratio
 rsdrms_plot[, , 4, ] <- rsdrms$p.val
-\donttest{
+\dontrun{
 PlotVsLTime(rsdrms_plot, toptitle = "Ratio ensemble spread / RMSE", ytitle = "", 
            monini = 11, limits = c(-1, 1.3), listexp = c('CMIP5 IC3'), 
-           listobs = c('ERSST'), biglab = FALSE, siglev = TRUE, 
-           fileout = 'tos_rsdrms.eps')
+           listobs = c('ERSST'), biglab = FALSE, siglev = TRUE) 
 }
 
 }
diff --git a/man/SignalNoiseRatio.Rd b/man/SignalNoiseRatio.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..f92cf0a1bff6262e0e886c47f153c9f22cfa9989
--- /dev/null
+++ b/man/SignalNoiseRatio.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/SignalNoiseRatio.R
+\name{SignalNoiseRatio}
+\alias{SignalNoiseRatio}
+\title{Calculate Signal-to-noise ratio}
+\usage{
+SignalNoiseRatio(
+  data,
+  time_dim = "year",
+  member_dim = "member",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{A numerical array with, at least, 'time_dim' and 'member_dim' 
+dimensions.}
+
+\item{time_dim}{A character string indicating the name of the time dimension
+in 'data'. The default value is 'year'.}
+
+\item{member_dim}{A character string indicating the name of the member 
+dimension in 'data'. The default value is 'member'.}
+
+\item{na.rm}{A logical value indicating whether to remove NA values during the
+computation. The default value is FALSE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+An array with of the signal-to-noise ratio. It has the same dimensions
+ as 'data' except 'time_dim' and 'member_dim' dimensions.
+}
+\description{
+This function computes the signal-to-noise ratio, where the signal is the 
+ensemble mean variance and the noise is the variance of the ensemble members 
+about the ensemble mean (Eade et al., 2014; Scaife and Smith, 2018).
+}
+\examples{
+exp <- array(data = runif(600), dim = c(year = 15, member = 10, lat = 2, lon = 2))
+SignalNoiseRatio(exp)
+
+}
diff --git a/man/Smoothing.Rd b/man/Smoothing.Rd
index 8d4a55871654d6159691f896bbed85434fe94da4..5769aa2b3763834651398bbecbd1268236264967 100644
--- a/man/Smoothing.Rd
+++ b/man/Smoothing.Rd
@@ -37,9 +37,8 @@ smooth_ano_exp <- Smoothing(ano_exp, time_dim = 'ftime', runmeanlen = 12)
 smooth_ano_obs <- Smoothing(ano_obs, time_dim = 'ftime', runmeanlen = 12)
 smooth_ano_exp <- Reorder(smooth_ano_exp, c(2, 3, 4, 1))
 smooth_ano_obs <- Reorder(smooth_ano_obs, c(2, 3, 4, 1))
- \donttest{
+ \dontrun{
 PlotAno(smooth_ano_exp, smooth_ano_obs, startDates, 
-       toptitle = "Smoothed Mediterranean mean SST", ytitle = "K",
-       fileout = "tos_smoothed_ano.png")
+       toptitle = "Smoothed Mediterranean mean SST", ytitle = "K")
  }
 }
diff --git a/man/Spread.Rd b/man/Spread.Rd
index 26e289e919178e5b945a2477e0f266992c7fcf08..e26bc14551fcbdbdfa5eaad2a9628119c9b65210 100644
--- a/man/Spread.Rd
+++ b/man/Spread.Rd
@@ -74,27 +74,27 @@ smooth_ano_exp_m_sub <- smooth_ano_exp - InsertDim(MeanDims(smooth_ano_exp, 'mem
                                                   name = 'member')
 spread <- Spread(smooth_ano_exp_m_sub, compute_dim = c('member', 'sdate'))
 
-\donttest{
+\dontrun{
 PlotVsLTime(Reorder(spread$iqr, c('dataset', 'stats', 'ftime')), 
            toptitle = "Inter-Quartile Range between ensemble members",
            ytitle = "K", monini = 11, limits = NULL, 
            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-           hlines = c(0), fileout = 'tos_iqr.png')
+           hlines = c(0))
 PlotVsLTime(Reorder(spread$maxmin, c('dataset', 'stats', 'ftime')), 
            toptitle = "Maximum minus minimum of the members", 
            ytitle = "K", monini = 11, limits = NULL, 
            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-           hlines = c(0), fileout = 'tos_maxmin.png')
+           hlines = c(0))
 PlotVsLTime(Reorder(spread$sd, c('dataset', 'stats', 'ftime')), 
            toptitle = "Standard deviation of the members", 
            ytitle = "K", monini = 11, limits = NULL, 
            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-           hlines = c(0), fileout = 'tos_sd.png')
+           hlines = c(0))
 PlotVsLTime(Reorder(spread$mad, c('dataset', 'stats', 'ftime')), 
            toptitle = "Median Absolute Deviation of the members",
            ytitle = "K", monini = 11, limits = NULL, 
            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-           hlines = c(0), fileout = 'tos_mad.png')
+           hlines = c(0))
 }
 
 }
diff --git a/s2dv-manual.pdf b/s2dv-manual.pdf
index de5dc4cb860f84e21a44781861245233d13c7fdf..b6faa89fbc92c8735e85d4fa7bc4d0f68d47030b 100644
Binary files a/s2dv-manual.pdf and b/s2dv-manual.pdf differ
diff --git a/tests/testdata/testIrregularGrid_1.rds b/tests/testdata/testIrregularGrid_1.rds
new file mode 100644
index 0000000000000000000000000000000000000000..21a46a726737566eee7880c86aacb45d7fa37ba2
Binary files /dev/null and b/tests/testdata/testIrregularGrid_1.rds differ
diff --git a/tests/testdata/testIrregularGrid_2.rds b/tests/testdata/testIrregularGrid_2.rds
new file mode 100644
index 0000000000000000000000000000000000000000..d014eafdf520077459c3bd6e49948ed7608c2eba
Binary files /dev/null and b/tests/testdata/testIrregularGrid_2.rds differ
diff --git a/tests/testdata/testIrregularGrid_3.rds b/tests/testdata/testIrregularGrid_3.rds
new file mode 100644
index 0000000000000000000000000000000000000000..3dc6acb70c4f69ece6d812992335391289a1c3f7
Binary files /dev/null and b/tests/testdata/testIrregularGrid_3.rds differ
diff --git a/tests/testthat/test-AMV.R b/tests/testthat/test-AMV.R
index 562605b6cf04efc348c27456fdbbc63250306b7b..9adfaefa51211858c4b475b5446b1ea4ef0298c2 100644
--- a/tests/testthat/test-AMV.R
+++ b/tests/testthat/test-AMV.R
@@ -134,16 +134,6 @@ test_that("1. Input checks", {
       month_dim = 'mth'),
   "Parameter 'month_dim' is not found in 'data' dimension."
   )
-  expect_error(
-  AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist',
-      member_dim = 1),
-  "Parameter 'member_dim' must be a character string."
-  )
-  expect_error(
-  AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist',
-      member_dim = 'ens'),
-  "Parameter 'member_dim' is not found in 'data' dimension."
-  )
 
 })
 
@@ -179,7 +169,7 @@ test_that("3. Output checks: dat2", {
  )
  expect_equal(
  mean(res2_1)*10^17,
- -3.549383,
+ -4.884981,
  tolerance = 0.00001
  ) 
 
@@ -187,7 +177,7 @@ test_that("3. Output checks: dat2", {
                indices_for_clim = 2:3)
  expect_equal(
  mean(res2_2)*10^16,
- 8.616545,
+ 8.837375,
  tolerance = 0.00001
  )
 
@@ -195,7 +185,7 @@ test_that("3. Output checks: dat2", {
                mask = mask)
  expect_equal(
  mean(res2_3)*10^17,
- -4.437162,
+ -1.554312,
  tolerance = 0.00001
  )
 
@@ -233,7 +223,7 @@ test_that("4. Output checks: dat3", {
                mask = mask)
  expect_equal(
  mean(res3_4)*10^17,
- 2.220446,
+ -1.554312,
  tolerance = 0.00001
  )
 
diff --git a/tests/testthat/test-Ano_CrossValid.R b/tests/testthat/test-Ano_CrossValid.R
index 60bff8c70a1752c795664cc2acebed6dc9c81413..b66fc5fbf54113b71826e7adf79d014890c003e0 100644
--- a/tests/testthat/test-Ano_CrossValid.R
+++ b/tests/testthat/test-Ano_CrossValid.R
@@ -79,8 +79,7 @@ test_that("1. Input checks", {
   # exp and obs (2)
   expect_error(
     Ano_CrossValid(exp1, array(1:20, dim = c(dataset = 1, member = 2, sdate = 4, ftime = 2))),
-    paste0("Parameter 'exp' and 'obs' must have same length of ",
-                "all dimension expect 'dat_dim'.")
+    "Parameter 'exp' and 'obs' must have the same length of all dimensions except 'dat_dim'."
   )
 
 
diff --git a/tests/testthat/test-CDORemap.R b/tests/testthat/test-CDORemap.R
new file mode 100644
index 0000000000000000000000000000000000000000..16da1b090a39a2dfa5659cfabb5c8c2af57aa62e
--- /dev/null
+++ b/tests/testthat/test-CDORemap.R
@@ -0,0 +1,271 @@
+context("s2dv::CDORemap tests")
+
+# data1: regular grid
+data1 <- array(1:360*181*2, dim = c(lon = 360, lat = 181, time = 2))
+lons1 <- seq(0, 359)
+lats1 <- seq(-90, 90)
+data1_1 <- .aperm2(data1, c(3, 1, 2))
+data1_2 <- InsertDim(data1, 1, 1, name = 'var')
+data1_3 <- InsertDim(data1_2, 5, 1, name = 'dat')
+
+# data2: irregular grid
+data2 <- readRDS('../testdata/testIrregularGrid_1.rds')
+lons2 <- attributes(data2)$Variables$common$nav_lon#[2:361, 2:331]
+lats2 <- attributes(data2)$Variables$common$nav_lat#[2:361, 2:331]
+
+data2_1 <- ClimProjDiags::Subset(data2, 1, 1, drop = 'selected')
+data2_2 <- ClimProjDiags::Subset(data2, 5, 1, drop = 'selected')
+data2_3 <- ClimProjDiags::Subset(data2, c(1, 2), list(1, 1), drop = 'selected')
+data2_4 <- ClimProjDiags::Subset(data2, c(1, 2, 5), list(1, 1, 1), drop = 'selected')
+data2_5 <- .aperm2(data2, c(3,4,1,2,5))
+data2_6 <- .aperm2(data2, c(3,4,5,2,1))
+data2_7 <- .aperm2(data2, c(4,3,1,2,5))
+data2_8 <- .aperm2(data2_4, c(2,1))
+
+# data3: irregular grid
+data3 <- readRDS('../testdata/testIrregularGrid_2.rds')
+lons3 <- attributes(data3)$Variables$common$longitude
+lats3 <- attributes(data3)$Variables$common$latitude
+
+data3_1 <- drop(data3)
+data3_2 <- .aperm2(data3_1, c(2, 1))
+
+# data4: irregular grid
+data4 <- readRDS('../testdata/testIrregularGrid_3.rds')
+lons4 <- attributes(data4)$Variables$common$longitude
+lats4 <- attributes(data4)$Variables$common$latitude
+
+data4_1 <- drop(data4)
+data4_2 <- ClimProjDiags::Subset(data4, c(1,2,3,4,8), list(1,1,1,1,1), drop = 'selected')
+data4_3 <- .aperm2(data4_2, c(1, 3, 2))
+
+##############################################
+
+test_that("1. Input checks", {
+
+expect_error(
+CDORemap(data_array = NULL, lons = 'lons', lats = 'lats'),
+"Expected numeric 'lons' and 'lats'."
+)
+expect_error(
+CDORemap(data_array = NULL, lons = c(NA, lons1), lats = lats1),
+"Found invalid values in 'lons'."
+)
+expect_error(
+CDORemap(data_array = NULL, lons = lons1, lats = c(NA, lats1)),
+"Found invalid values in 'lats'."
+)
+
+})
+
+test_that("2. regular regrid", {
+suppressWarnings(
+res <- CDORemap(data1, lons1, lats1, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res1 <- CDORemap(data1_1, lons1, lats1, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res2 <- CDORemap(data1_2, lons1, lats1, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res3 <- CDORemap(data1_3, lons1, lats1, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+
+expect_equal(
+as.vector(res),
+as.vector(.aperm2(res1, c(2,3,1)))
+)
+expect_equal(
+as.vector(res),
+as.vector(res2)
+)
+expect_equal(
+as.vector(res),
+as.vector(res3)
+)
+
+})
+
+##################################################################
+
+test_that("3. dat2: irregular regrid", {
+suppressWarnings(
+res <- CDORemap(data2, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res1 <- CDORemap(data2_1, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res2 <- CDORemap(data2_2, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res3 <- CDORemap(data2_3, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res4 <- CDORemap(data2_4, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res5 <- CDORemap(data2_5, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res6 <- CDORemap(data2_6, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res7 <- CDORemap(data2_7, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res8 <- CDORemap(data2_8, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+# no data
+suppressWarnings(
+res9 <- CDORemap(NULL, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)
+)
+
+expect_equal(
+drop(res)[95:100, 50],
+c(99.54009, 99.53733, 99.53453, 99.53169, 99.52792, 99.52304),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(res),
+as.vector(res1)
+)
+expect_equal(
+as.vector(res),
+as.vector(res2)
+)
+expect_equal(
+as.vector(res),
+as.vector(res3)
+)
+expect_equal(
+as.vector(res),
+as.vector(res4)
+)
+expect_equal(
+as.vector(res),
+as.vector(res5)
+)
+expect_equal(
+as.vector(res),
+as.vector(res6)
+)
+expect_equal(
+as.vector(res),
+as.vector(res7)
+)
+expect_equal(
+as.vector(res),
+as.vector(res8)
+)
+expect_equal(
+res9$data_array,
+NULL
+)
+expect_equal(
+as.vector(res9$lons)[1:5],
+c(0.0, 3.6, 7.2, 10.8, 14.4)
+)
+expect_equal(
+as.vector(res9$lats)[1:5],
+c(-88.2, -84.6, -81.0, -77.4, -73.8)
+)
+
+})
+
+#########################################################
+
+test_that("4. dat3: irregular regrid", {
+suppressWarnings(
+res <- CDORemap(data3, lons3, lats3, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res1 <- CDORemap(data3_1, lons3, lats3, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res2 <- CDORemap(data3_2, lons3, lats3, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+# no data
+suppressWarnings(
+res3 <- CDORemap(NULL, lons3, lats3, grid = 'r100x50', method = 'bil', crop = F)
+)
+
+expect_equal(
+as.vector(res),
+as.vector(res1)
+)
+expect_equal(
+as.vector(res),
+as.vector(res2)
+)
+expect_equal(
+res1[1:5, 50],
+c(-1.588068, -1.586225, -1.584405, -1.582579, -1.580713),
+tolerance = 0.0001
+)
+expect_equal(
+res3$data_array,
+NULL
+)
+expect_equal(
+as.vector(res3$lons)[1:5],
+c(0.0, 3.6, 7.2, 10.8, 14.4)
+)
+expect_equal(
+as.vector(res3$lats)[1:5],
+c(-88.2, -84.6, -81.0, -77.4, -73.8)
+)
+
+})
+
+############################################################
+
+test_that("5. dat4: irregular regrid", {
+suppressWarnings(
+res <- CDORemap(data4, lons4, lats4, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res1 <- CDORemap(data4_1, lons4, lats4, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res2 <- CDORemap(data4_2, lons4, lats4, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res3 <- CDORemap(data4_3, lons4, lats4, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+# no data
+suppressWarnings(
+res4 <- CDORemap(NULL, lons4, lats4, grid = 'r100x50', method = 'bil', crop = F)
+)
+expect_equal(
+as.vector(res),
+as.vector(res1)
+)
+expect_equal(
+as.vector(res),
+as.vector(res2)
+)
+expect_equal(
+as.vector(res),
+as.vector(res3)
+)
+expect_equal(
+res1[1:5, 50],
+c(-1.795211, -1.795057, -1.794898, -1.794707, -1.794533),
+tolerance = 0.0001
+)
+expect_equal(
+res4$data_array,
+NULL
+)
+expect_equal(
+as.vector(res4$lons)[1:5],
+c(0.0, 3.6, 7.2, 10.8, 14.4)
+)
+expect_equal(
+as.vector(res4$lats)[1:5],
+c(-88.2, -84.6, -81.0, -77.4, -73.8)
+)
+
+})
diff --git a/tests/testthat/test-Cluster.R b/tests/testthat/test-Cluster.R
index 9071cd8e5e479e0ef2ccee2b6b74208815e604c8..4e0a7b171d25ec559afbef8b97dafb45e87cff6e 100644
--- a/tests/testthat/test-Cluster.R
+++ b/tests/testthat/test-Cluster.R
@@ -61,17 +61,17 @@ test_that("1. Input checks", {
   )
   # nclusters
   expect_error(
-  Cluster(dat1, weights1, ncluster = 1),
+  Cluster(dat1, weights1, ncluster = 1, space_dim = 'space'),
   "Parameter 'nclusters' must be an integer bigger than 1."
   )
   # index
   expect_error(
-  Cluster(dat1, weights1, index = 1),
+  Cluster(dat1, weights1, index = 1, space_dim = 'space'),
   "Parameter 'index' should be a character strings accepted as 'index' by the function NbClust::NbClust."
   )
   # ncores
   expect_error(
-    Cluster(dat1, weights1, ncore = 0),
+    Cluster(dat1, weights1, ncore = 0, space_dim = 'space'),
     "Parameter 'ncores' must be a positive integer."
   )
 })
@@ -80,7 +80,7 @@ test_that("1. Input checks", {
 test_that("2. Output checks: dat1", {
 # The output is random. Only check dimensions.
   expect_equal(
-    length(Cluster(dat1, weights1)$cluster),
+    length(Cluster(dat1, weights1, space_dim = 'space')$cluster),
     50
   )
   expect_equal(
@@ -88,11 +88,11 @@ test_that("2. Output checks: dat1", {
     100
   )
   expect_equal(
-    dim(Cluster(dat1, weights1)$centers),
+    dim(Cluster(dat1, weights1, space_dim = 'space')$centers),
     c(8, 2)
   )
   expect_equal(
-    dim(Cluster(dat1, weights1, nclusters = 3)$centers),
+    dim(Cluster(dat1, weights1, nclusters = 3, space_dim = 'space')$centers),
     c(3, 2)
   )
 
@@ -103,7 +103,7 @@ test_that("2. Output checks: dat1", {
 test_that("3. Output checks: dat2", {
 
   expect_equal(
-    length(Cluster(dat2, weights2)$cluster),
+    length(Cluster(dat2, weights2, space_dim = c('lat', 'lon'))$cluster),
     50
   )
   expect_equal(
@@ -115,11 +115,11 @@ test_that("3. Output checks: dat2", {
     50
   )
   expect_equal(
-    dim(Cluster(dat2, weights2)$centers),
+    dim(Cluster(dat2, weights2, space_dim = c('lat', 'lon'))$centers),
     c(7, 6)
   )
   expect_equal(
-    dim(Cluster(dat2, weights2, nclusters = 5)$centers),
+    dim(Cluster(dat2, weights2, nclusters = 5, space_dim = c('lat', 'lon'))$centers),
     c(5, 6)
   )
 
diff --git a/tests/testthat/test-Consist_Trend.R b/tests/testthat/test-Consist_Trend.R
index 2dd2214a8ea31e851fae8564921c66266d54342c..aa66f45761c6aa6b8bca1f02d7fbf8649882b038 100644
--- a/tests/testthat/test-Consist_Trend.R
+++ b/tests/testthat/test-Consist_Trend.R
@@ -151,8 +151,10 @@ c(NA, -0.4826962, 1.2716524, -1.0952163, 0.3062600),
 tolerance = 0.0001
 )
 expect_equal(
-mean(Consist_Trend(exp2, obs2)$detrended_obs, na.rm = TRUE)*10^18,
-2.118364,
+as.vector(Consist_Trend(exp2, obs2)$detrended_obs),
+c(NA, -0.4826962, 1.2716524, -1.0952163, 0.3062600, -0.2100316, 0.2920924,
+  -0.7289708, 1.4217907, -0.7748807, 0.6108713, NA, -0.5417854, -1.3599143,
+  1.2908284),
 tolerance = 0.0001
 )
 expect_equal(
diff --git a/tests/testthat/test-EOF.R b/tests/testthat/test-EOF.R
index 7966518dec3472f37d267f36c9a9589ac089c43d..01428e6eb73bd736c9ca7b7b6c5da83fa6c4f678 100644
--- a/tests/testthat/test-EOF.R
+++ b/tests/testthat/test-EOF.R
@@ -205,20 +205,20 @@ test_that("2. dat1", {
 ##############################################
 test_that("3. dat2", {
   expect_equal(
-  dim(EOF(dat2, lon = lon2, lat = lat2)$EOFs),
+  dim(EOF(dat2, lon = lon2, lat = lat2, neofs = 12)$EOFs),
   c(mode = 12, lat = 6, lon = 2)
   )
   expect_equal(
-  dim(EOF(dat2, lon = lon2, lat = lat2)$PCs),
+  dim(EOF(dat2, lon = lon2, lat = lat2, neofs = 12)$PCs),
   c(sdate = 20, mode = 12)
   )
   expect_equal(
-  EOF(dat2, lon = lon2, lat = lat2)$EOFs[1:5],
+  EOF(dat2, lon = lon2, lat = lat2, neofs = 12)$EOFs[1:5],
   c(0.33197201, 0.18837900, -0.19697143, 0.08305805, -0.51297585),
   tolerance = 0.0001
   )
   expect_equal(
-  mean(EOF(dat2, lon = lon2, lat = lat2)$EOFs),
+  mean(EOF(dat2, lon = lon2, lat = lat2, neofs = 12)$EOFs),
   0.02720393,
   tolerance = 0.0001
   )
@@ -228,37 +228,37 @@ test_that("3. dat2", {
 ##############################################
 test_that("4. dat3", {
   expect_equal(
-  dim(EOF(dat3, lon = lon3, lat = lat3)$EOFs),
+  dim(EOF(dat3, lon = lon3, lat = lat3, neofs = 12)$EOFs),
   c(mode = 12, lat = 6, lon = 2, dat = 2)
   )
   expect_equal(
-  dim(EOF(dat3, lon = lon3, lat = lat3)$PCs),
+  dim(EOF(dat3, lon = lon3, lat = lat3, neofs = 12)$PCs),
   c(sdate = 20, mode = 12, dat = 2)
   )
   expect_equal(
-  dim(EOF(dat3, lon = lon3, lat = lat3)$var),
+  dim(EOF(dat3, lon = lon3, lat = lat3, neofs = 12)$var),
   c(mode = 12, dat = 2)
   )
   expect_equal(
-  dim(EOF(dat3, lon = lon3, lat = lat3)$mask),
+  dim(EOF(dat3, lon = lon3, lat = lat3, neofs = 12)$mask),
   c(lat = 6, lon = 2, dat = 2)
   )
   expect_equal(
-  dim(EOF(dat3, lon = lon3, lat = lat3)$wght),
+  dim(EOF(dat3, lon = lon3, lat = lat3, neofs = 12)$wght),
   c(lat = 6, lon = 2)
   )
   expect_equal(
-  mean(EOF(dat3, lon = lon3, lat = lat3)$EOFs),
+  mean(EOF(dat3, lon = lon3, lat = lat3, neofs = 12)$EOFs),
   0.01214845,
   tolerance = 0.0001
   )
   expect_equal(
-  EOF(dat3, lon = lon3, lat = lat3)$EOFs[1:5],
+  EOF(dat3, lon = lon3, lat = lat3, neofs = 12)$EOFs[1:5],
   c(0.3292733, 0.1787016, -0.3801986, 0.1957160, -0.4377031),
   tolerance = 0.0001
   )
   expect_equal(
-  EOF(dat3, lon = lon3, lat = lat3)$tot_var,
+  EOF(dat3, lon = lon3, lat = lat3, neofs = 12)$tot_var,
   array(c(213.2422, 224.4203), dim = c(dat = 2)),
   tolerance = 0.0001
   )
diff --git a/tests/testthat/test-GMST.R b/tests/testthat/test-GMST.R
index d6e3a4f882f02d38fa5d22164b7b7fc66398191b..01ab792c26a2145d7794659a81922d4de933ef39 100644
--- a/tests/testthat/test-GMST.R
+++ b/tests/testthat/test-GMST.R
@@ -171,18 +171,6 @@ test_that("1. Input checks", {
        month_dim = 'mth'),
   "Parameter 'month_dim' is not found in 'data_tas' or 'data_tos' dimension."
   )
-  expect_error(
-  GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon,
-       mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist',
-       member_dim = 1),
-  "Parameter 'member_dim' must be a character string."
-  )
-  expect_error(
-  GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon,
-       mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist',
-       member_dim = 'ens'),
-  "Parameter 'member_dim' is not found in 'data_tas' or 'data_tos' dimension."
-  )
 
 })
 
@@ -224,7 +212,7 @@ test_that("3. Output checks: dat2", {
  ) 
  expect_equal(
  range(res2_1),
- c(-2.028458, 2.028458),
+ c(-2, 2),
  tolerance = 0.00001
  )
 
@@ -241,7 +229,7 @@ test_that("3. Output checks: dat2", {
                 mask = mask)
  expect_equal(
  mean(res2_3)*10^16,
- 2.842084,
+ 0,
  tolerance = 0.00001
  )
 
@@ -281,7 +269,7 @@ test_that("4. Output checks: dat3", {
                 mask = mask)
  expect_equal(
  mean(res3_4)*10^16,
- -8.52686,
+ -5.684949,
  tolerance = 0.00001
  )
 
diff --git a/tests/testthat/test-GSAT.R b/tests/testthat/test-GSAT.R
index bd53e8ab062e232098cdcdb0d0db4e33dd912919..2d7d7e00c20f20d16288208e9ce72f12e2fbe454 100644
--- a/tests/testthat/test-GSAT.R
+++ b/tests/testthat/test-GSAT.R
@@ -134,16 +134,6 @@ test_that("1. Input checks", {
       month_dim = 'mth'),
   "Parameter 'month_dim' is not found in 'data' dimension."
   )
-  expect_error(
-  GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist',
-      member_dim = 1),
-  "Parameter 'member_dim' must be a character string."
-  )
-  expect_error(
-  GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist',
-      member_dim = 'ens'),
-  "Parameter 'member_dim' is not found in 'data' dimension."
-  )
 
 })
 
@@ -179,7 +169,7 @@ test_that("3. Output checks: dat2", {
  )
  expect_equal(
  mean(res2_1)*10^16,
- 2.841867,
+ 2.841954,
  tolerance = 0.00001
  ) 
 
@@ -195,7 +185,7 @@ test_that("3. Output checks: dat2", {
                mask = mask)
  expect_equal(
  mean(res2_3)*10^15,
- 1.136842,
+ 0.2841954,
  tolerance = 0.00001
  )
 
diff --git a/tests/testthat/test-InsertDim.R b/tests/testthat/test-InsertDim.R
index f3cfab801f33fd29e44e805ae65d6f8c52d9dad7..876e7e36c2e0c74f0a4a2017d35f4857e2109816 100644
--- a/tests/testthat/test-InsertDim.R
+++ b/tests/testthat/test-InsertDim.R
@@ -56,20 +56,24 @@ test_that("1. Input checks", {
 ##############################################
 test_that("2. Output checks: dat1", {
 
+  expect_warning(
+    InsertDim(dat1, posdim = 1, lendim = 2),
+    "The name of new dimension is not given. Set the name as 'new'."
+  )
   expect_equal(
-    dim(InsertDim(dat1, posdim = 1, lendim = 2)),
+    dim(InsertDim(dat1, posdim = 1, lendim = 2, name = 'new')),
     dim(array(dim = c(new = 2, dat = 1, sdate = 13, ftime = 2)))
   )
   expect_equal(
-    dim(InsertDim(dat1, posdim = 3, lendim = c(d = 2))),
+    dim(InsertDim(dat1, posdim = 3, lendim = 2, name = 'd')),
     c(dat = 1, sdate = 13, d = 2, ftime = 2)
   )
   expect_equal(
-    as.vector(InsertDim(dat1, posdim = 1, lendim = 2)[1,,,]),
+    as.vector(InsertDim(dat1, posdim = 1, lendim = 2, name = 'd')[1,,,]),
     as.vector(dat1)
   )
   expect_equal(
-    as.vector(InsertDim(dat1, posdim = 1, lendim = 2)[2,,,]),
+    as.vector(InsertDim(dat1, posdim = 1, lendim = 2, name = 'd')[2,,,]),
     as.vector(dat1)
   )
 
@@ -84,13 +88,13 @@ test_that("3. Output checks: dat2", {
   )
 
   expect_equal(
-    as.vector(InsertDim(dat2, posdim = 3, lendim = 1)[,,1,]),
+    as.vector(InsertDim(dat2, posdim = 3, lendim = 1, name = 'd')[,,1,]),
     as.vector(dat2)
   )
 
   expect_equal(
-    dim(InsertDim(InsertDim(dat2, posdim = 4, lendim = 2), posdim = 1, lendim = 4)),
-    dim(array(dim = c(new = 4, 2, 3, c = 4, new = 2)))
+    dim(InsertDim(InsertDim(dat2, posdim = 4, lendim = 2, name = 'd'), posdim = 1, lendim = 4, name = 'e')),
+    dim(array(dim = c(e = 4, 2, 3, c = 4, d = 2)))
   )
 
 })
diff --git a/tests/testthat/test-MeanDims.R b/tests/testthat/test-MeanDims.R
index 9c7c5666a6246759ea7abed3f478bd85f8749f2a..a431b995422d6cff9050d5f7a2248da0f15451bd 100644
--- a/tests/testthat/test-MeanDims.R
+++ b/tests/testthat/test-MeanDims.R
@@ -57,7 +57,10 @@ test_that("1. Input checks", {
   MeanDims(dat1, dims = 'ftime', na.rm = na.omit),
   "Parameter 'na.rm' must be one logical value."
   )
-
+  expect_error(
+  MeanDims(dat1, dims = 'ftime', drop = 'selected'),
+  "Parameter 'drop' must be one logical value."
+  )
 })
 
 ##############################################
@@ -79,6 +82,26 @@ test_that("2. Output checks: dat1", {
   MeanDims(dat1, dims = c('sdate'))[1:2],
   c(3, 8)
   )
+  expect_equal(
+  dim(MeanDims(dat1, dims = 1:2, drop = F)),
+  c(dat = 1, sdate = 1, ftime = 4)
+  )
+  expect_equal(
+  as.vector(drop(MeanDims(dat1, dims = 1:2, drop = F))),
+  as.vector(MeanDims(dat1, dims = 1:2, drop = T))
+  )
+  expect_equal(
+  dim(MeanDims(dat1, dims = 1:3, drop = F)),
+  c(dat = 1, sdate = 1, ftime = 1)
+  )
+  expect_equal(
+  dim(MeanDims(dat1, dims = 1:3, drop = T)),
+  1
+  )
+  expect_equal(
+  as.vector(drop(MeanDims(dat1, dims = 1:3, drop = F))),
+  as.vector(MeanDims(dat1, dims = 1:3, drop = T))
+  )
 
 })
 
@@ -126,6 +149,10 @@ test_that("5. Output checks: dat4", {
     length(MeanDims(dat4, dims = 1)),
     1
   )
+  expect_equal(
+    dim(MeanDims(dat4, dims = 1, drop = F)),
+    1
+  )
 })
 ##############################################
 
diff --git a/tests/testthat/test-Persistence.R b/tests/testthat/test-Persistence.R
index 2ce2dec47fe1872a183c46338b3fbcf946bef5b5..4eabe836734c3ecf8cd636facb8a9e4b69ffb6ff 100644
--- a/tests/testthat/test-Persistence.R
+++ b/tests/testthat/test-Persistence.R
@@ -8,7 +8,7 @@ dim(dat1) <- c(member = 1, time = 70, lat = 6, lon = 7)
 dates1 <- seq(1920, 1989, 1)
 start1 <- 1961
 end1 <- 1990
-res <- Persistence(obs1, dates = dates1, start = 1961, end = 1990, ft_start = 1,
+res <- Persistence(dat1, dates = dates1, start = 1961, end = 1990, ft_start = 1,
                    nmemb = 40)
 
 #dat2: day
diff --git a/tests/testthat/test-REOF.R b/tests/testthat/test-REOF.R
index 9f4bb480a5dc5bb3de121041fa9f72afe1cea72c..e118e1eac8284161c18fabdeee51c4ccf827101a 100644
--- a/tests/testthat/test-REOF.R
+++ b/tests/testthat/test-REOF.R
@@ -92,13 +92,17 @@ test_that("1. Input checks", {
 
 ##############################################
 test_that("2. dat1", {
-
+  expect_warning(
+  REOF(dat1, lon = lon1, lat = lat1),
+  "Parameter 'ntrunc' is changed to 10, the minimum among the length of time_dim, the production of the length of space_dim, and ntrunc.",
+  fixed = TRUE
+  )
   expect_equal(
-  names(REOF(dat1, lon = lon1, lat = lat1, ntrunc = 5)),
+  names(REOF(dat1, lon = lon1, lat = lat1, ntrunc = 10)),
   c("REOFs", "RPCs", "var", "wght")
   )
   expect_equal(
-  dim(REOF(dat1, lon = lon1, lat = lat1)$REOFs),
+  dim(REOF(dat1, lon = lon1, lat = lat1, ntrunc = 10)$REOFs),
   c(mode = 10, lat = 6, lon = 2)
   )
   expect_equal(
@@ -147,7 +151,7 @@ test_that("2. dat1", {
 ##############################################
 test_that("3. dat2", {
   expect_equal(
-  dim(REOF(dat2, lon = lon2, lat = lat2)$REOFs),
+  dim(REOF(dat2, lon = lon2, lat = lat2, ntrunc = 5)$REOFs),
   c(mode = 5, lat = 6, lon = 2, dat = 2)
   )
   expect_equal(
@@ -164,7 +168,7 @@ test_that("3. dat2", {
   tolerance = 0.0001
   )
   expect_equal(
-  mean(REOF(dat2, lon = lon2, lat = lat2)$REOFs),
+  mean(REOF(dat2, lon = lon2, lat = lat2, ntrunc = 5)$REOFs),
   0.01120786,
   tolerance = 0.0001
   )
diff --git a/tests/testthat/test-RatioPredictableComponents.R b/tests/testthat/test-RatioPredictableComponents.R
new file mode 100644
index 0000000000000000000000000000000000000000..b6a080835fa7fbfd2ccd51c756c2c1b9124cd1fb
--- /dev/null
+++ b/tests/testthat/test-RatioPredictableComponents.R
@@ -0,0 +1,99 @@
+context("RatioPredictableComponents test")
+
+##############################################
+  # dat1
+  set.seed(1)
+  exp1 <- array(rnorm(12), c(dat = 1, member = 2, year = 3, time = 2)) 
+  set.seed(2)
+  obs1 <- array(rnorm(6), c(dat = 1, year = 3, time = 2))
+
+  # dat2
+  exp2 <- exp1
+  exp2[1, 1, 2, 1] <- NA
+  obs2 <- obs1
+
+##############################################
+
+test_that("1. Input checks", {
+
+  expect_error(
+  RatioPredictableComponents(c()),
+  "Parameter 'exp' cannot be NULL."
+  )
+  expect_error(
+    RatioPredictableComponents(c(NA, NA)),
+    "Parameter 'exp' must be a numeric array."
+  )
+  expect_error(
+  RatioPredictableComponents(exp1, c()),
+  "Parameter 'obs' cannot be NULL."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, c(NA, NA)),
+    "Parameter 'obs' must be a numeric array."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, obs1, time_dim = 'a'),
+    "'exp' must have 'time_dim' dimension."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, obs1, member_dim = 'ens'),
+    "'exp' must have 'member_dim' dimension."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, array(rnorm(6), dim = c(sdate = 3, time = 2))),
+    "'obs' must have 'time_dim' dimension."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, array(rnorm(4), dim = c(year = 2, time = 2))),
+    "'exp' and 'obs' must have the same length for 'time_dim'."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, obs1, na.rm = 'TRUE'),
+    "Parameter 'na.rm' must be TRUE or FALSE."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, obs1, ncore = 0),
+    "Parameter 'ncores' must be a positive integer."
+  )
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+expect_equal(
+dim(RatioPredictableComponents(exp1, obs1)),
+c(dat = 1, time = 2)
+)
+expect_equal(
+dim(RatioPredictableComponents(exp1, obs1, time_dim = 'time')),
+c(dat = 1, year = 3)
+)
+expect_equal(
+as.vector(RatioPredictableComponents(exp1, obs1)),
+c(-0.29462414, 0.07955797),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RatioPredictableComponents(exp1, obs1, time_dim = 'time')),
+c(-1.054665, 6.843041, -1.000069),
+tolerance = 0.0001
+)
+
+})
+
+test_that("2. Output checks: dat2", {
+
+expect_equal(
+as.vector(RatioPredictableComponents(exp2, obs2, na.rm = T)),
+c(-0.07981636, 0.07955797),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RatioPredictableComponents(exp2, obs2, na.rm = F)),
+c(NA, 0.07955797),
+tolerance = 0.0001
+)
+})
+
diff --git a/tests/testthat/test-RatioRMS.R b/tests/testthat/test-RatioRMS.R
index b70d6fb5404998d7d737579a57012dd378383317..11df46ccf04d14e305059163e5c5d670906e57bc 100644
--- a/tests/testthat/test-RatioRMS.R
+++ b/tests/testthat/test-RatioRMS.R
@@ -114,12 +114,12 @@ names(RatioRMS(exp2_1, exp2_2, obs2)),
 c('ratiorms', 'p.val')
 )
 expect_equal(
-RatioRMS(exp2_1, exp2_2, obs2)$p.val,
+as.vector(RatioRMS(exp2_1, exp2_2, obs2)$p.val),
 0.7418331,
 tolerance = 0.0001
 )
 expect_equal(
-RatioRMS(exp2_1, exp2_2, obs2)$ratiorms,
+as.vector(RatioRMS(exp2_1, exp2_2, obs2)$ratiorms),
 0.8931399,
 tolerance = 0.0001
 )
diff --git a/tests/testthat/test-SPOD.R b/tests/testthat/test-SPOD.R
index 6d23a59745158a7d3c36807244cf96ceb453608c..8b56c721dc618409e5c14084262eb71aa07566af 100644
--- a/tests/testthat/test-SPOD.R
+++ b/tests/testthat/test-SPOD.R
@@ -134,16 +134,6 @@ test_that("1. Input checks", {
       month_dim = 'mth'),
   "Parameter 'month_dim' is not found in 'data' dimension."
   )
-  expect_error(
-  SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist',
-      member_dim = 1),
-  "Parameter 'member_dim' must be a character string."
-  )
-  expect_error(
-  SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist',
-      member_dim = 'ens'),
-  "Parameter 'member_dim' is not found in 'data' dimension."
-  )
 
 })
 
@@ -179,15 +169,15 @@ test_that("3. Output checks: dat2", {
  )
  expect_equal(
  mean(res2_1)*10^16,
- 2.754588,
+ 2.087219,
  tolerance = 0.00001
  ) 
 
  res2_2 <- SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist',
                indices_for_clim = 2:3)
  expect_equal(
- mean(res2_2)*10^16,
- 2.645063,
+ mean(res2_2)*10^18,
+ 8.881784,
  tolerance = 0.00001
  )
 
@@ -195,7 +185,7 @@ test_that("3. Output checks: dat2", {
                mask = mask)
  expect_equal(
  mean(res2_3)*10^17,
- 5.3553,
+ -12.4345,
  tolerance = 0.00001
  )
 
@@ -233,7 +223,7 @@ test_that("4. Output checks: dat3", {
                mask = mask)
  expect_equal(
  mean(res3_4)*10^17,
- 7.549517,
+ -8.881784,
  tolerance = 0.00001
  )
 
diff --git a/tests/testthat/test-SignalNoiseRatio.R b/tests/testthat/test-SignalNoiseRatio.R
new file mode 100644
index 0000000000000000000000000000000000000000..9fc5ce35478b9988d35523c4f19be5562f4185cf
--- /dev/null
+++ b/tests/testthat/test-SignalNoiseRatio.R
@@ -0,0 +1,73 @@
+
+context("SignalNoiseRatio test")
+
+##############################################
+  # dat1
+  set.seed(1)
+  data1 <- array(rnorm(12), c(dat = 1, member = 2, year = 3, time = 2)) 
+
+  # dat2
+  data2 <- data1
+  data2[1, 1, 2, 1] <- NA
+
+##############################################
+
+test_that("1. Input checks", {
+
+  expect_error(
+  SignalNoiseRatio(c()),
+  "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    SignalNoiseRatio(c(NA, NA)),
+    "Parameter 'data' must be a numeric array."
+  )
+  expect_error(
+    SignalNoiseRatio(data1, time_dim = 'a'),
+    "'data' must have 'time_dim' dimension."
+  )
+  expect_error(
+    SignalNoiseRatio(data1, member_dim = 'ens'),
+    "'data' must have 'member_dim' dimension."
+  )
+  expect_error(
+    SignalNoiseRatio(data1, na.rm = 'TRUE'),
+    "Parameter 'na.rm' must be TRUE or FALSE."
+  )
+  expect_error(
+    SignalNoiseRatio(data1, ncore = 0),
+    "Parameter 'ncores' must be a positive integer."
+  )
+
+})
+
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+expect_equal(
+dim(SignalNoiseRatio(data1)),
+c(dat = 1, time = 2)
+)
+expect_equal(
+as.vector(SignalNoiseRatio(data1)),
+c(0.1591161, 0.8003926),
+tolerance = 0.0001
+)
+
+})
+
+test_that("2. Output checks: dat2", {
+expect_equal(
+as.vector(SignalNoiseRatio(data2, na.rm = T)),
+c(4.5075526, 0.8003926),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(SignalNoiseRatio(data2, na.rm = F)),
+c(NA, 0.8003926),
+tolerance = 0.0001
+)
+
+})
+
diff --git a/tests/testthat/test-TPI.R b/tests/testthat/test-TPI.R
index 3b584c6f5a3eaf5d038f9f20a04c00fcdfe2e39f..b663c40bd1e36b4792da5adb4fa95e7ea1aae7a8 100644
--- a/tests/testthat/test-TPI.R
+++ b/tests/testthat/test-TPI.R
@@ -134,16 +134,6 @@ test_that("1. Input checks", {
       month_dim = 'mth'),
   "Parameter 'month_dim' is not found in 'data' dimension."
   )
-  expect_error(
-  TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist',
-      member_dim = 1),
-  "Parameter 'member_dim' must be a character string."
-  )
-  expect_error(
-  TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist',
-      member_dim = 'ens'),
-  "Parameter 'member_dim' is not found in 'data' dimension."
-  )
 
 })
 
@@ -184,7 +174,7 @@ test_that("3. Output checks: dat2", {
  )
  expect_equal(
  mean(res2_1)*10^17,
- -2.220099,
+ -2.664535,
  tolerance = 0.00001
  ) 
 
@@ -192,7 +182,7 @@ test_that("3. Output checks: dat2", {
                indices_for_clim = 2:3)
  expect_equal(
  mean(res2_2)*10^16,
- 3.840678,
+ 4.618528,
  tolerance = 0.00001
  )
 
@@ -200,7 +190,7 @@ test_that("3. Output checks: dat2", {
                mask = mask)
  expect_equal(
  mean(res2_3)*10^17,
- -2.288135,
+ 0.8881784,
  tolerance = 0.00001
  )
 
@@ -238,7 +228,7 @@ test_that("4. Output checks: dat3", {
                mask = mask)
  expect_equal(
  mean(res3_4)*10^17,
- -1.360023,
+ 0,
  tolerance = 0.00001
  )