diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000000000000000000000000000000000000..dff3121316a8d115d43408de3af7a904c4d2dede
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,12 @@
+Package: csdownscale
+Title: What the Package Does (One Line, Title Case)
+Version: 0.0.0.9000
+Authors@R: 
+    person("First", "Last", , "first.last@example.com", role = c("aut", "cre"),
+           comment = c(ORCID = "YOUR-ORCID-ID"))
+Description: What the package does (one paragraph).
+License: `use_mit_license()`, `use_gpl3_license()` or friends to pick a
+    license
+Encoding: UTF-8
+Roxygen: list(markdown = TRUE)
+RoxygenNote: 7.2.0
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000000000000000000000000000000000000..1ae93bd4a38fbc20a4098bdbb1d55852fe23e8c6
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,20 @@
+# Generated by roxygen2: do not edit by hand
+
+export(Analogs)
+export(CST_Analogs)
+export(CST_Intbc)
+export(CST_Interpolation)
+export(CST_Intlr)
+export(CST_LogisticReg)
+export(Intbc)
+export(Interpolation)
+export(Intlr)
+export(LogisticReg)
+import(CSTools)
+import(multiApply)
+import(nnet)
+import(plyr)
+importFrom(FNN,get.knnx)
+importFrom(laply,plyr)
+importFrom(s2dv,CDORemap)
+importFrom(s2dv,InsertDim)
diff --git a/R/Analogs.R b/R/Analogs.R
index 4e85d22a86beb37cf7e66f98cfd5a453f4c220c9..fd6e2f9776259d87745b7aeaaaf5554d0e3907cf 100644
--- a/R/Analogs.R
+++ b/R/Analogs.R
@@ -106,12 +106,17 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo
                  region = region, return_indices = return_indices, loocv_window = loocv_window, 
                  ncores = ncores)
 
-  res_s2dv <- list(exp = suppressWarnings(s2dv_cube(data = res$data, lon = res$lon, lat = res$lat)),
-                   obs = suppressWarnings(s2dv_cube(data = res$obs, lon = res$lon, lat = res$lat)))
+  # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
+  exp$data <- res$data
+  exp$lon <- res$lon
+  exp$lat <- res$lat
 
+  obs$data <- res$obs
+  obs$lat <- res$lat
+  obs$lon <- res$lon
 
+  res_s2dv <- list(exp = exp, obs = obs)
   return(res_s2dv)
-
 }
 
 #'@rdname Analogs
diff --git a/R/Intbc.R b/R/Intbc.R
index a819b76a729831001e3acb2275be02311232c6a0..86bb5a9c340afce14399417f3f656af872707eeb 100644
--- a/R/Intbc.R
+++ b/R/Intbc.R
@@ -84,10 +84,18 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point
                method_point_interp = method_point_interp, lat_dim = lat_dim, lon_dim = lon_dim, 
                sdate_dim = sdate_dim, member_dim = member_dim, region = region, ncores = ncores)
 
-  res_s2dv <- list(exp = suppressWarnings(s2dv_cube(data = res$data, lon = res$lon, lat = res$lat)),
-                   obs = suppressWarnings(s2dv_cube(data = res$obs, lon = res$lon, lat = res$lat)))
- 
+  # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
+  exp$data <- res$data
+  exp$lon <- res$lon
+  exp$lat <- res$lat
+  
+  obs$data <- res$obs
+  obs$lat <- res$lat
+  obs$lon <- res$lon
+
+  res_s2dv <- list(exp = exp, obs = obs)
   return(res_s2dv)
+  
 }
 
 #'@rdname Intbc
diff --git a/R/Interpolation.R b/R/Interpolation.R
index 7a81d8b76ecebd98bcce8240f0641bf11bd81cfb..e594691fc148cefdcba35c92044d5539e698de2a 100644
--- a/R/Interpolation.R
+++ b/R/Interpolation.R
@@ -72,9 +72,12 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr
                        method_remap = method_remap, target_grid = target_grid, lat_dim = lat_dim,
                        lon_dim = lon_dim, region = region, method_point_interp = method_point_interp)
 
-  res_s2dv <- list(exp = suppressWarnings(s2dv_cube(data = res$data, lon = res$lon, lat = res$lat)),
-                   obs  = NULL)
+  # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
+  exp$data <- res$data
+  exp$lon <- res$lon
+  exp$lat <- res$lat
 
+  res_s2dv <- list(exp = exp, obs  = NULL)
   return(res_s2dv)
 }
 
@@ -537,7 +540,7 @@ create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
       if(substr(method,2,2) == "E" & westtoeast == T) {seli <- ceiling(gridpoints$i) }
       if(substr(method,2,2) == "W" & westtoeast == T) {seli <- floor(gridpoints$i) }
       if(substr(method,2,2) == "E" & westtoeast == F) {seli <- floor(gridpoints$i) }
-      if(substr(method,2,2) == "W" & westtoeast == F) {seli <- ceil(gridpoints$i) }
+      if(substr(method,2,2) == "W" & westtoeast == F) {seli <- ceiling(gridpoints$i) }
       
       seli[seli == griddes$xsize + 1] <- 1  # close the longitudes
       
diff --git a/R/Intlr.R b/R/Intlr.R
index bf744b36bb8572011d706f9d850c64385119ab37..565e30468aeb0b5fcc7fe6eabc8c0d18bb7dff9a 100644
--- a/R/Intlr.R
+++ b/R/Intlr.R
@@ -107,14 +107,16 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in
                member_dim = member_dim, large_scale_predictor_dimname = large_scale_predictor_dimname, 
                loocv = loocv, region = region, ncores = ncores) 
 
-  #-----------------------------------
-  # Create an s2dv_cube object
-  #-----------------------------------
-  res_s2dv <- list(exp = suppressWarnings(s2dv_cube(data = res$data, lon = res$lon, lat = res$lat)),
-                   obs = suppressWarnings(s2dv_cube(data = res$obs, lon = res$lon, lat = res$lat)))
-
+  # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
+  exp$data <- res$data
+  exp$lon <- res$lon
+  exp$lat <- res$lat
 
+  obs$data <- res$obs
+  obs$lat <- res$lat
+  obs$lon <- res$lon
 
+  res_s2dv <- list(exp = exp, obs = obs)
   return(res_s2dv)
 }
 
diff --git a/R/LogisticReg.R b/R/LogisticReg.R
index d46fb6a396f8578328b6ce18c91e68b822d2880e..24be693603320058f1da85b60d9e1fc9c56d694b 100644
--- a/R/LogisticReg.R
+++ b/R/LogisticReg.R
@@ -110,13 +110,16 @@ CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_me
                      source_file = exp$source_files[1], region = region, loocv = loocv, 
                      ncores = ncores)
 
-  #-----------------------------------
-  # Create an s2dv_cube object
-  #-----------------------------------
-  res_s2dv <- list(exp = suppressWarnings(s2dv_cube(data = res$data, lon = res$lon, lat = res$lat)),
-                   obs = suppressWarnings(s2dv_cube(data = res$obs, lon = res$lon, lat = res$lat)))
+  # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
+  exp$data <- res$data
+  exp$lon <- res$lon
+  exp$lat <- res$lat
 
+  obs$data <- res$obs
+  obs$lat <- res$lat
+  obs$lon <- res$lon
 
+  res_s2dv <- list(exp = exp, obs = obs)
   return(res_s2dv)
 }
 
diff --git a/inst/examples/interpolation.R b/inst/examples/interpolation.R
index 96c446b04969e3503cea749e7a9bd498a4da35a7..59330b6aa12f82b26686093c951243c65b02cc67 100644
--- a/inst/examples/interpolation.R
+++ b/inst/examples/interpolation.R
@@ -24,11 +24,14 @@ sdates <- c('20000201','20010201','20020201','20030201','20040201','20050201','2
 
 exp1 <- startR::Start(dat = '/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc',
                      var = 'tas', time = indices(1), member = indices(1:3), sdate = sdates,
-                     lat = values(list(latmin, latmax)), lat_reorder = Sort(decreasing = FALSE),
+                     lat = values(list(latmin, latmax)), 
+                     lat_reorder = Sort(decreasing = FALSE),
                      lon = values(list(lonmin, lonmax)), 
                      lon_reorder = CircularSort(-180, 180),
-                     synonims = list(var = c('var','variable'), lon = c('lon', 'longitude'),
-                     lat = c('lat', 'latitude'), member = c('member','ensemble')),
+                     synonims = list(var = c('var','variable'), 
+                                     lon = c('lon', 'longitude'),
+                                     lat = c('lat', 'latitude'), 
+                                     member = c('member','ensemble')),
                      return_vars = list(lat = 'dat', lon = 'dat'),
                      num_procs = 1, retrieve = TRUE)
 exp2 <- s2dv_cube(exp1, lat = attr(exp1, "Variables")$dat1$lat, lon = attr(exp1, "Variables")$dat1$lon,
diff --git a/man/Analogs.Rd b/man/Analogs.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..644aeeb870f24f4f19135ca2ab5e546510db6272
--- /dev/null
+++ b/man/Analogs.Rd
@@ -0,0 +1,155 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Analogs.R
+\name{Analogs}
+\alias{Analogs}
+\title{Downscaling using Analogs based on large scale fields.}
+\usage{
+Analogs(
+  exp,
+  obs,
+  exp_lats,
+  exp_lons,
+  obs_lats,
+  obs_lons,
+  grid_exp,
+  obs2 = NULL,
+  obs2_lats = NULL,
+  obs2_lons = NULL,
+  nanalogs = 3,
+  fun_analog = NULL,
+  lat_dim = "lat",
+  lon_dim = "lon",
+  sdate_dim = "sdate",
+  time_dim = "time",
+  member_dim = "member",
+  region = NULL,
+  return_indices = FALSE,
+  loocv_window = TRUE,
+  ncores = 1
+)
+}
+\arguments{
+\item{exp}{an array with named dimensions containing the experimental field on the
+coarse scale for which the downscaling is aimed. The object must have, at least,
+the dimensions latitude, longitude, start date and time. The object is expected to be
+already subset for the desired region. Data can be in one or two integrated regions, e.g.,
+crossing the Greenwich meridian. To get the correct results in the latter case,
+the borders of the region should be specified in the parameter 'region'. See parameter
+'region'.}
+
+\item{obs}{an array with named dimensions containing the observational field. The object
+must have, at least, the dimensions latitude, longitude, start date and time. The object
+is expected to be already subset for the desired region. Optionally, 'obs' can have the
+dimension 'window', containing the sampled fields into which the function will look for
+the analogs. See function 'generate_window()'. Otherwise, the function will look for
+analogs using all the possible fields contained in obs.}
+
+\item{exp_lats}{a numeric vector containing the latitude values in 'exp'. Latitudes must
+range from -90 to 90.}
+
+\item{exp_lons}{a numeric vector containing the longitude values in 'exp'. Longitudes
+can range from -180 to 180 or from 0 to 360.}
+
+\item{obs_lats}{a numeric vector containing the latitude values in 'obs'. Latitudes must
+range from -90 to 90.}
+
+\item{obs_lons}{a numeric vector containing the longitude values in 'obs'. Longitudes
+can range from -180 to 180 or from 0 to 360.}
+
+\item{grid_exp}{a character vector with a path to an example file of the exp data.
+It can be either a path to another NetCDF file which to read the target grid from
+(a single grid must be defined in such file) or a character vector indicating the
+coarse grid to be passed to CDO, and it must be a grid recognised by CDO.}
+
+\item{obs2}{an 's2dv' object with named dimensions containing a different observational
+field to that in 'obs'. If provided, these observations will be used in the training,
+i.e. the searching for analogs, so that they should be in a coarser grid to those in
+'obs'. Training with observations on a grid with a spatial resolution closer to that
+in 'exp', will in principle ensure better results. The object must have, at least, the
+dimensions latitude, longitude and start date. The object is expected to be already
+subset for the desired region.}
+
+\item{nanalogs}{an integer indicating the number of analogs to be searched.}
+
+\item{fun_analog}{a function to be applied over the found analogs. Only these options
+are valid: "mean", "wmean", "max", "min", "median" or NULL. If set to NULL (default),
+the function returns the found analogs.}
+
+\item{lat_dim}{a character vector indicating the latitude dimension name in the element
+'data' in exp and obs. Default set to "lat".}
+
+\item{lon_dim}{a character vector indicating the longitude dimension name in the element
+'data' in exp and obs. Default set to "lon".}
+
+\item{sdate_dim}{a character vector indicating the start date dimension name in the
+element 'data' in exp and obs. Default set to "sdate".}
+
+\item{time_dim}{a character vector indicating the time dimension name in the element
+'data' in exp and obs. Default set to "time".}
+
+\item{member_dim}{a character vector indicating the member dimension name in the element
+'data' in exp and obs. Default set to "member".}
+
+\item{region}{a numeric vector indicating the borders of the region defined in exp.
+It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers
+to the left border, while lonmax refers to the right border. latmin indicates the lower
+border, whereas latmax indicates the upper border. If set to NULL (default), the function
+takes the first and last elements of the latitudes and longitudes.}
+
+\item{return_indices}{a logical vector indicating whether to return the indices of the
+analogs together with the downscaled fields. The indices refer to the position of the
+element in the vector time * start_date. If 'obs' contain the dimension 'window', it will
+refer to the position of the element in the dimension 'window'.  Default to FALSE.}
+
+\item{loocv_window}{a logical vector only to be used if 'obs' does not have the dimension
+'window'. It indicates whether to apply leave-one-out cross-validation in the creation
+of the window. It is recommended to be set to TRUE. Default to TRUE.}
+
+\item{ncores}{an integer indicating the number of cores to use in parallel computation.}
+}
+\value{
+A list of three elements. 'data' contains the dowscaled field, 'lat' the
+downscaled latitudes, and 'lon' the downscaled longitudes. If fun_analog is set to NULL
+(default), the output array in 'data' also contains the dimension 'analog' with the best
+analog days.
+}
+\description{
+This function performs a downscaling using Analogs. To compute
+the analogs given a coarse-scale field, the function looks for days with
+similar conditions in the historical observations. The coarse scale and
+observation data can be either global or regional. In the latter case, the
+region is defined by the user. In principle, the coarse and observation data
+should be of the same variable, although different variables can also be admitted.
+The analogs function will find the N best analogs based in Minimum Euclidean
+distance.
+
+The search of analogs must be done in the longest dataset posible, but might
+require high-memory computational resources. This is important since it is
+necessary to have a good representation of the possible states of the field in
+the past, and therefore, to get better analogs. The function can also look for
+analogs within a window of D days, but is the user who has to define that window.
+Otherwise, the function will look for analogs in the whole dataset. This function
+is intended to downscale climate prediction data (i.e., sub-seasonal, seasonal
+and decadal predictions) but can admit climate projections or reanalyses. It does
+not have constrains of specific region or variables to downscale.
+}
+\examples{
+exp <- rnorm(15000)
+dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 5, time  = 30)
+exp_lons <- 1:5
+exp_lats <- 1:4
+obs <- rnorm(27000)
+dim(obs) <- c(lat = 12, lon = 15, sdate = 5, time  = 30)
+obs_lons <- seq(0,6, 6/14)
+obs_lats <- seq(0,6, 6/11)
+downscaled_field <- Analogs(exp = exp, obs = obs, exp_lats = exp_lats, exp_lons = exp_lons, 
+obs_lats = obs_lats, obs_lons = obs_lons, grid_exp = 'r360x180')
+}
+\seealso{
+\code{\link[s2dverification]{CDORemap}}
+}
+\author{
+J. Ramon, \email{jaume.ramon@bsc.es}
+
+Ll. Lledó, \email{llorenc.lledo@ecmwf.int}
+}
diff --git a/man/CST_Analogs.Rd b/man/CST_Analogs.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..64fc4ffe42a367a87d9acf133955704d39f3e8e7
--- /dev/null
+++ b/man/CST_Analogs.Rd
@@ -0,0 +1,128 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Analogs.R
+\name{CST_Analogs}
+\alias{CST_Analogs}
+\title{Downscaling using Analogs based on large scale fields.}
+\usage{
+CST_Analogs(
+  exp,
+  obs,
+  grid_exp,
+  obs2 = NULL,
+  nanalogs = 3,
+  fun_analog = NULL,
+  lat_dim = "lat",
+  lon_dim = "lon",
+  sdate_dim = "sdate",
+  time_dim = "time",
+  member_dim = "member",
+  region = NULL,
+  return_indices = FALSE,
+  loocv_window = TRUE,
+  ncores = 1
+)
+}
+\arguments{
+\item{exp}{an 's2dv' object with named dimensions containing the experimental field on
+the coarse scale for which the downscaling is aimed. The object must have, at least,
+the dimensions latitude, longitude, start date and time. The object is expected to be
+already subset for the desired region. Data can be in one or two integrated regions, e.g.,
+crossing the Greenwich meridian. To get the correct results in the latter case,
+the borders of the region should be specified in the parameter 'region'. See parameter
+'region'.}
+
+\item{obs}{an 's2dv' object with named dimensions containing the observational field.
+The object must have, at least, the dimensions latitude, longitude and start date.
+The object is expected to be already subset for the desired region.}
+
+\item{grid_exp}{a character vector with a path to an example file of the exp data.
+It can be either a path to another NetCDF file which to read the target grid from
+(a single grid must be defined in such file) or a character vector indicating the
+coarse grid to be passed to CDO, and it must be a grid recognised by CDO.}
+
+\item{obs2}{an 's2dv' object with named dimensions containing a different observational
+field to that in 'obs'. If provided, these observations will be used in the training,
+i.e. the searching for analogs, so that they should be in a coarser grid to those in
+'obs'. Training with observations on a grid with a spatial resolution closer to that
+in 'exp', will in principle ensure better results. The object must have, at least, the
+dimensions latitude, longitude and start date. The object is expected to be already
+subset for the desired region.}
+
+\item{nanalogs}{an integer indicating the number of analogs to be searched}
+
+\item{fun_analog}{a function to be applied over the found analogs. Only these options
+are valid: "mean", "wmean", "max", "min", "median" or NULL. If set to NULL (default),
+the function returns the found analogs.}
+
+\item{lat_dim}{a character vector indicating the latitude dimension name in the element
+'data' in exp and obs. Default set to "lat".}
+
+\item{lon_dim}{a character vector indicating the longitude dimension name in the element
+'data' in exp and obs. Default set to "lon".}
+
+\item{sdate_dim}{a character vector indicating the start date dimension name in the
+element 'data' in exp and obs. Default set to "sdate".}
+
+\item{time_dim}{a character vector indicating the time dimension name in the element
+'data' in exp and obs. Default set to "time".}
+
+\item{member_dim}{a character vector indicating the member dimension name in the element
+'data' in exp and obs. Default set to "member".}
+
+\item{region}{a numeric vector indicating the borders of the region defined in exp.
+It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers
+to the left border, while lonmax refers to the right border. latmin indicates the lower
+border, whereas latmax indicates the upper border. If set to NULL (default), the function
+takes the first and last elements of the latitudes and longitudes.}
+
+\item{return_indices}{a logical vector indicating whether to return the indices of the
+analogs together with the downscaled fields. Default to FALSE.}
+
+\item{loocv_window}{a logical vector only to be used if 'obs' does not have the dimension
+'window'. It indicates whether to apply leave-one-out cross-validation in the creation
+of the window. It is recommended to be set to TRUE. Default to TRUE.}
+
+\item{ncores}{an integer indicating the number of cores to use in parallel computation.}
+}
+\value{
+An 's2dv' object. The element 'data' contains the dowscaled field, 'lat' the
+downscaled latitudes, and 'lon' the downscaled longitudes. If fun_analog is set to NULL
+(default), the output array in 'data' also contains the dimension 'analog' with the best
+analog days.
+}
+\description{
+This function performs a downscaling using Analogs. To compute
+the analogs given a coarse-scale field, the function looks for days with
+similar conditions in the historical observations. The coarse scale and
+observation data can be either global or regional. In the latter case, the
+region is defined by the user. In principle, the coarse and observation data
+should be of the same variable, although different variables can also be admitted.
+The analogs function will find the N best analogs based in Minimum Euclidean
+distance.
+
+The search of analogs must be done in the longest dataset posible, but might
+require high-memory computational resources. This is important since it is
+necessary to have a good representation of the possible states of the field in
+the past, and therefore, to get better analogs. The function can also look for
+analogs within a window of D days, but is the user who has to define that window.
+Otherwise, the function will look for analogs in the whole dataset. This function
+is intended to downscale climate prediction data (i.e., sub-seasonal, seasonal
+and decadal predictions) but can admit climate projections or reanalyses. It does
+not have constrains of specific region or variables to downscale.
+}
+\examples{
+exp <- rnorm(15000)
+dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 5, time  = 30)
+exp_lons <- 1:5
+exp_lats <- 1:4
+obs <- rnorm(27000)
+dim(obs) <- c(lat = 12, lon = 15, sdate = 5, time  = 30)
+obs_lons <- seq(0,6, 6/14)
+obs_lats <- seq(0,6, 6/11)
+exp <- s2dv_cube(data = exp, lat = exp_lats, lon = exp_lons)
+obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons)
+downscaled_field <- CST_Analogs(exp = exp, obs = obs, grid_exp = 'r360x180')
+}
+\author{
+J. Ramon, \email{jaume.ramon@bsc.es}
+}
diff --git a/man/CST_Intbc.Rd b/man/CST_Intbc.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..3bc4d1ae1e0370a6ffe03005674e24733ac3813d
--- /dev/null
+++ b/man/CST_Intbc.Rd
@@ -0,0 +1,104 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Intbc.R
+\name{CST_Intbc}
+\alias{CST_Intbc}
+\title{Downscaling using interpolation and bias adjustment.}
+\usage{
+CST_Intbc(
+  exp,
+  obs,
+  target_grid,
+  bc_method,
+  int_method = NULL,
+  points = NULL,
+  method_point_interp = NULL,
+  lat_dim = "lat",
+  lon_dim = "lon",
+  sdate_dim = "sdate",
+  member_dim = "member",
+  region = NULL,
+  ncores = 1
+)
+}
+\arguments{
+\item{exp}{an 's2dv object' containing the experimental field on the
+coarse scale for which the downscaling is aimed. The object must have, at least,
+the dimensions latitude, longitude, start date and member. The object is expected to be
+already subset for the desired region. Data can be in one or two integrated regions, e.g.,
+crossing the Greenwich meridian. To get the correct results in the latter case,
+the borders of the region should be specified in the parameter 'region'. See parameter
+'region'.}
+
+\item{obs}{an 's2dv object' containing the observational field. The object
+must have, at least, the dimensions latitude, longitude and start date. The object is
+expected to be already subset for the desired region.}
+
+\item{target_grid}{a character vector indicating the target grid to be passed to CDO.
+It must be a grid recognised by CDO or a NetCDF file.}
+
+\item{bc_method}{a character vector indicating the bias adjustment method to be applied after
+the interpolation. Accepted methods are 'simple_bias', 'calibration', 'quantile_mapping'. The
+abbreviations 'sbc', 'cal', 'qm' can also be used.}
+
+\item{int_method}{a character vector indicating the regridding method to be passed to CDORemap.
+Accepted methods are "con", "bil", "bic", "nn", "con2". If "nn" method is to be used, CDO_1.9.8
+or newer version is required.}
+
+\item{points}{a list of two elements containing the point latitudes and longitudes
+of the locations to downscale the model data. The list must contain the two elements
+named as indicated in the parameters 'lat_dim' and 'lon_dim'. If the downscaling is
+to a point location, only regular grids are allowed for exp and obs. Only needed if the
+downscaling is to a point location.}
+
+\item{method_point_interp}{a character vector indicating the interpolation method to interpolate
+model gridded data into the point locations. Accepted methods are "nearest", "bilinear", "9point",
+"invdist4nn", "NE", "NW", "SE", "SW". Only needed if the downscaling is to a point location.}
+
+\item{lat_dim}{a character vector indicating the latitude dimension name in the element 'data'
+in exp and obs. Default set to "lat".}
+
+\item{lon_dim}{a character vector indicating the longitude dimension name in the element 'data'
+in exp and obs. Default set to "lon".}
+
+\item{sdate_dim}{a character vector indicating the start date dimension name in the element
+'data' in exp and obs. Default set to "sdate".}
+
+\item{member_dim}{a character vector indicating the member dimension name in the element
+'data' in exp and obs. Default set to "member".}
+
+\item{region}{a numeric vector indicating the borders of the region defined in obs.
+It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers
+to the left border, while lonmax refers to the right border. latmin indicates the lower
+border, whereas latmax indicates the upper border. If set to NULL (default), the function
+takes the first and last elements of the latitudes and longitudes.}
+
+\item{ncores}{an integer indicating the number of cores to use in parallel computation.}
+}
+\value{
+An 's2dv' object. The element 'data' contains the dowscaled field, 'lat' the
+downscaled latitudes, and 'lon' the downscaled longitudes.
+}
+\description{
+This function performs a downscaling using an interpolation and a later bias
+adjustment. It is recommended that the observations are passed already in the target grid.
+Otherwise, the function will also perform an interpolation of the observed field into the
+target grid. The coarse scale and observation data can be either global or regional. In the
+latter case, the region is defined by the user. In principle, the coarse and observation data
+are intended to be of the same variable, although different variables can also be admitted.
+}
+\examples{
+exp <- rnorm(500)
+dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 5)
+exp_lons <- 1:5
+exp_lats <- 1:4
+obs <- rnorm(900)
+dim(obs) <- c(lat = 12, lon = 15, sdate = 5)
+obs_lons <- seq(1,5, 4/14)
+obs_lats <- seq(1,4, 3/11)
+exp <- s2dv_cube(data = exp, lat = exp_lats, lon = exp_lons)
+obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons)
+res <- CST_Intbc(exp = exp, obs = obs, target_grid = 'r1280x640', bc_method = 'simple_bias', int_method = 'conservative')
+}
+\author{
+J. Ramon, \email{jaume.ramon@bsc.es}
+}
diff --git a/man/CST_Interpolation.Rd b/man/CST_Interpolation.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..1b8b0cb0a7d388f935d0e4e5d9ffd04269a16cf0
--- /dev/null
+++ b/man/CST_Interpolation.Rd
@@ -0,0 +1,79 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Interpolation.R
+\name{CST_Interpolation}
+\alias{CST_Interpolation}
+\title{Regrid or interpolate gridded data to a point location.}
+\usage{
+CST_Interpolation(
+  exp,
+  points = NULL,
+  method_remap = NULL,
+  target_grid = NULL,
+  lat_dim = "lat",
+  lon_dim = "lon",
+  region = NULL,
+  method_point_interp = NULL
+)
+}
+\arguments{
+\item{exp}{s2dv object containing the experimental field on the
+coarse scale for which the downscaling is aimed. The object must have, at least,
+the dimensions latitude and longitude. The field data is expected to be already subset
+for the desired region. Data can be in one or two integrated regions, e.g.,
+crossing the Greenwich meridian. To get the correct results in the latter case,
+the borders of the region should be specified in the parameter 'region'. See parameter
+'region'.}
+
+\item{points}{a list of two elements containing the point latitudes and longitudes
+of the locations to downscale the model data. The list must contain the two elements
+named as indicated in the parameters 'lat_dim' and 'lon_dim'. If the downscaling is
+to a point location, only regular grids are allowed for exp and obs. Only needed if the
+downscaling is to a point location.}
+
+\item{method_remap}{a character vector indicating the regridding method to be passed
+to CDORemap. Accepted methods are "con", "bil", "bic", "nn", "con2". If "nn" method is
+to be used, CDO_1.9.8 or newer version is required.}
+
+\item{target_grid}{a character vector indicating the target grid to be passed to CDO.
+It must be a grid recognised by CDO or a NetCDF file.}
+
+\item{lat_dim}{a character vector indicating the latitude dimension name in the element
+'exp' and/or 'points'. Default set to "lat".}
+
+\item{lon_dim}{a character vector indicating the longitude dimension name in the element
+'exp' and/or 'points'. Default set to "lon".}
+
+\item{region}{a numeric vector indicating the borders of the region defined in exp.
+It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers
+to the left border, while lonmax refers to the right border. latmin indicates the lower
+border, whereas latmax indicates the upper border. If set to NULL (default), the function
+takes the first and last elements of the latitudes and longitudes.}
+
+\item{method_point_interp}{a character vector indicating the interpolation method to
+interpolate model gridded data into the point locations. Accepted methods are "nearest",
+"bilinear", "9point", "invdist4nn", "NE", "NW", "SE", "SW".}
+}
+\value{
+An s2dv object containing the dowscaled field.
+}
+\description{
+This function interpolates gridded model data from one grid to
+another (regrid) or interpolates gridded model data to a set of point locations.
+The gridded model data can be either global or regional. In the latter case, the
+region is defined by the user. It does not have constrains of specific region or
+variables to downscale.
+}
+\examples{
+exp <- rnorm(500)
+dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 5, time  = 1)
+lons <- 1:5
+lats <- 1:4
+exp <- s2dv_cube(data = exp, lat = lats, lon = lons) 
+res <- CST_Interpolation(exp = exp, method_remap = 'conservative', target_grid = 'r1280x640')
+}
+\seealso{
+\code{\link[s2dverification]{CDORemap}}
+}
+\author{
+J. Ramon, \email{jaume.ramon@bsc.es}
+}
diff --git a/man/CST_Intlr.Rd b/man/CST_Intlr.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..1a59b1204fdaf54e234004565b84009e66f1ab47
--- /dev/null
+++ b/man/CST_Intlr.Rd
@@ -0,0 +1,129 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Intlr.R
+\name{CST_Intlr}
+\alias{CST_Intlr}
+\title{Downscaling using interpolation and linear regression.}
+\usage{
+CST_Intlr(
+  exp,
+  obs,
+  lr_method,
+  target_grid = NULL,
+  points = NULL,
+  int_method = NULL,
+  method_point_interp = NULL,
+  predictors = NULL,
+  lat_dim = "lat",
+  lon_dim = "lon",
+  sdate_dim = "sdate",
+  time_dim = "time",
+  member_dim = "member",
+  large_scale_predictor_dimname = "vars",
+  loocv = FALSE,
+  region = NULL,
+  ncores = 1
+)
+}
+\arguments{
+\item{exp}{an 's2dv object' containing the experimental field on the
+coarse scale for which the downscaling is aimed. The object must have, at least,
+the dimensions latitude, longitude, start date and member. The object is expected to be
+already subset for the desired region. Data can be in one or two integrated regions, e.g.,
+crossing the Greenwich meridian. To get the correct results in the latter case,
+the borders of the region should be specified in the parameter 'region'. See parameter
+'region'.}
+
+\item{obs}{an 's2dv object' containing the observational field. The object
+must have, at least, the dimensions latitude, longitude and start date. The object is
+expected to be already subset for the desired region.}
+
+\item{lr_method}{a character vector indicating the linear regression method to be applied
+after the interpolation. Accepted methods are 'basic', 'large-scale' and '4nn'. The 'basic'
+method fits a linear regression using high resolution observations as predictands and the
+interpolated model data as predictor. Then, the regression equation is to the interpolated
+model data to correct the interpolated values. The 'large-scale' method fits a linear
+regression with large-scale predictors from the same model (e.g. teleconnection indices)
+as predictors and high-resolution observations as predictands. This equation is then
+applied to the interpolated model values. Finally, the '4nn' method uses a linear
+regression with the four nearest neighbours as predictors and high-resolution observations
+as predictands. It is then applied to model data to correct the interpolated values.}
+
+\item{target_grid}{a character vector indicating the target grid to be passed to CDO.
+It must be a grid recognised by CDO or a NetCDF file.}
+
+\item{points}{a list of two elements containing the point latitudes and longitudes
+of the locations to downscale the model data. The list must contain the two elements
+named as indicated in the parameters 'lat_dim' and 'lon_dim'. If the downscaling is
+to a point location, only regular grids are allowed for exp and obs. Only needed if the
+downscaling is to a point location.}
+
+\item{int_method}{a character vector indicating the regridding method to be passed
+to CDORemap. Accepted methods are "con", "bil", "bic", "nn", "con2". If "nn" method is
+to be used, CDO_1.9.8 or newer version is required.}
+
+\item{method_point_interp}{a character vector indicating the interpolation method to
+interpolate model gridded data into the point locations. Accepted methods are "nearest",
+"bilinear", "9point", "invdist4nn", "NE", "NW", "SE", "SW".}
+
+\item{predictors}{an array with large-scale data to be used in the 'large-scale' method.
+Only needed if the linear regression method is set to 'large-scale'. It must have, at
+least the dimension start date and another dimension whose name has to be specified in
+the parameter 'large_scale_predictor_dimname'. It should contain as many elements as the
+number of large-scale predictors.}
+
+\item{lat_dim}{a character vector indicating the latitude dimension name in the element 'data'
+in exp and obs. Default set to "lat".}
+
+\item{lon_dim}{a character vector indicating the longitude dimension name in the element 'data'
+in exp and obs. Default set to "lon".}
+
+\item{sdate_dim}{a character vector indicating the start date dimension name in the element
+'data' in exp and obs. Default set to "sdate".}
+
+\item{time_dim}{a character vector indicating the time dimension name in the element
+'data' in exp and obs. Default set to "time".}
+
+\item{large_scale_predictor_dimname}{a character vector indicating the name of the
+dimension in 'predictors' that contain the predictor variables. See parameter 'predictors'.}
+
+\item{loocv}{a logical indicating whether to apply leave-one-out cross-validation when
+generating the linear regressions. Default to FALSE.}
+
+\item{region}{a numeric vector indicating the borders of the region defined in exp.
+It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers
+to the left border, while lonmax refers to the right border. latmin indicates the lower
+border, whereas latmax indicates the upper border. If set to NULL (default), the function
+takes the first and last elements of the latitudes and longitudes.}
+
+\item{ncores}{an integer indicating the number of cores to use in parallel computation.}
+}
+\value{
+A list of three elements. 'data' contains the dowscaled field, 'lat' the
+downscaled latitudes, and 'lon' the downscaled longitudes.
+}
+\description{
+This function performs a downscaling using an interpolation and a linear
+regression. Different methodologies that employ linear regressions are available. See
+parameter 'lr_method' for more information. It is recommended that the observations
+are passed already in the target grid. Otherwise, the function will also perform an
+interpolation of the observed field into the target grid. The coarse scale and
+observation data can be either global or regional. In the latter case, the region is
+defined by the user. In principle, the coarse and observation data are intended to
+be of the same variable, although different variables can also be admitted.
+}
+\examples{
+exp <- rnorm(500) 
+dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 5) 
+exp_lons <- 1:5 
+exp_lats <- 1:4 
+obs <- rnorm(900) 
+dim(obs) <- c(lat = 12, lon = 15, sdate = 5) 
+obs_lons <- seq(1,5, 4/14) 
+obs_lats <- seq(1,4, 3/11) 
+exp <- s2dv_cube(data = exp, lat = exp_lats, lon = exp_lons)
+obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons)
+res <- CST_Intlr(exp = exp, obs = obs, target_grid = 'r1280x640', lr_method = 'basic', int_method = 'conservative')
+}
+\author{
+J. Ramon, \email{jaume.ramon@bsc.es}
+}
diff --git a/man/CST_LogisticReg.Rd b/man/CST_LogisticReg.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..5ddcae02977f4523e5f4618363c5de7eb8b3c9a7
--- /dev/null
+++ b/man/CST_LogisticReg.Rd
@@ -0,0 +1,129 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/LogisticReg.R
+\name{CST_LogisticReg}
+\alias{CST_LogisticReg}
+\title{Downscaling using interpolation and logistic regression.}
+\usage{
+CST_LogisticReg(
+  exp,
+  obs,
+  target_grid,
+  int_method = NULL,
+  log_reg_method = "ens_mean",
+  probs_cat = c(1/3, 2/3),
+  return_most_likely_cat = FALSE,
+  points = NULL,
+  method_point_interp = NULL,
+  lat_dim = "lat",
+  lon_dim = "lon",
+  sdate_dim = "sdate",
+  member_dim = "member",
+  region = NULL,
+  loocv = FALSE,
+  ncores = 1
+)
+}
+\arguments{
+\item{exp}{an 's2dv object' with named dimensions containing the experimental field on
+the coarse scale for which the downscaling is aimed. The object must have, at least,
+the dimensions latitude, longitude, start date and member. The object is expected to be
+already subset for the desired region. Data can be in one or two integrated regions, e.g.,
+crossing the Greenwich meridian. To get the correct results in the latter case,
+the borders of the region should be specified in the parameter 'region'. See parameter
+'region'.}
+
+\item{obs}{an 's2dv object' with named dimensions containing the observational field.
+The object must have, at least, the dimensions latitude, longitude and start date. The
+object is expected to be already subset for the desired region.}
+
+\item{target_grid}{a character vector indicating the target grid to be passed to CDO.
+It must be a grid recognised by CDO or a NetCDF file.}
+
+\item{int_method}{a character vector indicating the regridding method to be passed to CDORemap.
+Accepted methods are "con", "bil", "bic", "nn", "con2". If "nn" method is to be used, CDO_1.9.8
+or newer version is required.}
+
+\item{log_reg_method}{a character vector indicating the logistic regression method to be used.
+Accepted methods are "ens_mean", "ens_mean_sd", "sorted_members". "ens_mean" uses the ensemble
+mean anomalies as predictors in the logistic regression, "ens_mean_sd" uses the ensemble
+mean anomalies and the ensemble spread (computed as the standard deviation of all the members)
+as predictors in the logistic regression, and "sorted_members" considers all the members
+ordered decreasingly as predictors in the logistic regression. Default method is "ens_mean".}
+
+\item{probs_cat}{a numeric vector indicating the percentile thresholds separating the
+climatological distribution into different classes (categories). Default to c(1/3, 2/3). See
+\code{\link[easyVerification]{convert2prob}}.}
+
+\item{return_most_likely_cat}{if TRUE, the function returns the most likely category. If
+FALSE, the function returns the probabilities for each category. Default to FALSE.}
+
+\item{points}{a list of two elements containing the point latitudes and longitudes
+of the locations to downscale the model data. The list must contain the two elements
+named as indicated in the parameters 'lat_dim' and 'lon_dim'. If the downscaling is
+to a point location, only regular grids are allowed for exp and obs. Only needed if the
+downscaling is to a point location.}
+
+\item{method_point_interp}{a character vector indicating the interpolation method to interpolate
+model gridded data into the point locations. Accepted methods are "nearest", "bilinear", "9point",
+"invdist4nn", "NE", "NW", "SE", "SW". Only needed if the downscaling is to a point location.}
+
+\item{lat_dim}{a character vector indicating the latitude dimension name in the element 'data'
+in exp and obs. Default set to "lat".}
+
+\item{lon_dim}{a character vector indicating the longitude dimension name in the element 'data'
+in exp and obs. Default set to "lon".}
+
+\item{sdate_dim}{a character vector indicating the start date dimension name in the element
+'data' in exp and obs. Default set to "sdate".}
+
+\item{member_dim}{a character vector indicating the member dimension name in the element
+'data' in exp and obs. Default set to "member".}
+
+\item{region}{a numeric vector indicating the borders of the region defined in obs.
+It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers
+to the left border, while lonmax refers to the right border. latmin indicates the lower
+border, whereas latmax indicates the upper border. If set to NULL (default), the function
+takes the first and last elements of the latitudes and longitudes.}
+
+\item{loocv}{a logical vector indicating whether to perform leave-one-out cross-validation
+in the fitting of the logistic regression. Default to FALSE.}
+
+\item{ncores}{an integer indicating the number of cores to use in parallel computation.}
+
+\item{source_file}{a character vector with a path to an example file of the exp data.
+Only needed if the downscaling is to a point location.}
+}
+\value{
+An list of three elements. 'data' contains the dowscaled data, that could be either
+in the form of probabilities for each category or the most likely category. 'lat' contains the
+downscaled latitudes, and 'lon' the downscaled longitudes.
+}
+\description{
+This function performs a downscaling using an interpolation and a logistic
+regression. See \code{\link[nnet]{multinom}} for further details. It is recommended that
+the observations are passed already in the target grid. Otherwise, the function will also
+perform an interpolation of the observed field into the target grid. The coarse scale and
+observation data can be either global or regional. In the latter case, the region is
+defined by the user. In principle, the coarse and observation data are intended to be of
+the same variable, although different variables can also be admitted.
+}
+\examples{
+exp <- rnorm(1500) 
+dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 15) 
+exp_lons <- 1:5 
+exp_lats <- 1:4 
+obs <- rnorm(2700) 
+dim(obs) <- c(lat = 12, lon = 15, sdate = 15) 
+obs_lons <- seq(1,5, 4/14) 
+obs_lats <- seq(1,4, 3/11) 
+exp <- s2dv_cube(data = exp, lat = exp_lats, lon = exp_lons)
+obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons)
+res <- CST_LogisticReg(exp = exp, obs = obs, int_method = 'bil', target_grid = 'r1280x640', 
+probs_cat = c(1/3, 2/3))
+}
+\seealso{
+\code{\link[nnet]{multinom}}
+}
+\author{
+J. Ramon, \email{jaume.ramon@bsc.es}
+}
diff --git a/man/Intbc.Rd b/man/Intbc.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..55815d512070174156fa2338b28f7950fa31fc7f
--- /dev/null
+++ b/man/Intbc.Rd
@@ -0,0 +1,132 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Intbc.R
+\name{Intbc}
+\alias{Intbc}
+\title{Downscaling using interpolation and bias adjustment.}
+\usage{
+Intbc(
+  exp,
+  obs,
+  exp_lats,
+  exp_lons,
+  obs_lats,
+  obs_lons,
+  target_grid,
+  bc_method,
+  int_method = NULL,
+  points = NULL,
+  method_point_interp = NULL,
+  lat_dim = "lat",
+  lon_dim = "lon",
+  sdate_dim = "sdate",
+  time_dim = "time",
+  member_dim = "member",
+  source_file = NULL,
+  region = NULL,
+  ncores = 1,
+  ...
+)
+}
+\arguments{
+\item{exp}{an array with named dimensions containing the experimental field on the
+coarse scale for which the downscaling is aimed. The object must have, at least,
+the dimensions latitude, longitude, start date and member. The object is expected to be
+already subset for the desired region. Data can be in one or two integrated regions, e.g.,
+crossing the Greenwich meridian. To get the correct results in the latter case,
+the borders of the region should be specified in the parameter 'region'. See parameter
+'region'.}
+
+\item{obs}{an array with named dimensions containing the observational field. The object
+must have, at least, the dimensions latitude, longitude and start date. The object is
+expected to be already subset for the desired region.}
+
+\item{exp_lats}{a numeric vector containing the latitude values in 'exp'. Latitudes must
+range from -90 to 90.}
+
+\item{exp_lons}{a numeric vector containing the longitude values in 'exp'. Longitudes
+can range from -180 to 180 or from 0 to 360.}
+
+\item{obs_lats}{a numeric vector containing the latitude values in 'obs'. Latitudes must
+range from -90 to 90.}
+
+\item{obs_lons}{a numeric vector containing the longitude values in 'obs'. Longitudes
+can range from -180 to 180 or from 0 to 360.}
+
+\item{target_grid}{a character vector indicating the target grid to be passed to CDO.
+It must be a grid recognised by CDO or a NetCDF file.}
+
+\item{bc_method}{a character vector indicating the bias adjustment method to be applied after
+the interpolation. Accepted methods are 'simple_bias', 'calibration', 'quantile_mapping'. The
+abbreviations 'sbc', 'cal', 'qm' can also be used.}
+
+\item{int_method}{a character vector indicating the regridding method to be passed to CDORemap.
+Accepted methods are "con", "bil", "bic", "nn", "con2". If "nn" method is to be used, CDO_1.9.8
+or newer version is required.}
+
+\item{points}{a list of two elements containing the point latitudes and longitudes
+of the locations to downscale the model data. The list must contain the two elements
+named as indicated in the parameters 'lat_dim' and 'lon_dim'. If the downscaling is
+to a point location, only regular grids are allowed for exp and obs. Only needed if the
+downscaling is to a point location.}
+
+\item{method_point_interp}{a character vector indicating the interpolation method to interpolate
+model gridded data into the point locations. Accepted methods are "nearest", "bilinear", "9point",
+"invdist4nn", "NE", "NW", "SE", "SW". Only needed if the downscaling is to a point location.}
+
+\item{lat_dim}{a character vector indicating the latitude dimension name in the element 'data'
+in exp and obs. Default set to "lat".}
+
+\item{lon_dim}{a character vector indicating the longitude dimension name in the element 'data'
+in exp and obs. Default set to "lon".}
+
+\item{sdate_dim}{a character vector indicating the start date dimension name in the element
+'data' in exp and obs. Default set to "sdate".}
+
+\item{member_dim}{a character vector indicating the member dimension name in the element
+'data' in exp and obs. Default set to "member".}
+
+\item{source_file}{a character vector with a path to an example file of the exp data.
+Only needed if the downscaling is to a point location.}
+
+\item{region}{a numeric vector indicating the borders of the region defined in obs.
+It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers
+to the left border, while lonmax refers to the right border. latmin indicates the lower
+border, whereas latmax indicates the upper border. If set to NULL (default), the function
+takes the first and last elements of the latitudes and longitudes.}
+
+\item{ncores}{an integer indicating the number of cores to use in parallel computation.}
+}
+\value{
+An list of three elements. 'data' contains the dowscaled field, 'lat' the
+downscaled latitudes, and 'lon' the downscaled longitudes.
+}
+\description{
+This function performs a downscaling using an interpolation and a later bias
+adjustment. It is recommended that the observations are passed already in the target grid.
+Otherwise, the function will also perform an interpolation of the observed field into the
+target grid. The coarse scale and observation data can be either global or regional. In the
+latter case, the region is defined by the user. In principle, the coarse and observation data
+are intended to be of the same variable, although different variables can also be admitted.
+}
+\examples{
+exp <- rnorm(500) 
+dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 5) 
+exp_lons <- 1:5 
+exp_lats <- 1:4 
+obs <- rnorm(900) 
+dim(obs) <- c(lat = 12, lon = 15, sdate = 5) 
+obs_lons <- seq(1,5, 4/14) 
+obs_lats <- seq(1,4, 3/11) 
+res <- Intbc(exp = exp, obs = obs, exp_lats = exp_lats, exp_lons = exp_lons, obs_lats = obs_lats, 
+obs_lons = obs_lons, target_grid = 'r1280x640', bc_method = 'simple_bias', int_method = 'conservative')
+}
+\seealso{
+\code{\link[CSTools]{BiasCorrection}}
+
+\code{\link[CSTools]{Calibration}}
+
+\code{\link[CSTools]{QuantileMapping}}
+}
+\author{
+J. Ramon, \email{jaume.ramon@bsc.es}
+}
diff --git a/man/Interpolation.Rd b/man/Interpolation.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..c33b1b642580aa3e96af5e1c41bdb3156f04d922
--- /dev/null
+++ b/man/Interpolation.Rd
@@ -0,0 +1,94 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Interpolation.R
+\name{Interpolation}
+\alias{Interpolation}
+\title{Regrid or interpolate gridded data to a point location.}
+\usage{
+Interpolation(
+  exp,
+  lats,
+  lons,
+  points = NULL,
+  source_file = NULL,
+  method_remap = NULL,
+  target_grid = NULL,
+  lat_dim = "lat",
+  lon_dim = "lon",
+  region = NULL,
+  method_point_interp = NULL
+)
+}
+\arguments{
+\item{exp}{an array with named dimensions containing the experimental field on the
+coarse scale for which the downscaling is aimed. The object must have, at least,
+the dimensions latitude and longitude. The object is expected to be already subset
+for the desired region. Data can be in one or two integrated regions, e.g.,
+crossing the Greenwich meridian. To get the correct results in the latter case,
+the borders of the region should be specified in the parameter 'region'. See parameter
+'region'.}
+
+\item{lats}{a numeric vector containing the latitude values. Latitudes must range from
+-90 to 90.}
+
+\item{lons}{a numeric vector containing the longitude values. Longitudes can range from
+-180 to 180 or from 0 to 360.}
+
+\item{points}{a list of two elements containing the point latitudes and longitudes
+of the locations to downscale the model data. The list must contain the two elements
+named as indicated in the parameters 'lat_dim' and 'lon_dim'. If the downscaling is
+to a point location, only regular grids are allowed for exp and obs. Only needed if the
+downscaling is to a point location.}
+
+\item{source_file}{a character vector with a path to an example file of the exp data.
+Only needed if the downscaling is to a point location.}
+
+\item{method_remap}{a character vector indicating the regridding method to be passed
+to CDORemap. Accepted methods are "con", "bil", "bic", "nn", "con2". If "nn" method is
+to be used, CDO_1.9.8 or newer version is required.}
+
+\item{target_grid}{a character vector indicating the target grid to be passed to CDO.
+It must be a grid recognised by CDO or a NetCDF file.}
+
+\item{lat_dim}{a character vector indicating the latitude dimension name in the element
+'exp' and/or 'points'. Default set to "lat".}
+
+\item{lon_dim}{a character vector indicating the longitude dimension name in the element
+'exp' and/or 'points'. Default set to "lon".}
+
+\item{region}{a numeric vector indicating the borders of the region defined in exp.
+It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers
+to the left border, while lonmax refers to the right border. latmin indicates the lower
+border, whereas latmax indicates the upper border. If set to NULL (default), the function
+takes the first and last elements of the latitudes and longitudes.}
+
+\item{method_point_interp}{a character vector indicating the interpolation method to
+interpolate model gridded data into the point locations. Accepted methods are "nearest",
+"bilinear", "9point", "invdist4nn", "NE", "NW", "SE", "SW". Only needed if the downscaling
+is to a point location.}
+}
+\value{
+An list of three elements. 'data' contains the dowscaled field, 'lat' the
+downscaled latitudes, and 'lon' the downscaled longitudes.
+}
+\description{
+This function interpolates gridded model data from one grid to
+another (regrid) or interpolates gridded model data to a set of point locations.
+The gridded model data can be either global or regional. In the latter case, the
+region is defined by the user. It does not have constrains of specific region or
+variables to downscale.
+}
+\examples{
+exp <- rnorm(500) 
+dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 5, time  = 1) 
+lons <- 1:5 
+lats <- 1:4 
+res <- Interpolation(exp = exp, lats = lats, lons = lons, method_remap = 'conservative', target_grid = 'r1280x640')
+}
+\seealso{
+\code{\link[s2dverification]{CDORemap}}
+}
+\author{
+J. Ramon, \email{jaume.ramon@bsc.es}
+
+Ll. Lledó, \email{llorenc.lledo@ecmwf.int}
+}
diff --git a/man/Intlr.Rd b/man/Intlr.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..305ad0846d37ddd0097f701726fc300b27e01091
--- /dev/null
+++ b/man/Intlr.Rd
@@ -0,0 +1,152 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Intlr.R
+\name{Intlr}
+\alias{Intlr}
+\title{Downscaling using interpolation and linear regression.}
+\usage{
+Intlr(
+  exp,
+  obs,
+  exp_lats,
+  exp_lons,
+  obs_lats,
+  obs_lons,
+  lr_method,
+  target_grid = NULL,
+  points = NULL,
+  int_method = NULL,
+  method_point_interp = NULL,
+  source_file_exp = NULL,
+  source_file_obs = NULL,
+  predictors = NULL,
+  lat_dim = "lat",
+  lon_dim = "lon",
+  sdate_dim = "sdate",
+  time_dim = "time",
+  member_dim = "member",
+  region = NULL,
+  large_scale_predictor_dimname = "vars",
+  loocv = FALSE,
+  ncores = 1
+)
+}
+\arguments{
+\item{exp}{an array with named dimensions containing the experimental field on the
+coarse scale for which the downscaling is aimed. The object must have, at least,
+the dimensions latitude, longitude and start date. The object is expected to be
+already subset for the desired region. Data can be in one or two integrated regions, e.g.,
+crossing the Greenwich meridian. To get the correct results in the latter case,
+the borders of the region should be specified in the parameter 'region'. See parameter
+'region'.}
+
+\item{obs}{an array with named dimensions containing the observational field. The object
+must have, at least, the dimensions latitude, longitude and start date. The object is
+expected to be already subset for the desired region.}
+
+\item{exp_lats}{a numeric vector containing the latitude values in 'exp'. Latitudes must
+range from -90 to 90.}
+
+\item{exp_lons}{a numeric vector containing the longitude values in 'exp'. Longitudes
+can range from -180 to 180 or from 0 to 360.}
+
+\item{obs_lats}{a numeric vector containing the latitude values in 'obs'. Latitudes must
+range from -90 to 90.}
+
+\item{obs_lons}{a numeric vector containing the longitude values in 'obs'. Longitudes
+can range from -180 to 180 or from 0 to 360.}
+
+\item{lr_method}{a character vector indicating the linear regression method to be applied
+after the interpolation. Accepted methods are 'basic', 'large-scale' and '4nn'. The 'basic'
+method fits a linear regression using high resolution observations as predictands and the
+interpolated model data as predictor. Then, the regression equation is to the interpolated
+model data to correct the interpolated values. The 'large-scale' method fits a linear
+regression with large-scale predictors from the same model (e.g. teleconnection indices)
+as predictors and high-resolution observations as predictands. This equation is then
+applied to the interpolated model values. Finally, the '4nn' method uses a linear
+regression with the four nearest neighbours as predictors and high-resolution observations
+as predictands. It is then applied to model data to correct the interpolated values.}
+
+\item{target_grid}{a character vector indicating the target grid to be passed to CDO.
+It must be a grid recognised by CDO or a NetCDF file.}
+
+\item{points}{a list of two elements containing the point latitudes and longitudes
+of the locations to downscale the model data. The list must contain the two elements
+named as indicated in the parameters 'lat_dim' and 'lon_dim'. If the downscaling is
+to a point location, only regular grids are allowed for exp and obs. Only needed if the
+downscaling is to a point location.}
+
+\item{int_method}{a character vector indicating the regridding method to be passed
+to CDORemap. Accepted methods are "con", "bil", "bic", "nn", "con2". If "nn" method is
+to be used, CDO_1.9.8 or newer version is required.}
+
+\item{method_point_interp}{a character vector indicating the interpolation method to
+interpolate model gridded data into the point locations. Accepted methods are "nearest",
+"bilinear", "9point", "invdist4nn", "NE", "NW", "SE", "SW".}
+
+\item{source_file_exp}{a character vector with a path to an example file of the exp data.
+Only needed if the downscaling is to a point location.}
+
+\item{source_file_obs}{a character vector with a path to an example file of the obs data.
+Only needed if the downscaling is to a point location.}
+
+\item{predictors}{an array with large-scale data to be used in the 'large-scale' method.
+Only needed if the linear regression method is set to 'large-scale'. It must have, at
+least the dimension start date and another dimension whose name has to be specified in
+the parameter 'large_scale_predictor_dimname'. It should contain as many elements as the
+number of large-scale predictors.}
+
+\item{lat_dim}{a character vector indicating the latitude dimension name in the element 'data'
+in exp and obs. Default set to "lat".}
+
+\item{lon_dim}{a character vector indicating the longitude dimension name in the element 'data'
+in exp and obs. Default set to "lon".}
+
+\item{sdate_dim}{a character vector indicating the start date dimension name in the element
+'data' in exp and obs. Default set to "sdate".}
+
+\item{time_dim}{a character vector indicating the time dimension name in the element
+'data' in exp and obs. Default set to "time".}
+
+\item{region}{a numeric vector indicating the borders of the region defined in exp.
+It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers
+to the left border, while lonmax refers to the right border. latmin indicates the lower
+border, whereas latmax indicates the upper border. If set to NULL (default), the function
+takes the first and last elements of the latitudes and longitudes.}
+
+\item{large_scale_predictor_dimname}{a character vector indicating the name of the
+dimension in 'predictors' that contain the predictor variables. See parameter 'predictors'.}
+
+\item{loocv}{a logical indicating whether to apply leave-one-out cross-validation when
+generating the linear regressions. Default to FALSE.}
+
+\item{ncores}{an integer indicating the number of cores to use in parallel computation.}
+}
+\value{
+A list of three elements. 'data' contains the dowscaled field, 'lat' the
+downscaled latitudes, and 'lon' the downscaled longitudes.
+}
+\description{
+This function performs a downscaling using an interpolation and a linear
+regression. Different methodologies that employ linear regressions are available. See
+parameter 'lr_method' for more information. It is recommended that the observations
+are passed already in the target grid. Otherwise, the function will also perform an
+interpolation of the observed field into the target grid. The coarse scale and
+observation data can be either global or regional. In the latter case, the region is
+defined by the user. In principle, the coarse and observation data are intended to
+be of the same variable, although different variables can also be admitted.
+}
+\examples{
+exp <- rnorm(500) 
+dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 5) 
+exp_lons <- 1:5 
+exp_lats <- 1:4 
+obs <- rnorm(900) 
+dim(obs) <- c(lat = 12, lon = 15, sdate = 5) 
+obs_lons <- seq(1,5, 4/14) 
+obs_lats <- seq(1,4, 3/11) 
+res <- Intlr(exp = exp, obs = obs, exp_lats = exp_lats, exp_lons = exp_lons, obs_lats = obs_lats, 
+obs_lons = obs_lons, target_grid = 'r1280x640', lr_method = 'basic', int_method = 'conservative')
+}
+\author{
+J. Ramon, \email{jaume.ramon@bsc.es}
+}
diff --git a/man/LogisticReg.Rd b/man/LogisticReg.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..55146faca8d576bb3b1ccb0db5e85801976bf230
--- /dev/null
+++ b/man/LogisticReg.Rd
@@ -0,0 +1,145 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/LogisticReg.R
+\name{LogisticReg}
+\alias{LogisticReg}
+\title{Downscaling using interpolation and logistic regression.}
+\usage{
+LogisticReg(
+  exp,
+  obs,
+  exp_lats,
+  exp_lons,
+  obs_lats,
+  obs_lons,
+  target_grid,
+  int_method = NULL,
+  log_reg_method = "ens_mean",
+  probs_cat = c(1/3, 2/3),
+  return_most_likely_cat = FALSE,
+  points = NULL,
+  method_point_interp = NULL,
+  lat_dim = "lat",
+  lon_dim = "lon",
+  sdate_dim = "sdate",
+  member_dim = "member",
+  source_file = NULL,
+  region = NULL,
+  loocv = FALSE,
+  ncores = 1
+)
+}
+\arguments{
+\item{exp}{an array with named dimensions containing the experimental field on the
+coarse scale for which the downscaling is aimed. The object must have, at least,
+the dimensions latitude, longitude, start date and member. The object is expected to be
+already subset for the desired region. Data can be in one or two integrated regions, e.g.,
+crossing the Greenwich meridian. To get the correct results in the latter case,
+the borders of the region should be specified in the parameter 'region'. See parameter
+'region'.}
+
+\item{obs}{an array with named dimensions containing the observational field. The object
+must have, at least, the dimensions latitude, longitude and start date. The object is
+expected to be already subset for the desired region.}
+
+\item{exp_lats}{a numeric vector containing the latitude values in 'exp'. Latitudes must
+range from -90 to 90.}
+
+\item{exp_lons}{a numeric vector containing the longitude values in 'exp'. Longitudes
+can range from -180 to 180 or from 0 to 360.}
+
+\item{obs_lats}{a numeric vector containing the latitude values in 'obs'. Latitudes must
+range from -90 to 90.}
+
+\item{obs_lons}{a numeric vector containing the longitude values in 'obs'. Longitudes
+can range from -180 to 180 or from 0 to 360.}
+
+\item{target_grid}{a character vector indicating the target grid to be passed to CDO.
+It must be a grid recognised by CDO or a NetCDF file.}
+
+\item{int_method}{a character vector indicating the regridding method to be passed to CDORemap.
+Accepted methods are "con", "bil", "bic", "nn", "con2". If "nn" method is to be used, CDO_1.9.8
+or newer version is required.}
+
+\item{log_reg_method}{a character vector indicating the logistic regression method to be used.
+Accepted methods are "ens_mean", "ens_mean_sd", "sorted_members". "ens_mean" uses the ensemble
+mean anomalies as predictors in the logistic regression, "ens_mean_sd" uses the ensemble
+mean anomalies and the ensemble spread (computed as the standard deviation of all the members)
+as predictors in the logistic regression, and "sorted_members" considers all the members
+ordered decreasingly as predictors in the logistic regression. Default method is "ens_mean".}
+
+\item{probs_cat}{a numeric vector indicating the percentile thresholds separating the
+climatological distribution into different classes (categories). Default to c(1/3, 2/3). See
+\code{\link[easyVerification]{convert2prob}}.}
+
+\item{return_most_likely_cat}{if TRUE, the function returns the most likely category. If
+FALSE, the function returns the probabilities for each category. Default to FALSE.}
+
+\item{points}{a list of two elements containing the point latitudes and longitudes
+of the locations to downscale the model data. The list must contain the two elements
+named as indicated in the parameters 'lat_dim' and 'lon_dim'. If the downscaling is
+to a point location, only regular grids are allowed for exp and obs. Only needed if the
+downscaling is to a point location.}
+
+\item{method_point_interp}{a character vector indicating the interpolation method to interpolate
+model gridded data into the point locations. Accepted methods are "nearest", "bilinear", "9point",
+"invdist4nn", "NE", "NW", "SE", "SW". Only needed if the downscaling is to a point location.}
+
+\item{lat_dim}{a character vector indicating the latitude dimension name in the element 'data'
+in exp and obs. Default set to "lat".}
+
+\item{lon_dim}{a character vector indicating the longitude dimension name in the element 'data'
+in exp and obs. Default set to "lon".}
+
+\item{sdate_dim}{a character vector indicating the start date dimension name in the element
+'data' in exp and obs. Default set to "sdate".}
+
+\item{member_dim}{a character vector indicating the member dimension name in the element
+'data' in exp and obs. Default set to "member".}
+
+\item{source_file}{a character vector with a path to an example file of the exp data.
+Only needed if the downscaling is to a point location.}
+
+\item{region}{a numeric vector indicating the borders of the region defined in obs.
+It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers
+to the left border, while lonmax refers to the right border. latmin indicates the lower
+border, whereas latmax indicates the upper border. If set to NULL (default), the function
+takes the first and last elements of the latitudes and longitudes.}
+
+\item{loocv}{a logical vector indicating whether to perform leave-one-out cross-validation
+in the fitting of the logistic regression. Default to FALSE.}
+
+\item{ncores}{an integer indicating the number of cores to use in parallel computation.}
+}
+\value{
+An list of three elements. 'data' contains the dowscaled data, that could be either
+in the form of probabilities for each category or the most likely category. 'lat' contains the
+downscaled latitudes, and 'lon' the downscaled longitudes.
+}
+\description{
+This function performs a downscaling using an interpolation and a logistic
+regression. See \code{\link[nnet]{multinom}} for further details. It is recommended that
+the observations are passed already in the target grid. Otherwise, the function will also
+perform an interpolation of the observed field into the target grid. The coarse scale and
+observation data can be either global or regional. In the latter case, the region is
+defined by the user. In principle, the coarse and observation data are intended to be of
+the same variable, although different variables can also be admitted.
+}
+\examples{
+exp <- rnorm(1500) 
+dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 15) 
+exp_lons <- 1:5 
+exp_lats <- 1:4 
+obs <- rnorm(2700) 
+dim(obs) <- c(lat = 12, lon = 15, sdate = 15) 
+obs_lons <- seq(1,5, 4/14) 
+obs_lats <- seq(1,4, 3/11) 
+res <- LogisticReg(exp = exp, obs = obs, exp_lats = exp_lats, exp_lons = exp_lons, 
+obs_lats = obs_lats, obs_lons = obs_lons, int_method = 'bil', target_grid = 'r1280x640', 
+probs_cat = c(1/3, 2/3))
+}
+\seealso{
+\code{\link[nnet]{multinom}}
+}
+\author{
+J. Ramon, \email{jaume.ramon@bsc.es}
+}
diff --git a/tests/testthat/test_CST_Intbc.R b/tests/testthat/test_CST_Intbc.R
index 122c50b12c8e7eeca730907793399c5ba8927c06..61d9270cdf0eca44c845cbcec9208e880ab1d0ff 100644
--- a/tests/testthat/test_CST_Intbc.R
+++ b/tests/testthat/test_CST_Intbc.R
@@ -19,25 +19,24 @@ test_that("Sanity checks", {
   exp_lats <- 1:4
   exp <- suppressWarnings(s2dv_cube(data = exp, lat = exp_lats, lon = exp_lons)) 
 
-  expect_error(CST_Intbc(exp = exp), "argument \"obs\" is missing, with no default")
-
   obs <- rnorm(900)
   dim(obs) <- c(lat = 12, lon = 15, sdate = 5, time = 1)
   obs_lons <- seq(1,5, 4/14)
   obs_lats <- seq(1,4, 3/11)
-  obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons)
+  obs <- suppressWarnings(s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons))
 
   expect_error(
-    CST_Intbc(exp = exp, obs = obs, bc_method = 'cal'),
+    supressWarnings(CST_Intbc(exp = exp, obs = obs, bc_method = 'cal'),
     paste0("Parameter 'method_remap' must be a character vector indicating the ",
-      "interpolation method. Accepted methods are con, bil, bic, nn, con2"))
+      "interpolation method. Accepted methods are con, bil, bic, nn, con2")))
 
   obs2 <- rnorm(180)
   dim(obs2) <- c(lat = 12, lon = 15, sdate = 1)
-  obs2 <- s2dv_cube(data = obs2, lat = obs_lats, lon = obs_lons)  
+  obs2 <- suppressWarnings(s2dv_cube(data = obs2, lat = obs_lats, lon = obs_lons))  
 
   expect_warning(
-    CST_Intbc(exp = exp, obs = obs2, int_method = 'bil', bc_method = 'sbc', target_grid = 'r1280x640'),
+    CST_Intbc(exp = exp, obs = obs2, int_method = 'bil', bc_method = 'sbc', 
+              target_grid = 'r1280x640'),
     'Simple Bias Correction should not be used with only one observation. Returning NA.')
  
   exp2 <- rnorm(100)
@@ -45,17 +44,21 @@ test_that("Sanity checks", {
   exp2 <- suppressWarnings(s2dv_cube(data = exp2, lat = exp_lats, lon = exp_lons))
 
   expect_error( 
-    CST_Intbc(exp = exp2, obs = obs, int_method = 'bil', bc_method = 'cal', target_grid = 'r1280x640'),	       
-    'Calibration must not be used with only one ensemble member.')	 
+    suppressWarnings(CST_Intbc(exp = exp2, obs = obs, int_method = 'bil', bc_method = 'cal', 
+                              target_grid = 'r1280x640'), 'Calibration must not be used 
+                              with only one ensemble member.'))	 
 
-  d1 <- CST_Intbc(exp = exp, obs = obs, int_method = 'bil', bc_method = 'sbc', target_grid = 'r640x320')
-  expect_equal(round(d1$data[, 1, 1, 1, 1], 2), c(-0.40, 0.08, -0.23, 0.38, -0.03))
+  d1 <- suppressWarnings(CST_Intbc(exp = exp, obs = obs, int_method = 'bil', bc_method = 'sbc', 
+                                   target_grid = 'r640x320'))
+  expect_equal(round(d1$exp$data[, 1, 1, 1, 1], 2), c(-0.40, 0.08, -0.23, 0.38, -0.03))
 
-  d2 <- CST_Intbc(exp = exp, obs = obs, int_method = 'bil', bc_method = 'cal', target_grid = 'r640x320')
-  expect_equal(round(d2$data[, 1, 1, 1, 1], 2), c(-0.47, -0.03, -0.31, 0.24, -0.13))
+  d2 <- suppressWarnings(CST_Intbc(exp = exp, obs = obs, int_method = 'bil', bc_method = 'cal', 
+                                   target_grid = 'r640x320'))
+  expect_equal(round(d2$exp$data[, 1, 1, 1, 1], 2), c(-0.47, -0.03, -0.31, 0.24, -0.13))
 
-  d3 <- CST_Intbc(exp = exp, obs = obs, int_method = 'bil', bc_method = 'qm', target_grid = 'r640x320')
-  expect_equal(round(d3$data[, 1, 1, 1, 1], 2), c(0.00, 0.00, 0.00, 0.33, 0.00))
+  d3 <- suppressWarnings(CST_Intbc(exp = exp, obs = obs, int_method = 'bil', bc_method = 'qm', 
+                        target_grid = 'r640x320'))
+  expect_equal(round(d3$exp$data[, 1, 1, 1, 1], 2), c(0.00, 0.00, 0.00, 0.33, 0.00))
   })
     
 
diff --git a/tests/testthat/test_CST_Interpolation.R b/tests/testthat/test_CST_Interpolation.R
index ae9dcc5a9a887a92ffebc02b56136d23f25b5bc3..5d8f3d292d6e28316cd823c87cc51078c6018217 100644
--- a/tests/testthat/test_CST_Interpolation.R
+++ b/tests/testthat/test_CST_Interpolation.R
@@ -1,6 +1,7 @@
 source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Interpolation.R')
 require(testthat)
 require(CSTools)
+require(s2dv)
 
 context("Generic tests")
 test_that("Sanity checks", {
@@ -39,10 +40,12 @@ test_that("Sanity checks", {
 
   d <- CST_Interpolation(exp = exp, method_remap = 'bil', target_grid = 'r1280x640')
 
-  expect_equal(length(d), 8)
-  expect_equal(class(d), "s2dv_cube")
-  expect_equal(as.numeric(dim(d$data)[names(dim(d$data)) == 'lat']), length(d$lat))
-  expect_equal(as.numeric(dim(d$data)[names(dim(d$data)) == 'lon']), length(d$lon))
+  expect_equal(length(d$exp), 9)
+  expect_equal(d$obs, NULL)
+  expect_equal(class(d), "list")
+  expect_equal(class(d$exp), "s2dv_cube")
+  expect_equal(as.numeric(dim(d$exp$data)[names(dim(d$exp$data)) == 'lat']), length(d$exp$lat))
+  expect_equal(as.numeric(dim(d$exp$data)[names(dim(d$exp$data)) == 'lon']), length(d$exp$lon))
 
   expect_error(
     CST_Interpolation(exp = exp, points = 1),
@@ -76,7 +79,7 @@ test_that("Sanity checks", {
   lats <- 1:3
   exp <- suppressWarnings(s2dv_cube(data = exp, lat = lats, lon = lons))
 
-  d2 <- CST_Interpolation(exp = exp, method_remap = 'bil', target_grid = 'r640x320')$data
-  expect_equal(d2, array(c(2, 2, 3, 3, 4, 5), dim = c(lat = 3, lon = 2)))
+  d2 <- CST_Interpolation(exp = exp, method_remap = 'bil', target_grid = 'r640x320')
+  expect_equal(d2$exp$data, array(c(2, 2, 3, 3, 4, 5), dim = c(lat = 3, lon = 2)))
   })