diff --git a/DESCRIPTION b/DESCRIPTION
index 0492186bda3b08f64f7554927f06942a8ae43bbc..09b9cf06450e0e4b680739590b7975f895565646 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -47,6 +47,7 @@ Depends:
     maps
 Imports:
     s2dverification,
+    s2dv,
     rainfarmr,
     multiApply (>= 2.1.1),
     qmap,
diff --git a/NAMESPACE b/NAMESPACE
index 342847f35714e9473b71200e7c58f7f19515e5ac..e070ff78e4d57862b7f3fa934184e1c1d9d2ffea 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -37,6 +37,7 @@ export(RFSlope)
 export(RFTemp)
 export(RainFARM)
 export(RegimesAssign)
+export(SaveExp)
 export(SplitDim)
 export(WeatherRegime)
 export(as.s2dv_cube)
@@ -77,6 +78,7 @@ importFrom(maps,map)
 importFrom(plyr,.)
 importFrom(plyr,dlply)
 importFrom(reshape2,melt)
+importFrom(s2dv,Reorder)
 importFrom(utils,glob2rx)
 importFrom(utils,head)
 importFrom(utils,tail)
diff --git a/NEWS.md b/NEWS.md
index 19caa068dc1f95e7f3a710ae5d05e45292f78317..1a9069e5925f3a2869dc8c18f823b2a5e0a5baa7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -20,6 +20,7 @@
     + DESCRIPTION specifies the minimum multiApply version required
     + EnsClustering has a fixed 'closest_member' output
     + PlotCombinedMap handles masks correctly
+    + CST_SaveExp uses multiApply and save time dimension correctly
 
 ### CSTools 3.0.0
 **Submission date to CRAN: 10-02-2020**
diff --git a/R/CST_SaveExp.R b/R/CST_SaveExp.R
index 1b3b4b1011179ec8a529c3914b40d64287f651a5..c3e04c2a2b2d474fad708649a24775f3f5c8af44 100644
--- a/R/CST_SaveExp.R
+++ b/R/CST_SaveExp.R
@@ -16,9 +16,10 @@
 #'
 #'@seealso \code{\link{CST_Load}}, \code{\link{as.s2dv_cube}} and \code{\link{s2dv_cube}}
 #'
-#'@import s2dverification
+#'@importFrom ClimProjDiags Subset
 #'@import ncdf4
-#'@import abind
+#'@importFrom s2dv Reorder
+#'@import multiApply
 #'
 #'@examples
 #'\dontrun{
@@ -39,33 +40,97 @@ CST_SaveExp <- function(data, destination = "./CST_Data") {
     stop("Parameter 'data' must be of the class 's2dv_cube', ",
          "as output by CSTools::CST_Load.")
   }
-dimname <- names(dim(data$data))
-  if (any(dimname == "time")) {
-    dimname[which(dimname == "time")] <- "ftime"
-    names(dim(data$data))[which(dimname == "time")] <- "ftime"
+  sdates <- lapply(1:length(data$Datasets), function(x) {
+                          data$Datasets[[x]]$InitializationDates[[1]]})[[1]]
+  if (!is.character(attributes(data$Variable)$units)) {
+      units <- attributes(data$Variable)$variable$units
+  } else {
+      units <- attributes(data$Variable)$units
+  }
+  cdo_grid_name = attr(data$lon, 'cdo_grid_name')
+  projection = attr(data$lon, 'projection')
+  var_name <- data$Variable$varName
+  time_values <- data$Dates$start
+  dim(time_values) <- c(time = length(time_values) / length(sdates), 
+                        sdate = length(sdates)) 
+  SaveExp(data = data$data, lon = data$lon, lat = data$lat,
+          Dataset = names(data$Datasets), var_name = var_name,
+          units = units, cdo_grid_name = cdo_grid_name, projection = projection,
+          startdates = sdates, Dates = time_values, destination) 
+} 
+#'Save an experiment in a format compatible with CST_Load
+#'@description This function is created for compatibility with CST_Load/Load for saving post-processed datasets such as those calibrated of downscaled with CSTools functions
+#'
+#'@author Perez-Zanon Nuria, \email{nuria.perez@bsc.es}
+#'
+#'@param data an multi-dimensional array with named dimensions (longitude, latitude, time, member, sdate)
+#'@param lon vector of logitud corresponding to the longitudinal dimension in data
+#'@param lat vector of latitud corresponding to the latitudinal dimension in data
+#'@param Dataset a vector of character string indicating the names of the datasets
+#'@param var_name a character string indicating the name of the variable to be saved
+#'@param units a character string indicating the units of the variable
+#'@param startdates a vector of dates indicating the initialization date of each simulations
+#'@param Dates a matrix of dates with two dimension 'time' and 'sdate'.
+#'@param cdo_grid_name a character string indicating the name of the grid e.g.: 'r360x181'
+#'@param projection a character string indicating the projection name
+#'@param destination a character string indicating the path where to store the NetCDF files
+#'
+#'@return the function creates as many files as sdates per dataset. Each file could contain multiple members
+#' The path will be created with the name of the variable and each Datasets.
+#' 
+#'@importFrom ClimProjDiags Subset
+#'@import ncdf4
+#'@importFrom s2dv Reorder
+#'@import multiApply
+#'
+#'@examples
+#'\dontrun{
+#'data <- lonlat_data$exp$data
+#'lon <- lonlat_data$exp$lon
+#'lat <- lonlat_data$exp$lat
+#'Dataset <- 'XXX'
+#'var_name <- 'tas'
+#'units <- 'k'
+#'startdates <- lapply(1:length(lonlat_data$exp$Datasets),
+#'                     function(x) {
+#'                         lonlat_data$exp$Datasets[[x]]$InitializationDates[[1]]})[[1]]
+#'Dates <- lonlat_data$exp$Dates$start
+#'dim(Dates) <- c(time = length(Dates)/length(startdates), sdate = length(startdates))
+#'cdo_grid_name = attr(lonlat_data$exp$lon, 'cdo_grid_name')
+#'projection = attr(lonlat_data$exp$lon, 'projection')
+#'destination = './path/'
+#'SaveExp(data, lon, lat, Dataset, var_name, units, startdates, Dates,
+#'                    cdo_grid_name, projection, destination)
+#'}
+#'@export
+SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates,
+                    cdo_grid_name, projection, destination) {  
+  dimname <- names(dim(data))
+  if (any(dimname == "ftime")) {
+    dimname[which(dimname == "ftime")] <- "time"
+    names(dim(data))[which(dimname == "ftime")] <- "time"
   }
   if (any(dimname == "memb")) {
     dimname[which(dimname == "memb")] <- "member"
-    names(dim(data$data))[which(dimname == "memb")] <- "member"
+    names(dim(data))[which(dimname == "memb")] <- "member"
   }
   if (any(dimname == "ensemble")) {
     dimname[which(dimname == "ensemble")] <- "member"
-    names(dim(data$data))[which(dimname == "ensemble")] <- "member"
+    names(dim(data))[which(dimname == "ensemble")] <- "member"
   }
-  if (any(dimname == "longitude")) {
-    dimname[which(dimname == "longitude")] <- "lon"
-    names(dim(data$data))[which(dimname == "longitude")] <- "lon"
+  if (any(dimname == "lon")) {
+    dimname[which(dimname == "lon")] <- "longitude"
+    names(dim(data))[which(dimname == "lon")] <- "longitude"
   }
-  if (any(dimname == "latitude")) {
-    dimname[which(dimname == "latitude")] <- "lat"
-    names(dim(data$data))[which(dimname == "latitude")] <- "lat"
+  if (any(dimname == "lat")) {
+    dimname[which(dimname == "lat")] <- "latitude"
+    names(dim(data))[which(dimname == "lat")] <- "latitude"
   }
-names(dim(data$data)) <- dimname
-
+  names(dim(data)) <- dimname
   if (is.null(dimname)) {
     stop("Element 'data' in parameter 'data' must have named dimensions.")
   }
-sdate_pos <- which(dimname == "sdate")
+  sdate_pos <- which(dimname == "sdate")
 
   if (length(sdate_pos) == 0) {
     stop("Element 'data' in parameter 'data' hasn't 'sdate' dimension.")
@@ -73,24 +138,23 @@ sdate_pos <- which(dimname == "sdate")
     stop("Element 'data' in parameter 'data' has more than one 'sdate'",
          " dimension.")
   }
-n_sdates <- dim(data$data)[sdate_pos] # number of files to create
 
-dataset_pos <- which(dimname == "dataset" | dimname == "dat")
-dims <- dim(data$data)
+  dataset_pos <- which(dimname == "dataset" | dimname == "dat")
+  dims <- dim(data)
   if (length(dataset_pos) == 0) {
     warning("Element 'data' in parameter 'data' hasn't 'dataset' dimension. ",
             "All data is stored in the same 'dataset' folder.")
-    data$data <- InsertDim(var = data$data, posdim = 1, lendim = 1)
-    names(dim(data$data))[1] <- "dataset"
+    data$data <- InsertDim(var = data, posdim = 1, lendim = 1)
+    names(dim(data))[1] <- "dataset"
     dimname <- c("dataset", dimname)
     dataset_pos = 1
   } else if (length(dataset_pos) > 1) {
     stop("Element 'data' in parameter 'data' has more than one 'dataset'",
          " dimension.")
   }
-n_datasets <- dim(data$data)[dataset_pos] # number of folder by dataset
-# dataset names:
-datasets <- names(data$Datasets)
+  n_datasets <- dim(data)[dataset_pos] # number of folder by dataset
+  # dataset names:
+  datasets <- Dataset
   if (n_datasets > length(datasets)) {
     warning("Dimension 'dataset' in element 'data' from parameter 'data' ",
             "is greater than those listed in element 'Datasets' and the ",
@@ -102,15 +166,14 @@ datasets <- names(data$Datasets)
             " first element will be used.")
     datasets <- datasets[1 : n_datasets]
   }
-# var names:
+  # var names:
   if ('var' %in% dimname) {
     var_pos <- which(dimname == 'var')
       if (dims[var_pos] == 1) {
-          data$data <- adrop(data$data, drop = var_pos)
-          dimname <- names(dim(data$data))
+          data <- adrop(data, drop = var_pos)
+          dimname <- names(dim(data))
       }
   }
-var_name <- data$Variable$varName
   if (length(var_name) != 1) {
     stop("One variable name must be included in element 'Variable$varName' ",
          "of parameter 'data'.") 
@@ -120,20 +183,20 @@ var_name <- data$Variable$varName
          "must be a character string.")
   }
 
-known_dim_names <- c("var", "lat", "latitude", "lon", "longitude", "time", 
-                     "ftime", "sdate", "dataset", "dat", "nlevel", "levels")
-dims_var <- NULL
-list_pos <- 1
+  known_dim_names <- c("var", "lat", "latitude", "lon", "longitude", "time", 
+                       "ftime", "sdate", "dataset", "dat", "nlevel", "levels")
+  dims_var <- NULL
+  list_pos <- 1
 
   if (any(dimname == 'longitude') | any(dimname == 'lon')) {
     dim_lon <- ncdim_def(name = 'lon', units = 'degrees',
-                         vals = as.vector(data$lon), longname = 'longitude')
+                         vals = as.vector(lon), longname = 'longitude')
     dims_var[[list_pos]] <- dim_lon
     list_pos <- list_pos + 1
   }
     if (any(dimname == 'latitude') | any(dimname == 'lat')) {
     dim_lat <- ncdim_def(name = 'lat', units = 'degrees_north',
-                         vals = as.vector(data$lat), longname = 'latitude')
+                         vals = as.vector(lat), longname = 'latitude')
     dims_var[[list_pos]] <- dim_lat
     list_pos <- list_pos + 1
   }
@@ -143,15 +206,13 @@ list_pos <- 1
       stop("Ask for saving realizations or further dimensions to the mantainer.")
     } else {
       dim_memb <- ncdim_def(name = 'ensemble', units = "adim",
-                          vals = 1 : dim(data$data)[which(dimname == 'member')],
+                          vals = 1 : dim(data)[which(dimname == 'member')],
                           longname = 'ensemble', create_dimvar = TRUE)
     dims_var[[list_pos]] <- dim_memb
     list_pos <- list_pos + 1
     }
   }
 
-  # Lead-time depends on the start date
-  nlt <- length(data$Dates$start)/n_sdates
 
   if (any(dimname == 'level')) {
     stop("Ask for saving 3Dim fields to the mantainer.")
@@ -160,49 +221,57 @@ list_pos <- 1
   for (i in 1 : n_datasets) {
     path <- file.path(destination, datasets[i], var_name)
     dir.create(path, recursive = TRUE)
-    startdate <- gsub("-", "", data$Datasets[[i]]$InitializationDates[[1]])
-    file_name <- paste0(var_name, "_", startdate, ".nc")
-    full_filename <- file.path(path, file_name)
-
-    data_dataset <- Subset(data$data, along = which(dimname == 'dataset'), indices = i)
-    standard_order <- c("lon", "lat", "member", "ftime")
-    change_names <- c("lon", "lat", "ensemble", "ftime")
-      for (j in 1 : n_sdates) {
-        n_data <- s2dverification::Subset(data_dataset,
-                                          along = which(dimname == 'sdate'),
-                                          indices = j, drop = TRUE)
-        pos_standard_order <- match( standard_order, names(dim(n_data)))
-        n_data <- aperm(n_data, pos_standard_order)
-
-        names(dim(n_data)) <- change_names
+    startdate <- gsub("-", "", startdates)
 
-    # Lead-time depends on the start date
-    # The correct time should be selected from $Dates$start
-    time_values <- as.Date(substr(data$Dates$start[(j * nlt - nlt + 1):(j * nlt)],
-                           1, 10))
-
-    if (any(dimname == 'time') | any(dimname == 'ftime')) {
-        dim_time <- ncdim_def(name = 'time', units = 'days since 1970-01-01',
-                              vals = as.numeric(time_values),
-                              longname = 'time', unlim = TRUE)
-          if (i == 1 & j == 1) {
-             dims_var[[list_pos]] <- dim_time
-             list_pos <- list_pos + 1
-          }
-    }
-    if (!is.character(attributes(data$Variable)$units)) {
-        units = attributes(data$Variable)$variable$units
-    } else {
-        units = attributes(data$Variable)$units
-    }
-        datanc  <- ncvar_def(name = var_name,
-                             units = units,
-                             dim = dims_var, missval = -99999)
-        file_nc <- nc_create(full_filename[j], datanc)
-        ncvar_put(file_nc, datanc, n_data)
-        ncatt_put(file_nc, datanc, 'coordinates', attr(data$lon, 'cdo_grid_name'))
-        ncatt_put(file_nc, datanc, 'projection', attr(data$lon, 'projection'))
-        nc_close(file_nc)
-      }
+    dim(startdate) <- c(sdate = length(startdate)) 
+    Apply(list(data, startdate, Dates),
+          target_dims = list(c('member', 'time', 'latitude', 'longitude'), 
+                             NULL, 'time'),
+          fun = .saveExp, var_name = var_name, units = units,
+          dims_var = dims_var, cdo_grid_name = cdo_grid_name, projection = projection,
+          destination = path)
   }
 }
+
+# data is an array with dimensions: member, time, lat, lon:
+# Dates is a vector of the dates for the time dimension
+# dims_var is a list with the ncdim_def of common variables in dataset: member, lat and lon:
+# data <- 1:(3 * 4 * 5 * 6)
+# dim(data) <- c(longitude = 3, latitude = 4, time = 5, member = 6)
+# var_name <- 'tas'
+# units <- 'K'
+# lon <- 1:3
+# lat <- 1:4
+# sdate = '19001101'
+# destination = '/esarchive/scratch/nperez/git/Flor/cstools/'
+# dims_var = list(ncdim_def(name = 'lon', units = 'degrees',
+#                          vals = as.vector(lon), longname = 'longitude'),
+#                ncdim_def(name = 'lat', units = 'degrees_north',
+#                         vals = as.vector(lat), longname = 'latitude'),
+#                ncdim_def(name = 'ensemble', units = "adim",
+#                          vals = 1 : 6,
+#                          longname = 'ensemble', create_dimvar = TRUE))
+#Dates <- as.Date(c("1900-11-01", "1900-12-01", "1901-01-01", "1901-02-01", "1901-03-01"))
+#.saveExp(data, sdate, Dates, var_name, units, dims_var, cdo_grid_name = 'r360x181', projection = 'none', destination)
+.saveExp <- function(data, sdate, Dates, var_name, units, dims_var, 
+                     cdo_grid_name, projection, destination) {
+  dim_names <- names(dim(data))
+  if (any(dim_names != c('longitude', 'latitude', 'member', 'time'))) {
+    data <- Reorder(data, c('longitude', 'latitude', 'member', 'time'))
+  } 
+  dim_time <- ncdim_def(name = 'time', units = 'days since 1970-01-01',
+                        vals = as.numeric(Dates),
+                        longname = 'time', unlim = TRUE)
+  list_pos = length(dims_var) + 1
+  dims_var[[list_pos]] <- dim_time
+  datanc  <- ncvar_def(name = var_name,
+                       units = units,
+                       dim = dims_var, missval = -99999)
+  file_name <- paste0(var_name, "_", sdate, ".nc")
+  full_filename <- file.path(destination, file_name)
+  file_nc <- nc_create(full_filename, datanc)
+  ncvar_put(file_nc, datanc, data)
+  ncatt_put(file_nc, datanc, 'coordinates', cdo_grid_name)
+  ncatt_put(file_nc, datanc, 'projection', projection)
+  nc_close(file_nc)
+}
diff --git a/man/SaveExp.Rd b/man/SaveExp.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..40ace2dbab6cc9bb1c84292ce7387f6a3e991dc0
--- /dev/null
+++ b/man/SaveExp.Rd
@@ -0,0 +1,73 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_SaveExp.R
+\name{SaveExp}
+\alias{SaveExp}
+\title{Save an experiment in a format compatible with CST_Load}
+\usage{
+SaveExp(
+  data,
+  lon,
+  lat,
+  Dataset,
+  var_name,
+  units,
+  startdates,
+  Dates,
+  cdo_grid_name,
+  projection,
+  destination
+)
+}
+\arguments{
+\item{data}{an multi-dimensional array with named dimensions (longitude, latitude, time, member, sdate)}
+
+\item{lon}{vector of logitud corresponding to the longitudinal dimension in data}
+
+\item{lat}{vector of latitud corresponding to the latitudinal dimension in data}
+
+\item{Dataset}{a vector of character string indicating the names of the datasets}
+
+\item{var_name}{a character string indicating the name of the variable to be saved}
+
+\item{units}{a character string indicating the units of the variable}
+
+\item{startdates}{a vector of dates indicating the initialization date of each simulations}
+
+\item{Dates}{a matrix of dates with two dimension 'time' and 'sdate'.}
+
+\item{cdo_grid_name}{a character string indicating the name of the grid e.g.: 'r360x181'}
+
+\item{projection}{a character string indicating the projection name}
+
+\item{destination}{a character string indicating the path where to store the NetCDF files}
+}
+\value{
+the function creates as many files as sdates per dataset. Each file could contain multiple members
+The path will be created with the name of the variable and each Datasets.
+}
+\description{
+This function is created for compatibility with CST_Load/Load for saving post-processed datasets such as those calibrated of downscaled with CSTools functions
+}
+\examples{
+\dontrun{
+data <- lonlat_data$exp$data
+lon <- lonlat_data$exp$lon
+lat <- lonlat_data$exp$lat
+Dataset <- 'XXX'
+var_name <- 'tas'
+units <- 'k'
+startdates <- lapply(1:length(lonlat_data$exp$Datasets),
+                    function(x) {
+                        lonlat_data$exp$Datasets[[x]]$InitializationDates[[1]]})[[1]]
+Dates <- lonlat_data$exp$Dates$start
+dim(Dates) <- c(time = length(Dates)/length(startdates), sdate = length(startdates))
+cdo_grid_name = attr(lonlat_data$exp$lon, 'cdo_grid_name')
+projection = attr(lonlat_data$exp$lon, 'projection')
+destination = './path/'
+SaveExp(data, lon, lat, Dataset, var_name, units, startdates, Dates,
+                   cdo_grid_name, projection, destination)
+}
+}
+\author{
+Perez-Zanon Nuria, \email{nuria.perez@bsc.es}
+}