diff --git a/.Rbuildignore b/.Rbuildignore
index aa8227b142f77c731283cc23cb40abe0e7fdc490..19a1af65ac5ae3e96e68397e9d01cc76b0c9e8ac 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -9,3 +9,4 @@ README\.md$
 \..*\.RData$
 vignettes
 .gitlab-ci.yml
+^tests$
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index b1cf9976e0a017b4076e024242f9c1a9a7e61487..39a31ecb9326dc181d17a8f0f7ad1da5a68bce8f 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -4,7 +4,7 @@ build:
   stage: build
   script:
     - module load R
-#    - module load CDO
+    - module load CDO
     - R CMD build --resave-data .
     - R CMD check --as-cran --no-manual --run-donttest s2dverification_*.tar.gz
     - R -e 'covr::package_coverage()'
diff --git a/DESCRIPTION b/DESCRIPTION
index 0045265cfeeefd668ffd0fc634ae9e42cac29915..99565b9b212ab5ba7e422819baab425f8c4da651 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: s2dverification
 Title: Set of Common Tools for Forecast Verification
-Version: 2.8.6
+Version: 2.9.0
 Authors@R: c(
     person("BSC-CNS", role = c("aut", "cph")),
     person("Virginie", "Guemas", , "virginie.guemas@bsc.es", role = "aut"),
@@ -35,7 +35,8 @@ Description: Set of tools to verify forecasts through the computation of typical
     equations from a model, not a pure observational dataset). Intended for seasonal
     to decadal climate forecasts although can be useful to verify other kinds of
     forecasts. The package can be helpful in climate sciences for other purposes
-    than forecasting.
+    than forecasting. To find more details, see the review paper Manubens, N.et al. (2018)
+    <doi:10.1016/j.envsoft.2018.01.018>.
 Depends:
     maps,
     methods,
@@ -55,8 +56,8 @@ Suggests:
     easyVerification,
     testthat
 License: LGPL-3
-URL: https://earth.bsc.es/gitlab/es/s2dverification/wikis/home
-BugReports: https://earth.bsc.es/gitlab/es/s2dverification/issues
+URL: https://earth.bsc.es/gitlab/es/s2dverification/-/wikis/home
+BugReports: https://earth.bsc.es/gitlab/es/s2dverification/-/issues
 LazyData: true
 SystemRequirements: cdo
 Encoding: UTF-8
diff --git a/NAMESPACE b/NAMESPACE
index 4e8557898444268933dfffceb4f5d7d3c8460183..8561a4f02ff2318d1146b35d48cd4102ada054ce 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -93,6 +93,7 @@ import(methods)
 import(ncdf4)
 import(parallel)
 import(plyr)
+importFrom(abind,abind)
 importFrom(grDevices,bmp)
 importFrom(grDevices,col2rgb)
 importFrom(grDevices,colorRampPalette)
diff --git a/NEWS.md b/NEWS.md
index 6d81252db9d51a1ac86b6fa5bd5fca14ead4dc89..5e12423059279e6f35db1e21f75cdca9a692ca7e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,14 @@
-# s2dverification 2.8.6 (Release date: )
+# s2dverification 2.9.0 (Release date: 2020-10-30)
+- Correct the time retrieval in Load() when start date and the first lead time in netCDF file do not match. 
+In addition, when the lead time in each data set is not consistent, the retrieved data should change according to
+ the first data set in the same Load call.
+- One new parameter 'K' is added in Composite() to indicate the maximum number of composites.
+- Revise the per-paired method in Clim() when NA exists.
+- Correct the Corr() p-value.
+- Bugfix for CDO version reading. The error occurred when the CDO version number is followed by letters.
+- Bugfix for Ano() when obs and exp have inconsistent dimensions.
+
+# s2dverification 2.8.6 (Release date: 2019-10-17)
 - Apply Roxygen2 format to all the files.
 - Bug fix in Composite().
 - Bug fix in Ano(). Recommend to assign the dimensions by name to avoid confusion when the dimensions have same length.
diff --git a/R/Ano.R b/R/Ano.R
index 1c8ffee9059e29aac9b9c20496424b234d917263..ded1e9ac5249e5ec51ef873be837d67211a21b90 100644
--- a/R/Ano.R
+++ b/R/Ano.R
@@ -48,45 +48,113 @@ Ano <- function(var, clim) {
     clim <- InsertDim(clim, 2, dimvar[2])
   }
 
-  if ((length(dimvar) > length(dim(clim))) & (dim(clim)[2] != dimvar[2])) {
-    clim <- InsertDim(clim, 2, dimvar[2])
-  }
-  if ((length(dimvar) > length(dim(clim))) & (dim(clim)[2] == dimvar[2])) {
-    if (is.null(names(dim(clim))) | is.null(names(dimvar))) {
-      stop('Provide dimension names on parameter \'var\' and \'clim\' to avoid ambiguity.')
-    } else {
-      if (names(dim(clim)[2]) != names(dimvar[2])) {
-        clim <- InsertDim(clim, 2, dimvar[2])
+  # Check 2nd dim
+  if (length(dimvar) > length(dim(clim))) {
+
+    if (dim(clim)[2] != dimvar[2]) {
+      clim <- InsertDim(clim, 2, dimvar[2])
+    } else {  # clim[2] == dimvar[2]
+      if (is.null(names(dim(clim))) | is.null(names(dimvar))) {
+        stop('Provide dimension names on parameter \'var\' and \'clim\' to avoid ambiguity.')
+      } else {
+        if (names(dim(clim)[2]) != names(dimvar[2])) {
+          clim <- InsertDim(clim, 2, dimvar[2])
+        }
       }
     }
+
   }
 
-  if ((length(dimvar) > length(dim(clim))) & (dim(clim)[3] != dimvar[3])) {
+  # Check if clim has 3rd dim. If yes, directly go to the next if. If not, add 3rd dim.
+  if (length(dimvar) > length(dim(clim)) & length(dim(clim)) == 2) {
     clim <- InsertDim(clim, 3, dimvar[3])
   }
-  if ((length(dimvar) > length(dim(clim))) & (dim(clim)[3] == dimvar[3])) {
-    if (is.null(names(dim(clim))) | is.null(names(dimvar))) {
-      stop('Provide dimension names on parameter \'var\' and \'clim\' to avoid ambiguity.')
+
+  # Check 3rd dim
+  if (length(dimvar) > length(dim(clim))) {
+
+    if (dim(clim)[3] != dimvar[3]) {
+      clim <- InsertDim(clim, 3, dimvar[3])
+    } else {  # clim[3] == dimvar[3]
+      if (is.null(names(dim(clim))) | is.null(names(dimvar))) {
+        stop('Provide dimension names on parameter \'var\' and \'clim\' to avoid ambiguity.')
+      } else {
+        if (names(dim(clim)[3]) != names(dimvar[3])) {
+          clim <- InsertDim(clim, 3, dimvar[3])
+        } 
+      }
+    }
+
+  }
+  # Check if clim has 4th dim. If yes, directly go to the next if. If not, add 4th dim.
+  if (length(dimvar) > length(dim(clim)) & length(dim(clim)) == 3) {
+    clim <- InsertDim(clim, 4, dimvar[4])
+  }
+
+  # Check 4th dim
+  if (length(dimvar) > length(dim(clim))) {
+
+    if (dim(clim)[4] != dimvar[4]) {
+      clim <- InsertDim(clim, 4, dimvar[4])
     } else {
-      if (names(dim(clim)[3]) != names(dimvar[3])) {
-        clim <- InsertDim(clim, 3, dimvar[3])
+      if (is.null(names(dim(clim))) | is.null(names(dimvar))) {
+        stop('Provide dimension names on parameter \'var\' and \'clim\' to avoid ambiguity.')
+      } else {
+        if (names(dim(clim)[4]) != names(dimvar[4])) {
+          clim <- InsertDim(clim, 4, dimvar[4])
+        }
       }
     }
+
+  }
+  # Check if clim has 5th dim. If yes, directly go to the next if. If not, add 5th dim.
+  if (length(dimvar) > length(dim(clim)) & length(dim(clim)) == 4) {
+    clim <- InsertDim(clim, 5, dimvar[5])
   }
 
-  if ((length(dimvar) > length(dim(clim))) & (dim(clim)[4] != dimvar[4])) {
-    clim <- InsertDim(clim, 4, dimvar[4])
+  # Check 5th dim
+  if (length(dimvar) > length(dim(clim))) {
+
+    if (dim(clim)[5] != dimvar[5]) {
+      clim <- InsertDim(clim, 5, dimvar[5])
+    } else {
+      if (is.null(names(dim(clim))) | is.null(names(dimvar))) {
+        stop('Provide dimension names on parameter \'var\' and \'clim\' to avoid ambiguity.')
+      } else {
+        if (names(dim(clim)[5]) != names(dimvar[5])) {
+          clim <- InsertDim(clim, 5, dimvar[5])
+        }
+      }
+    }
+
+  }
+  # Check if clim has 6th dim. If yes, directly go to the next if. If not, add 6th dim.
+  if (length(dimvar) > length(dim(clim)) & length(dim(clim)) == 5) {
+    clim <- InsertDim(clim, 6, dimvar[6])
   }
-  if ((length(dimvar) > length(dim(clim))) & (dim(clim)[4] == dimvar[4])) {
-    if (is.null(names(dim(clim))) | is.null(names(dimvar))) {
-      stop('Provide dimension names on parameter \'var\' and \'clim\' to avoid ambiguity.')
+
+  # Check 6th dim
+  if (length(dimvar) > length(dim(clim))) {
+
+    if (dim(clim)[6] != dimvar[6]) {
+      clim <- InsertDim(clim, 6, dimvar[6])
     } else {
-      if (names(dim(clim)[4]) != names(dimvar[4])) {
-        clim <- InsertDim(clim, 4, dimvar[4])
+      if (is.null(names(dim(clim))) | is.null(names(dimvar))) {
+        stop('Provide dimension names on parameter \'var\' and \'clim\' to avoid ambiguity.')
+      } else {
+        if (names(dim(clim)[6]) != names(dimvar[6])) {
+          clim <- InsertDim(clim, 6, dimvar[6])
+        }
       }
     }
+
+  }
+  # Check if clim has 7th dim. If yes, directly go to the next if. If not, add 7th dim.
+  if (length(dimvar) > length(dim(clim)) & length(dim(clim)) == 6) {
+    clim <- InsertDim(clim, 7, dimvar[7])
   }
 
+
   #
   #  Raw anomalies
   # ~~~~~~~~~~~~~~~
diff --git a/R/CDORemap.R b/R/CDORemap.R
index a04d88a536447a2488283d5e5c3991d079490012..ea6ff1374333384d7b7298b4d5bb19f1e529b181 100644
--- a/R/CDORemap.R
+++ b/R/CDORemap.R
@@ -406,6 +406,7 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
   }
   if (is.logical(crop)) {
     if (crop) {
+            warning("Parameter 'crop' = 'TRUE'. The output grid range will follow the input lons and lats.")
       if (length(lons) == 1 || length(lats) == 1) {
         stop("CDORemap cannot remap if crop = TRUE and values for only one ",
              "longitude or one latitude are provided. Either a) provide ",
@@ -516,6 +517,8 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
         }
 ###---
       }
+    }  else if (crop == FALSE) {
+       warning("Parameter 'crop' = 'FALSE'. The output grid range will follow parameter 'grid'.")
     }
   } else if (is.numeric(crop)) {
     if (length(crop) != 4) {
diff --git a/R/Clim.R b/R/Clim.R
index 8445e46277bcba517dd1256d957dbe3347c1b0fb..a60106a3ffc3ca4e70ce850eff4b93d90665caab 100644
--- a/R/Clim.R
+++ b/R/Clim.R
@@ -70,20 +70,50 @@ Clim <- function(var_exp, var_obs, memb = TRUE, kharin = FALSE, NDV = FALSE) {
   #  Find common points to compute climatologies 
   # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   #
-  tmp <- MeanListDim(var_obs, dims = 5:7, narm = TRUE)
-  tmp2 <- MeanListDim(var_exp, dims = 5:7, narm = TRUE)
-  nan <- MeanListDim(tmp, dims = 1:2, narm = FALSE) + Mean1Dim(Mean1Dim(tmp2, 2, 
-                                                                    narm = TRUE), 
-                                                           1, narm = FALSE)
-  for (jdate in 1:dimexp[3]) {
-    for (jtime in 1:dimexp[4]) {
-      if (is.na(nan[jdate, jtime])) {
-        var_exp[, , jdate, jtime, , , ] <- NA
-        var_obs[, , jdate, jtime, , , ] <- NA
-      }
+
+  na_array <- array(0, dim = dim(var_exp)[3:7])
+
+  for (i_dat in 1:dimexp[1]) {
+    for (i_memb in 1:dimexp[2]) {
+      na_array <- na_array + as.numeric(is.na(var_exp[i_dat, i_memb, , , , , ]))
+    }
+  }
+  for (i_dat in 1:dimobs[1]) {
+    for (i_memb in 1:dimobs[2]) {
+      na_array <- na_array + as.numeric(is.na(var_obs[i_dat, i_memb, , , , , ]))
+    }
+  }
+  na_array <- as.logical(na_array)  # TRUE is NA, FALSE is not
+
+  for (i_dat in 1:dimexp[1]) {
+    for (i_memb in 1:dimexp[2]) {      
+      asd <- var_exp[i_dat, i_memb, , , , , ]
+      asd[which(na_array)] <- NA
+      var_exp[i_dat, i_memb, , , , , ] <- asd
+    }
+  }
+  for (i_dat in 1:dimobs[1]) {
+    for (i_memb in 1:dimobs[2]) {
+      asd <- var_obs[i_dat, i_memb, , , , , ]
+      asd[which(na_array)] <- NA
+      var_obs[i_dat, i_memb, , , , , ] <- asd
     }
   }
 
+#  nan <- MeanListDim(var_exp, dims = c(1:2, 5:7), narm = FALSE) + 
+#         MeanListDim(var_obs, dims = c(1:2, 5:7), narm = FALSE)
+#
+#  for (jdate in 1:dimexp[3]) {
+#    for (jtime in 1:dimexp[4]) {
+#      if (is.na(nan[jdate, jtime])) {
+#        var_exp[, , jdate, jtime, , , ] <- NA
+#        var_obs[, , jdate, jtime, , , ] <- NA
+#      }
+#    }
+#  }
+
+
+
   #
   #  Compute climatologies 
   # ~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/R/Composite.R b/R/Composite.R
index 865501d6c91a3e5d672093471cbdbbdff606e080..f68604971373b1fc9cdadd1af01e5640557ff86c 100644
--- a/R/Composite.R
+++ b/R/Composite.R
@@ -17,6 +17,7 @@
 #'  use +2 occurrences (i.e., shifted 2 time steps forward). Default is lag = 0.
 #'@param eno For using the effective sample size (TRUE) or the total sample 
 #'  size (FALSE that is the default) for the number of degrees of freedom.
+#'@param K numeric value indicating the maximum number of composites. By default (NULL), it takes the maximum value provided in occ.
 #'@param fileout Name of the .sav output file (NULL is the default).
 #'
 #'@return 
@@ -67,18 +68,27 @@
 #'occ2[c(3, 9, 15, 21)] <- 1
 #'
 #'filled.contour(Composite(var=f1, occ=occ2)$composite[,,1])
+#'
+#'# Example with one missing composite (#3) in occ:
+#'data <- 1 : (4 * 5 * 6)
+#'dim(data) <- c(lon = 4, lat = 5, case = 6)
+#'occ <- c(1, 1, 2, 2, 3, 3) 
+#'res <- Composite(data, occ, K = 4)
+#'
 #'@importFrom stats sd pt
 #'@export
-Composite <- function(var, occ, lag = 0, eno = FALSE, fileout = NULL) {
+Composite <- function(var, occ, lag = 0, eno = FALSE, K = NULL, fileout = NULL) {
 
   if ( dim(var)[3] != length(occ) ) {
      stop("Temporal dimension of var is not equal to length of occ.")
   }
-  K         <- max(occ)
-  composite <- array(dim = c(dim(var)[1:2], K))
+  if (is.null(K)) {
+     K <- max(occ)
+  }
+  composite <- array(dim = c(dim(var)[1:2], composite = K))
   tvalue    <- array(dim = dim(var)[1:2])
   dof       <- array(dim = dim(var)[1:2])
-  pvalue    <- array(dim = c(dim(var)[1:2], K))
+  pvalue    <- array(dim = c(dim(var)[1:2], composite = K))
 
   if (eno == TRUE) { 
     n_tot <- Eno(var, posdim = 3)
@@ -89,34 +99,34 @@ Composite <- function(var, occ, lag = 0, eno = FALSE, fileout = NULL) {
   stdv_tot <- apply(var, c(1, 2), sd, na.rm = TRUE) 
 
   for (k in 1 : K) {
+    if (length(which(occ == k)) >= 1) {
+      indices <- which(occ == k) + lag
+      toberemoved <-  which(0 > indices | indices > dim(var)[3])
 
-    indices <- which(occ == k) + lag
-    toberemoved <-  which(0 > indices | indices > dim(var)[3])
-
-    if (length(toberemoved) > 0) {
+      if (length(toberemoved) > 0) {
         indices <- indices[-toberemoved]
-    }
-    if (eno == TRUE) {
+      }
+      if (eno == TRUE) {
         n_k <- Eno(var[, , indices], posdim = 3)
-    } else {
+      } else {
         n_k <- length(indices)
-    }
-    if (length(indices) == 1) {
+      }
+      if (length(indices) == 1) {
         composite[, , k] <- var[, , indices] 
         warning(paste("Composite", k, "has length 1 and pvalue is NA."))
-    }  else {
+      }  else {
         composite[, , k] <- Mean1Dim(var[, , indices], posdim = 3, narm = TRUE)
-    }
-    stdv_k <- apply(var[, , indices], c(1, 2), sd, na.rm = TRUE)
+      }
+      stdv_k <- apply(var[, , indices], c(1, 2), sd, na.rm = TRUE)
     
-    tvalue <- (mean_tot - composite[, , k]) / 
-               sqrt(stdv_tot ^ 2 / n_tot + stdv_k ^ 2 / n_k)
-    dof <- (stdv_tot ^ 2 / n_tot + stdv_k ^ 2 / n_k) ^ 2 / 
-           ((stdv_tot ^ 2 / n_tot) ^ 2 / (n_tot - 1) +
-           (stdv_k ^ 2 / n_k) ^ 2 / (n_k - 1))
-    pvalue[, , k] <- 2 * pt(-abs(tvalue), df = dof)
-  }  
-
+      tvalue <- (mean_tot - composite[, , k]) / 
+                 sqrt(stdv_tot ^ 2 / n_tot + stdv_k ^ 2 / n_k)
+      dof <- (stdv_tot ^ 2 / n_tot + stdv_k ^ 2 / n_k) ^ 2 / 
+             ((stdv_tot ^ 2 / n_tot) ^ 2 / (n_tot - 1) +
+             (stdv_k ^ 2 / n_k) ^ 2 / (n_k - 1))
+      pvalue[, , k] <- 2 * pt(-abs(tvalue), df = dof)
+    }  
+  }
   if (is.null(fileout) == FALSE) { 
     output <- list(composite = composite, pvalue = pvalue)   
     save(output, file = paste(fileout, '.sav', sep = ''))
diff --git a/R/ConfigFileOpen.R b/R/ConfigFileOpen.R
index 182b205e31f898fbec7ab32b4a6e9c699de14849..d0a4ada0d0477217529591653ca67fd3d2d74840 100644
--- a/R/ConfigFileOpen.R
+++ b/R/ConfigFileOpen.R
@@ -147,7 +147,7 @@
 #'  ConfigEditEntry, ConfigFileOpen, ConfigShowSimilarEntries, ConfigShowTable
 #'@references
 #'[1] \url{https://stat.ethz.ch/R-manual/R-devel/library/base/html/regex.html}\cr
-#'[2] \url{http://tldp.org/LDP/abs/html/globbingref.html}
+#'[2] \url{https://tldp.org/LDP/abs/html/globbingref.html}
 #'@author History:
 #'  0.1 - 2015-05 (N. Manubens, \email{nicolau.manubens@@ic3.cat}) - First version
 #'  1.0 - 2015-11 (N. Manubens, \email{nicolau.manubens@@ic3.cat}) - Removed grid column and storage formats
diff --git a/R/Corr.R b/R/Corr.R
index 419002305586a9920628244c6445ee7f7ef20774..fe8e5716814fcaeb61252163e0237a51ee6ceaf0 100644
--- a/R/Corr.R
+++ b/R/Corr.R
@@ -198,9 +198,13 @@ Corr <- function(var_exp, var_obs, posloop = 1, poscor = 2, compROW = NULL,
                         }
                         if (pval) {
                           #t <- qt(0.95, eno - 2)
-                          t <- qt(siglev, eno - 2)
-                          CORR[jexp, jobs, dim_pval, j3, j4, j5, j6, j7, j8, j9, 
-                               j10] <- sqrt((t * t) / ((t * t) + eno - 2))
+#                          t <- qt(siglev, eno - 2)
+#                          CORR[jexp, jobs, dim_pval, j3, j4, j5, j6, j7, j8, j9, 
+#                               j10] <- sqrt((t * t) / ((t * t) + eno - 2))
+                          t <-sqrt(toto * toto * (eno - 2) / (1 - (toto ^ 2)))
+                          CORR[jexp, jobs, dim_pval, j3, j4, j5, j6, j7, j8, j9,
+                               j10] <- pt(t, eno - 2, lower.tail = FALSE)
+
                         }
                         if (conf) {
                           CORR[jexp, jobs, 1, j3, j4, j5, j6, j7, j8, j9, 
diff --git a/R/Load.R b/R/Load.R
index e17f3cfa683dbaf839d5f3cb022e091bec0352ff..9feec5d37bb6aa9cdd00cd2687675c7fca10d63a 100644
--- a/R/Load.R
+++ b/R/Load.R
@@ -840,6 +840,7 @@
 #'  } 
 #'@import parallel bigmemory methods
 #'@importFrom stats ts window na.omit
+#'@importFrom abind abind
 #'@export
 Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, 
                  nmemberobs = NULL, nleadtime = NULL, leadtimemin = 1, 
@@ -850,6 +851,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
                  configfile = NULL, varmin = NULL, varmax = NULL, 
                  silent = FALSE, nprocs = NULL, dimnames = NULL, 
                  remapcells = 2, path_glob_permissive = 'partial') {
+
   #library(parallel)
   #library(bigmemory)
 
@@ -1052,7 +1054,6 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
   if (!is.character(storefreq) || !(storefreq %in% c('monthly', 'daily'))) {
     stop("Error: parameter 'storefreq' is wrong, can take value 'daily' or 'monthly'.")
   }
-
   # sampleperiod
   if (is.null(sampleperiod) || !is.numeric(sampleperiod)) {
     stop("Error: parameter 'sampleperiod' is wrong. It should be numeric.")
@@ -1382,6 +1383,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
   dims2define <- TRUE
   is_file_per_member_exp <- rep(nmod, FALSE)
   exp_work_pieces <- list()
+  first_time_step_list <- NULL
   jmod <- 1
   while (jmod <= nmod) {
     first_dataset_file_found <- FALSE
@@ -1427,6 +1429,21 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
       if (is_file_per_member_exp[jmod]) {
         replace_values[["MEMBER_NUMBER"]] <- '*'
       }
+      if (jsdate == 1) {
+        work_piecetime <- list(dataset_type = dataset_type,
+                           filename = .ConfigReplaceVariablesInString(quasi_final_path,
+                                                                      replace_values),
+                           namevar = namevar, grid = grid, remap = remap, 
+                           remapcells = remapcells,
+                           is_file_per_member = is_file_per_member_exp[jmod],
+                           is_file_per_dataset = FALSE,
+                           lon_limits = c(lonmin, lonmax),
+                           lat_limits = c(latmin, latmax), dimnames = exp[[jmod]][['dimnames']],
+                           single_dataset = single_dataset)
+        looking_time <- .LoadDataFile(work_piecetime, explore_dims = TRUE, 
+                                      silent = silent)
+        first_time_step_list <- c(first_time_step_list, list(looking_time$time_dim))
+      }
       # If the dimensions of the output matrices are still to define, we try to read
       # the metadata of the data file that corresponds to the current iteration
       if (dims2define) {
@@ -1531,6 +1548,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
       jsdate <- jsdate + 1
     }
     replace_values[extra_vars] <- NULL
+    #first_dataset_file_found <- FALSE
     jmod <- jmod + 1
   }
   if (dims2define && length(exp) > 0) {
@@ -1553,7 +1571,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
   # the current date.
   if (is.null(exp) || dims == 0) {
     if (is.null(leadtimemax)) {
-      diff <- Sys.time() - as.POSIXct(sdates[1], format = '%Y%m%d')
+      diff <- Sys.time() - as.POSIXct(sdates[1], format = '%Y%m%d', tz = "UTC")
       if (diff > 0) {
         .warning("Loading observations only and no 'leadtimemax' specified. Data will be loaded from each starting date to current time.")
         if (storefreq == 'monthly') { 
@@ -1571,7 +1589,52 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
     }
     leadtimes <- seq(leadtimemin, leadtimemax, sampleperiod)
   }
-  
+  # If there are differences in the first time stamp in exp files:
+  if (!is.null(exp)) {
+    in_date <- lapply(first_time_step_list, function(x) {
+        origin <- as.POSIXct(
+          paste(strsplit(x$time_units, " ")[[1]][c(3,4)],
+                collapse = " "), tz = 'UTC')
+        units <- strsplit(x$time_units, " ")[[1]][1]
+        if (units == 'hours') {
+          exp_first_time_step <- as.POSIXct(
+            x$first_time_step_in_file *
+            3600, origin = origin, tz = 'UTC')
+        } else if (units == 'days') {
+          exp_first_time_step <- as.POSIXct(
+            x$first_time_step_in_file *
+            86400, origin = origin, tz = 'UTC')
+        }
+        day <- as.numeric(format(exp_first_time_step, "%d"))
+        return(day)
+    })
+    exp_first_time_step <- min(unlist(in_date))
+    if (max(unlist(in_date)) > 1) {
+      leadtimes <- seq(exp_first_time_step, leadtimemax + max(unlist(in_date)) - 1, 
+                       sampleperiod)
+    }
+    if (leadtimemin > 1 & length(in_date) > 1) {
+      lags <- lapply(in_date, function(x) {x - in_date[[1]]})
+      new_leadtimemin <- lapply(lags, function(x) {leadtimemin - x})
+      new_leadtimemax <- lapply(lags, function(x) {leadtimemax - x})
+      jmod <- 2
+      npieces <- length(exp_work_pieces)/nmod
+      while (jmod <= nmod) {
+        jpiece <- 1
+        while (jpiece <= npieces) {
+          exp_work_pieces[[npieces * (jmod - 1) + jpiece]]$leadtimes <- 
+            seq(new_leadtimemin[[jmod]], new_leadtimemax[[jmod]], sampleperiod)
+          jpiece <- jpiece + 1
+        }
+        jmod <- jmod + 1
+      }
+    }
+    lag <- 1 - in_date[[1]]
+    leadtimes <- seq(leadtimemin - lag, leadtimemax #+ max(unlist(in_date)) + lag,
+                     - lag,
+                     sampleperiod)
+    exp_first_time_step <- leadtimemin - lag
+  }
   # Now we start iterating over observations. We try to find the output matrix
   # dimensions and we build anyway the work pieces corresponding to the observational
   # data that time-corresponds the experimental data or the time-steps until the
@@ -1635,6 +1698,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
         found_data <- .LoadDataFile(work_piece, explore_dims = TRUE, silent = silent)
         found_dims <- found_data$dims
         var_long_name <- found_data$var_long_name
+        first_time_step_list <- c(first_time_step_list, list(found_data$time_dim))
         units <- found_data$units
         if (!is.null(found_dims)) {
           is_2d_var <- found_data$is_2d_var
@@ -1709,7 +1773,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
           }
           day <- as.integer(day)
           startdate <- as.POSIXct(paste(substr(sdate, 1, 4), '-',
-                       substr(sdate, 5, 6), '-', day, ' 12:00:00', sep = '')) + 
+                       substr(sdate, 5, 6), '-', day, ' 12:00:00', sep = ''), tz = "UTC") + 
                        (leadtimemin - 1) * 86400
           year <- as.integer(substr(startdate, 1, 4))
           month <- as.integer(substr(startdate, 6, 7))
@@ -1732,7 +1796,18 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
             ## This condition must be fulfilled to put all the month time steps
             ## in the dimension of length nleadtimes. Otherwise it must be cut:
             #(length(leadtimes) - 1)*sampleperiod + 1 - (jleadtime - 1)*sampleperiod >= days_in_month - day + 1
-            obs_file_indices <- seq(day, min(days_in_month, (length(leadtimes) - jleadtime) * sampleperiod + day), sampleperiod)
+            
+            ## The first time step in exp could be different from sdate:
+            if (jleadtime == 1 & !is.null(exp)) {
+             if (is.null(first_time_step_list[[1]])) {
+                stop("Check 'time' variable in the experimental files ",
+                     "since not units or first time step have been found.")
+              } else {
+                  day <- leadtimes[1]
+              }
+            }
+          obs_file_indices <- seq(day, min(days_in_month,
+                  (length(leadtimes) - jleadtime) * sampleperiod + day), sampleperiod)
           } else {
             obs_file_indices <- 1
           }
@@ -1828,7 +1903,8 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
           }
           
           if (storefreq == 'daily') {
-            startdate <- startdate + 86400 * sampleperiod * length(obs_file_indices)
+            startdate <- startdate + 86400 * sampleperiod * 
+                         max(obs_file_indices)
             year <- as.integer(substr(startdate, 1, 4))
             month <- as.integer(substr(startdate, 6, 7))
             day <- as.integer(substr(startdate, 9, 10))
@@ -2142,49 +2218,54 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
     attr(variable, 'daily_agg_cellfun') <- 'none'
     attr(variable, 'monthly_agg_cellfun') <- 'none'
     attr(variable, 'verification_time') <- 'none'
-
-    number_ftime <- NULL 
+    
+    number_ftime <- NULL
     if (is.null(var_exp)) {
       mod_data <- NULL
-    } else { 
+    } else {
       dim_reorder <- length(dim_exp):1
       dim_reorder[2:3] <- dim_reorder[3:2]
       old_dims <- dim_exp
       dim_exp <- dim_exp[dim_reorder]
-      mod_data <- aperm(array(bigmemory::as.matrix(var_exp), dim = old_dims), dim_reorder)
+      mod_data <-
+        aperm(array(bigmemory::as.matrix(var_exp), dim = old_dims), dim_reorder)
       attr(mod_data, 'dimensions') <- names(dim_exp)
       names(dim(mod_data)) <- names(dim_exp)
       number_ftime <- dim_exp[["ftime"]]
     }
-
+    
     if (is.null(var_obs)) {
       obs_data <- NULL
-    } else {                     
+    } else {
       dim_reorder <- length(dim_obs):1
       dim_reorder[2:3] <- dim_reorder[3:2]
       old_dims <- dim_obs
       dim_obs <- dim_obs[dim_reorder]
-      obs_data <- aperm(array(bigmemory::as.matrix(var_obs), dim = old_dims), dim_reorder)
+      obs_data <-
+        aperm(array(bigmemory::as.matrix(var_obs), dim = old_dims), dim_reorder)
       attr(obs_data, 'dimensions') <- names(dim_obs)
       names(dim(obs_data)) <- names(dim_obs)
       if (is.null(number_ftime)) {
         number_ftime <- dim_obs[["ftime"]]
       }
     }
-
     if (is.null(latitudes)) {
       lat <- 0
       attr(lat, 'cdo_grid_name') <- 'none'
       attr(lat, 'first_lat') <- 'none'
-      attr(lat, 'last_lat') <- 'none' 
+      attr(lat, 'last_lat') <- 'none'
     } else {
       lat <- latitudes
-      attr(lat, 'cdo_grid_name') <- if (is.null(grid)) 'none' else grid
+      attr(lat, 'cdo_grid_name') <-
+        if (is.null(grid))
+          'none'
+      else
+        grid
       attr(lat, 'first_lat') <- tail(lat, 1)
       attr(lat, 'last_lat') <- head(lat, 1)
     }
     attr(lat, 'projection') <- 'none'
-
+    
     if (is.null(longitudes)) {
       lon <- 0
       attr(lon, 'cdo_grid_name') <- 'none'
@@ -2194,14 +2275,18 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
       attr(lon, 'last_lon') <- 'none'
     } else {
       lon <- longitudes
-      attr(lon, 'cdo_grid_name') <- if (is.null(grid)) 'none' else grid
+      attr(lon, 'cdo_grid_name') <-
+        if (is.null(grid))
+          'none'
+      else
+        grid
       attr(lon, 'data_across_gw') <- data_across_gw
       attr(lon, 'array_across_gw') <- array_across_gw
       attr(lon, 'first_lon') <- lon[which.min(abs(lon - lonmin))]
       attr(lon, 'last_lon') <- lon[which.min(abs(lon - lonmax))]
     }
     attr(lon, 'projection') <- 'none'
-
+    
     dates <- list()
     ## we must put a start and end time for each prediction c(start date, forecast time)
     if (storefreq == 'minutely') {
@@ -2209,26 +2294,46 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
     } else if (storefreq == 'hourly') {
       store_period <- 'hour'
     } else if (storefreq == 'daily') {
-      store_period <- 'day'
+      store_period <- 'DSTday'
     } else if (storefreq == 'monthly') {
       store_period <- 'month'
     }
-
+    
     addTime <- function(date, period, n = 1) {
       seq(date, by = paste(n, period), length = 2)[2]
     }
-
+    
     # We build dates, a list with components start and end.
     # Start is a list with as many components as start dates.
     # Each component is a vector of the initial POSIXct date of each
     # forecast time step
-    dates[["start"]] <- do.call(c, lapply(sdates,
+    if (!is.null(exp)) {
+      if (storefreq == 'daily' & leadtimes[[1]] > 1) {
+        origin <- leadtimes[[1]] - 1
+        leadtimemin <- 1
+      } else {
+        origin <- 0
+      }
+      dates[["start"]] <- do.call(c, lapply(sdates,
+          function(x) {
+            do.call(c, lapply((origin:(origin + number_ftime - 1)) * sampleperiod,
+              function(y) {
+                 addTime(as.POSIXct(x, format = "%Y%m%d", tz = "UTC"),
+                                    store_period, y + leadtimemin - 1)
+              }))
+           }))
+    } else {
+      origin <- 0
+      dates[["start"]] <- do.call(c, lapply(sdates,
       function(x) {
         do.call(c, lapply((0:(number_ftime - 1)) * sampleperiod,
           function(y) {
-            addTime(as.POSIXct(x, format = "%Y%m%d"), store_period, y + leadtimemin - 1)
+            addTime(as.POSIXct(x, format = "%Y%m%d", tz = "UTC"),
+                    store_period, y + leadtimemin - 1)
           }))
       }))
+    }
+    attr(dates[["start"]], "tzone") <- "UTC"
     # end is similar to start, but contains the end dates of each forecast 
     # time step
     dates[["end"]] <- do.call(c, lapply(dates[["start"]],
@@ -2238,6 +2343,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
             addTime(y, store_period)
           }))
       }))
+    attr(dates[["end"]], "tzone") <- "UTC"
 
     tags_to_find <- c('START_DATE', 'MEMBER_NUMBER', 'YEAR', 'MONTH', 'DAY')
     position_of_tags <- na.omit(match(tags_to_find, names(replace_values)))
@@ -2249,14 +2355,19 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
       models <- list()
       for (jmod in 1:length(exp)) {
         member_names <- paste0("Member_", 1:nmember[jmod])
-        models[[exp[[jmod]][["name"]]]] <- list(
-          InitializationDates = lapply(member_names,
+
+        InitializationDates <- do.call(c, lapply(member_names,
             function(x) {
               do.call(c, lapply(sdates, function(y) {
-                as.POSIXct(y, format = "%Y%m%d")
+                as.POSIXct(y, format = "%Y%m%d", tz = "UTC")
               }))
-            }),
+            }))
+        attr(InitializationDates, "tzone") <- "UTC"
+
+        models[[exp[[jmod]][["name"]]]] <- list(
+          InitializationDates = InitializationDates,
           Members = member_names)
+
         names(models[[exp[[jmod]][["name"]]]]$InitializationDates) <- member_names
         attr(models[[exp[[jmod]][["name"]]]], 'dataset') <- exp[[jmod]][["name"]]
         attr(models[[exp[[jmod]][["name"]]]], 'source') <- {
diff --git a/R/Utils.R b/R/Utils.R
index 4460018139d3bf55bfd8cc3f9b1b2a4f6f17f850..1dbf48727cb3e9f5b1e05e8fc57d53231bd62e5d 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -206,7 +206,11 @@
       if (Sys.which("cdo")[[1]] == "") {
         stop("Error: CDO libraries not available")
       }
-      cdo_version <- as.numeric_version(strsplit(suppressWarnings(system2("cdo", args = '-V', stderr = TRUE))[[1]], ' ')[[1]][5])
+
+     cdo_version <- strsplit(suppressWarnings(system2("cdo", args = '-V', stderr = TRUE))[[1]], ' ')[[1]][5]
+
+      cdo_version <- as.numeric_version(unlist(strsplit(cdo_version, "[A-Za-z]", fixed = FALSE))[[1]])
+
     }
     # If the variable to load is 2-d, we need to determine whether:
     #  - interpolation is needed
@@ -591,6 +595,9 @@
         nltime <- fnc$var[[namevar]][['dim']][[match(time_dimname, var_dimnames)]]$len
         expected_dims <- c(expected_dims, time_dimname)
         dim_matches <- match(expected_dims, var_dimnames)
+        first_time_step_in_file <- fnc$var[[namevar]][['dim']][[match(time_dimname,
+                                        var_dimnames)]]$vals[1]
+        time_units <- fnc$var[[namevar]][['dim']][[match(time_dimname, var_dimnames)]]$units
       } else {
         if (!is.null(old_members_dimname)) {
           expected_dims[which(expected_dims == 'lev')] <- old_members_dimname
@@ -631,7 +638,7 @@
         time_indices <- ts(time_indices, start = c(years[1], mons[1]), 
                            end = c(years[length(years)], mons[length(mons)]),
                            frequency = 12)
-        ltimes_list <- list()
+                ltimes_list <- list()
         for (sdate in work_piece[['startdates']]) {
           selected_time_indices <- window(time_indices, start = c(as.numeric(
                                    substr(sdate, 1, 4)), as.numeric(substr(sdate, 5, 6))), 
@@ -938,7 +945,9 @@
   if (explore_dims) {
     list(dims = dims, is_2d_var = is_2d_var, grid = grid_name, 
          units = units, var_long_name = var_long_name, 
-         data_across_gw = data_across_gw, array_across_gw = array_across_gw)
+         data_across_gw = data_across_gw, array_across_gw = array_across_gw,
+         time_dim = list(first_time_step_in_file = first_time_step_in_file,
+                        time_units = time_units))
   } else {
     ###if (!silent && !is.null(progress_connection) && !is.null(work_piece[['progress_amount']])) {
     ###  foobar <- writeBin(work_piece[['progress_amount']], progress_connection)
diff --git a/R/s2dverification.R b/R/s2dverification.R
index d332f38dea348be4a5071dc174f46d5a6fc43c53..f22f06eb235b15e656ee0fc9813cbe8c31bb3b65 100644
--- a/R/s2dverification.R
+++ b/R/s2dverification.R
@@ -11,12 +11,12 @@
 #'  \tabular{ll}{ 
 #'Package: \tab s2dverification\cr 
 #'Type: \tab Package\cr 
-#'Version: \tab 2.8.4\cr 
-#'Date: \tab 2019-01-30\cr 
+#'Version: \tab 2.9.0\cr 
+#'Date: \tab 2020-10-30\cr 
 #'License: \tab LGPLv3\cr 
 #'  } 
 #'Check an overview of the package functionalities and its modules at  
-#'\url{https://earth.bsc.es/gitlab/es/s2dverification/wikis/home}. 
+#'\url{https://earth.bsc.es/gitlab/es/s2dverification/-/wikis/home}. 
 #'For more information load the package and check the help for each function  
 #'or the documentation attached to the package. 
 #' 
diff --git a/inst/config/BSC.conf b/inst/config/BSC.conf
index d9d617eb1a4f505013fd9af4bf7b998e28176015..aca91dbb44fbb0fc2212fbcecfa06ef3cd2a0f0a 100644
--- a/inst/config/BSC.conf
+++ b/inst/config/BSC.conf
@@ -22,7 +22,7 @@ DEFAULT_DIM_NAME_MEMBERS = ensemble
 # Helper variables
 EXP_FILE = $VAR_NAME$_*$START_DATE$*.nc
 OBS_FILE = $VAR_NAME$_$YEAR$$MONTH$*.nc
-RECON_LIST = (20thcr_v2|copernicus012|ds083.2|ecco|era40|era40-erainterim|eraclim|erainterim|erainterim-lores|eraland|gecco_v2|gfs|glorys2_v1|glorys2_v3|glosea5|ishii-kimoto|jra55|merra|merra_v2|ncep-reanalysis|oaflux|oi_v2|orap5|piomas|seaice-lim2|sst|tropflux|nemovar_system4)
+RECON_LIST = (20thcr_v2|copernicus012|ds083.2|ecco|era40|era40-erainterim|eraclim|erainterim|erainterim-lores|eraland|gecco_v2|gfs|glorys2_v1|glorys2_v3|glosea5|ishii-kimoto|jra55|merra|merra_v2|ncep-reanalysis|oaflux|oi_v2|orap5|piomas|seaice-lim2|sst|tropflux|nemovar_system4|era5)
 
 # Defaults for extra variables from Load
 DEFAULT_FILE_GRID = 
diff --git a/man/Composite.Rd b/man/Composite.Rd
index 4e990020d18f83c62e7104c49cd39c3985ab62ef..3d9a3fd8eba3c149af352eb0027b8688dce3e8c4 100644
--- a/man/Composite.Rd
+++ b/man/Composite.Rd
@@ -4,7 +4,7 @@
 \alias{Composite}
 \title{Computes composites}
 \usage{
-Composite(var, occ, lag = 0, eno = FALSE, fileout = NULL)
+Composite(var, occ, lag = 0, eno = FALSE, K = NULL, fileout = NULL)
 }
 \arguments{
 \item{var}{3-dimensional array (x, y, time) containing the variable to 
@@ -23,6 +23,8 @@ use +2 occurrences (i.e., shifted 2 time steps forward). Default is lag = 0.}
 \item{eno}{For using the effective sample size (TRUE) or the total sample 
 size (FALSE that is the default) for the number of degrees of freedom.}
 
+\item{K}{numeric value indicating the maximum number of composites. By default (NULL), it takes the maximum value provided in occ.}
+
 \item{fileout}{Name of the .sav output file (NULL is the default).}
 }
 \value{
@@ -76,6 +78,13 @@ occ2 <- rep(0, 30)
 occ2[c(3, 9, 15, 21)] <- 1
 
 filled.contour(Composite(var=f1, occ=occ2)$composite[,,1])
+
+# Example with one missing composite (#3) in occ:
+data <- 1 : (4 * 5 * 6)
+dim(data) <- c(lon = 4, lat = 5, case = 6)
+occ <- c(1, 1, 2, 2, 3, 3) 
+res <- Composite(data, occ, K = 4)
+
 }
 \author{
 History:
diff --git a/man/ConfigFileOpen.Rd b/man/ConfigFileOpen.Rd
index cff7427e181a886f1ba2822c9aa9923059e5b26a..3630c155cc7f9db9f56bd1537df613496370318e 100644
--- a/man/ConfigFileOpen.Rd
+++ b/man/ConfigFileOpen.Rd
@@ -193,7 +193,7 @@ History:
 }
 \references{
 [1] \url{https://stat.ethz.ch/R-manual/R-devel/library/base/html/regex.html}\cr
-[2] \url{http://tldp.org/LDP/abs/html/globbingref.html}
+[2] \url{https://tldp.org/LDP/abs/html/globbingref.html}
 }
 \seealso{
 ConfigApplyMatchingEntries, ConfigEditDefinition, 
diff --git a/man/s2dverification.Rd b/man/s2dverification.Rd
index 698fa4edbc4e4d6268427931832bd6cdc7f036b3..47acf220b8bc3e319180fc219cb23a0cd24acd74 100644
--- a/man/s2dverification.Rd
+++ b/man/s2dverification.Rd
@@ -18,12 +18,12 @@ purposes than forecasting.
 \tabular{ll}{ 
 Package: \tab s2dverification\cr 
 Type: \tab Package\cr 
-Version: \tab 2.8.4\cr 
-Date: \tab 2019-01-30\cr 
+Version: \tab 2.9.0\cr 
+Date: \tab 2020-10-30\cr 
 License: \tab LGPLv3\cr 
  } 
 Check an overview of the package functionalities and its modules at  
-\url{https://earth.bsc.es/gitlab/es/s2dverification/wikis/home}. 
+\url{https://earth.bsc.es/gitlab/es/s2dverification/-/wikis/home}. 
 For more information load the package and check the help for each function  
 or the documentation attached to the package.
 }
diff --git a/s2dverification-manual.pdf b/s2dverification-manual.pdf
index 0392b213206b1f57fe404dee8b31e4b81d31694c..67373123d612604dcdd3746f955af9a30fe4f665 100644
Binary files a/s2dverification-manual.pdf and b/s2dverification-manual.pdf differ
diff --git a/tests/testthat/test-Composite.R b/tests/testthat/test-Composite.R
index f8edf1cd4596e844e78cc8edff14aa9441323e73..6c7eb1def99c80809dbefdd628f288dc42fd0690 100644
--- a/tests/testthat/test-Composite.R
+++ b/tests/testthat/test-Composite.R
@@ -11,7 +11,7 @@ test_that("Sanity checks", {
 
   var <- array(rep(c(1, 3, 2, 1, 2), 8), dim = c(x = 2, y = 4, time = 5))
   occ <- c(1, 2, 2, 2, 1)
-  output <- c(x = 2, y = 4, 2)  #dim(asd$composite)
+  output <- c(x = 2, y = 4, composite = 2)  #dim(asd$composite)
   expect_equal(
     dim(Composite(var, occ)$composite),
     output
@@ -29,6 +29,20 @@ test_that("Sanity checks", {
     Composite(var, occ)$composite[, , 2],
     output
   )
+  expect_equal(
+    dim(Composite(var, occ, K = 2)$composite),
+    c(x = 2, y = 4, composite = 2)
+  )
+  expect_equal(
+    Composite(var, occ, K = 3),
+    Composite(var, occ)
+  )
+  expect_equal(
+    Composite(var, occ, K = 4)$composite[1, 1, ],
+    c(1.5, 1.5, 2.0, NA),
+    tolerance = 0.01
+  )
+
 
 })
 
diff --git a/tests/testthat/test-Load.R b/tests/testthat/test-Load.R
new file mode 100644
index 0000000000000000000000000000000000000000..645ceaae964b81e6bf30591598ab273c0cf7881a
--- /dev/null
+++ b/tests/testthat/test-Load.R
@@ -0,0 +1,274 @@
+context("Testing Load with ECMWF System5c3s daily data")
+test_that("First time step is correctly interpreted:", {
+  system5 <- Load(var = 'prlr',
+                  exp = list(list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")),
+                  sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+                  nmember = 5,leadtimemin = 1,
+                  leadtimemax = NULL,
+                  storefreq = "daily", sampleperiod = 1,
+                  latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+                  output = 'lonlat', nprocs = 1,
+                  method = 'conservative', grid = 'r360x181',
+                  maskmod = vector("list", 15), maskobs = vector("list", 15),
+                  configfile = NULL, varmin = NULL, varmax = NULL,
+                  silent = FALSE, dimnames = NULL,
+                  remapcells = 2, path_glob_permissive = 'partial',
+                  nmemberobs = NULL, nleadtime = NULL)
+
+  
+  data <- Load(var = 'prlr',
+               exp = list(list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")),
+               obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 5,leadtimemin = 1,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+  
+  expect_equal(dim(data$mod), c(dataset = 1, member = 5, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$obs), c(dataset = 1, member = 1, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$mod)), 0)
+  expect_equal(sum(is.na(data$obs)), 0)
+  dates <- c(seq(as.POSIXct("1993-11-02", tz = 'UTC'), 
+                                       as.POSIXct("1994-06-03", tz = 'UTC'), "d"),
+                                   seq(as.POSIXct("1994-11-02", tz = 'UTC'),
+                                       as.POSIXct("1995-06-03", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$start, dates)
+  dates <- c(seq(as.POSIXct("1993-11-03", tz = 'UTC'),
+                 as.POSIXct("1994-06-04", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-03", tz = 'UTC'),
+                 as.POSIXct("1995-06-04", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$end, dates)
+  
+  expect_equal(data$mod, system5$mod)
+  
+  era5 <- Load(var = 'prlr', exp = NULL,
+              obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'),
+               nmember = 1, leadtimemin = 1,
+               leadtimemax = 215,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = 1, nleadtime = NULL)
+  # The values of observation when loading exp and obs simultaneously corresponds to the second day of the observations in the file:
+  expect_equal(data$obs[1,1,1,,,], era5$obs[1,1,1,2:215,,])
+  expect_equal(data$obs[1,1,2,,,], era5$obs[1,1,2,2:215,,]) 
+
+ # Test when 2 observational datasets are requested:
+  data <- Load(var = 'prlr',
+             exp = list(list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")),
+             obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc'), list(name = 'erainterim')),
+             sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+             nmember = 5, leadtimemin = 1,
+             leadtimemax = NULL,
+             storefreq = "daily", sampleperiod = 1,
+             latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+             output = 'lonlat', nprocs = 1,
+             method = 'conservative', grid = 'r360x181',
+             maskmod = vector("list", 15), maskobs = vector("list", 15),
+             configfile = NULL, varmin = NULL, varmax = NULL,
+             silent = FALSE, dimnames = NULL,
+             remapcells = 2, path_glob_permissive = 'partial',
+             nmemberobs = NULL, nleadtime = NULL)  
+  expect_equal(dim(data$obs), c(dataset = 2, member = 1, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$mod), c(dataset = 1, member = 5, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$mod)), 0)
+  expect_equal(sum(is.na(data$obs)), 0)
+  expect_equal(data$obs[1,1,,,,], era5$obs[1,1,,2:215,,])
+  expect_equal(system5$mod, data$mod)
+  dates <- c(seq(as.POSIXct("1993-11-02", tz = 'UTC'),
+                 as.POSIXct("1994-06-03", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-02", tz = 'UTC'),
+                 as.POSIXct("1995-06-03", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(dates, data$Dates$start)
+  dates <- c(seq(as.POSIXct("1993-11-03", tz = 'UTC'),
+                 as.POSIXct("1994-06-04", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-03", tz = 'UTC'),
+                 as.POSIXct("1995-06-04", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'  
+  expect_equal(dates, data$Dates$end)
+  
+  eraint <- Load(var = 'prlr', exp = NULL,
+              obs = list(list(name = 'erainterim')),
+              sdates = c('19931101', '19941101'),
+              nmember = 1, leadtimemin = 1,
+              leadtimemax = 215,
+              storefreq = "daily", sampleperiod = 1,
+              latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+              output = 'lonlat', nprocs = 1,
+              method = 'conservative', grid = 'r360x181',
+              maskmod = vector("list", 15), maskobs = vector("list", 15),
+              configfile = NULL, varmin = NULL, varmax = NULL,
+              silent = FALSE, dimnames = NULL,
+              remapcells = 2, path_glob_permissive = 'partial',
+              nmemberobs = 1, nleadtime = NULL)
+  expect_equal(data$obs[2,1,,,,], eraint$obs[1,1,,2:215,,])
+
+  # Test for 2 experimental datasets
+  data <- Load(var = 'prlr',
+               exp = list(list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc"),
+                          list(name = 'ecmwfS4', path = "/esarchive/exp/ecmwf/system4_m1/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")),
+               obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 5,leadtimemin = 1,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+  expect_equal(dim(data$obs), c(dataset = 1, member = 1, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$mod), c(dataset = 2, member = 5, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$obs)), 0)
+  expect_equal(sum(is.na(data$mod)), 0)
+  dates <- c(seq(as.POSIXct("1993-11-02", tz = 'UTC'),
+                                       as.POSIXct("1994-06-03", tz = 'UTC'), "d"),
+                                   seq(as.POSIXct("1994-11-02", tz = 'UTC'),
+                                       as.POSIXct("1995-06-03", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$start, dates)
+  dates <- c(seq(as.POSIXct("1993-11-03", tz = 'UTC'),
+                 as.POSIXct("1994-06-04", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-03", tz = 'UTC'),
+                 as.POSIXct("1995-06-04", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$end, dates)
+  expect_equal(data$obs[,,,,,], era5$obs[,,,2:215,,])
+  expect_equal(system5$exp, data$exp[1,,,,,])
+  
+  system4 <- Load(var = 'prlr',
+                  exp = list(list(name = 'ecmwfS4', path = "/esarchive/exp/ecmwf/system4_m1/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 1, leadtimemin = 1,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+
+  expect_equal(system4$exp, data$exp[2,,,,,])
+  expect_equal(system4$Dates$start[c(-1,-216)], data$Dates$start)
+  
+  
+  # Test for 2 experimental datasets oposite order
+  data <- Load(var = 'prlr',
+               exp = list(list(name = 'ecmwfS4', path = "/esarchive/exp/ecmwf/system4_m1/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc"), 
+                          list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")),
+               obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 5,leadtimemin = 1,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+  expect_equal(dim(data$mod), c(dataset = 2, member = 5, sdate = 2, ftime = 215,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$obs), c(dataset = 1, member = 1, sdate = 2, ftime = 215,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$mod[1,,,,,])), 0)
+  expect_equal(sum(is.na(data$mod[2,,,-215,,])),0)
+  expect_equal(sum(is.na(data$mod[2,,,215,,])), 640)
+  expect_equal(data$Dates$start, system4$Dates$start)
+  expect_equal(data$obs, era5$obs)
+  # Test leadtimemin > 1
+  data <- Load(var = 'prlr',
+               exp = list(list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc"),
+                          list(name = 'ecmwfS4', path = "/esarchive/exp/ecmwf/system4_m1/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")),
+               obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 5,leadtimemin = 2,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+  expect_equal(dim(data$mod), c(dataset = 2, member = 5, sdate = 2, ftime = 213,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$obs), c(dataset = 1, member = 1, sdate = 2,ftime = 213,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$mod)), 0)
+  expect_equal(sum(is.na(data$obs)), 0)
+  expect_equal(system4$mod[,1,,3:215,,], data$mod[2,1,,,,])
+  expect_equal(system5$mod[,,,2:214,,], data$mod[1,,,,,]) 
+  expect_equal(era5$obs[,,,3:215,,], data$obs[1,,,,,], tolerance = 0.001)
+  dates <- c(seq(as.POSIXct("1993-11-03", tz = 'UTC'),
+                 as.POSIXct("1994-06-03", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-03", tz = 'UTC'),
+                 as.POSIXct("1995-06-03", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$start, dates)
+  
+  data <- Load(var = 'prlr',
+               exp = list(list(name = 'ecmwfS4', path = "/esarchive/exp/ecmwf/system4_m1/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc"),
+                          list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")),
+               obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 5,leadtimemin = 2,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+  expect_equal(system4$mod[1,1,,2:215,,], data$mod[1,1,,1:214,,])
+  expect_equal(system5$mod[1,1,,1:214,,], data$mod[2,1,,1:214,,])
+  expect_equal(dim(data$obs), c(dataset = 1, member = 1, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$mod), c(dataset = 2, member = 5, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$obs)), 0)
+  expect_equal(sum(is.na(data$mod)), 0)
+  expect_equal(data$obs[,,,1,,], era5$obs[,,,2,,])
+  dates <- c(seq(as.POSIXct("1993-11-02", tz = 'UTC'),
+                 as.POSIXct("1994-06-03", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-02", tz = 'UTC'),
+                 as.POSIXct("1995-06-03", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$start, dates)
+})