diff --git a/example_scripts/exec_ecvs_seasonal_oper.R b/example_scripts/exec_ecvs_seasonal_oper.R
new file mode 100644
index 0000000000000000000000000000000000000000..18f2e4938cd8d48734782717b028c4775c891316
--- /dev/null
+++ b/example_scripts/exec_ecvs_seasonal_oper.R
@@ -0,0 +1,35 @@
+rm(list=ls())
+gc()
+setwd("/esarchive/scratch/nperez/git/auto-s2s")
+
+source("modules/Loading/Loading.R")
+source("modules/Calibration/Calibration.R")
+source("modules/Anomalies/Anomalies.R")
+source("modules/Skill/Skill.R")
+source("modules/Saving/Saving.R")
+source("modules/Visualization/Visualization.R")
+source("tools/prepare_outputs.R")
+
+# Read recipe
+args = commandArgs(trailingOnly = TRUE)
+recipe_file <- args[1]
+recipe <- read_atomic_recipe(recipe_file)
+## to test a single recipe:
+#recipe_file <- "recipe_ecvs_seasonal_oper.yml"
+#recipe <- prepare_outputs(recipe_file)
+
+# Load datasets
+data <- load_datasets(recipe)
+# Calibrate datasets
+data <- calibrate_datasets(recipe, data)
+# Compute skill metrics 
+skill_metrics <- compute_skill_metrics(recipe, data)
+# Compute percentiles and probability bins
+probabilities <- compute_probabilities(recipe, data)
+# Export all data to netCDF
+## TODO: Fix plotting
+# save_data(recipe, data, skill_metrics, probabilities)
+# Plot data
+plot_data(recipe, data, skill_metrics, probabilities, significance = T)
+
+
diff --git a/example_scripts/exec_units.R b/example_scripts/exec_units.R
new file mode 100644
index 0000000000000000000000000000000000000000..bcd443c5da0427f553c38cf4dec234873a9a16f1
--- /dev/null
+++ b/example_scripts/exec_units.R
@@ -0,0 +1,42 @@
+rm(list=ls())
+gc()
+setwd("/esarchive/scratch/nperez/git/auto-s2s")
+
+source("modules/Loading/Loading.R")
+source("modules/Calibration/Calibration.R")
+source("modules/Anomalies/Anomalies.R")
+source("modules/Skill/Skill.R")
+source("modules/Saving/Saving.R")
+source("modules/Visualization/Visualization.R")
+source("tools/prepare_outputs.R")
+
+# Read recipe
+#args = commandArgs(trailingOnly = TRUE)
+#recipe_file <- args[1]
+#recipe <- read_atomic_recipe(recipe_file)
+## to test a single recipe:
+ #  recipe_file <- "recipes/examples/recipe_tas_seasonal_units.yml"
+ #  recipe_file <- "recipes/examples/recipe_prlr_seasonal_units.yml"
+
+recipe <- prepare_outputs(recipe_file)
+
+# Load datasets
+data <- load_datasets(recipe)
+# Units transformation
+source("modules/Units/Units.R")
+test <- Units(recipe, data)
+# Calibrate datasets
+data <- calibrate_datasets(recipe, test)
+# Compute skill metrics 
+skill_metrics <- compute_skill_metrics(recipe, data)
+# Compute percentiles and probability bins
+probabilities <- compute_probabilities(recipe, data)
+# Export all data to netCDF
+## TODO: Fix plotting
+# save_data(recipe, data, skill_metrics, probabilities)
+# Plot data
+plot_data(recipe, data, skill_metrics, probabilities, significance = T)
+
+
+
+
diff --git a/example_scripts/test_seasonal.R b/example_scripts/test_seasonal.R
index 0c1c1d2b9ac751b7be30d2e911122dedd886dc4d..9734aeafb6629664e18a2eeb907e0d42bec5a78a 100644
--- a/example_scripts/test_seasonal.R
+++ b/example_scripts/test_seasonal.R
@@ -6,6 +6,7 @@
 
 # Load modules
 source("modules/Loading/Loading.R")
+source("modules/Units/Units.R")
 source("modules/Calibration/Calibration.R")
 source("modules/Anomalies/Anomalies.R")
 source("modules/Skill/Skill.R")
@@ -18,6 +19,8 @@ recipe <- prepare_outputs(recipe_file)
 
 # Load datasets
 data <- load_datasets(recipe)
+# Change units
+data <- Units(recipe, data)
 # Calibrate datasets
 data <- calibrate_datasets(recipe, data)
 # Compute anomalies
diff --git a/exec_units.R b/exec_units.R
new file mode 100644
index 0000000000000000000000000000000000000000..fa95767b12fcc68eec4f69d69b45e4a21bd66c27
--- /dev/null
+++ b/exec_units.R
@@ -0,0 +1,42 @@
+rm(list=ls())
+gc()
+setwd("/esarchive/scratch/nperez/git/auto-s2s")
+
+source("modules/Loading/Loading.R")
+source("modules/Calibration/Calibration.R")
+source("modules/Anomalies/Anomalies.R")
+source("modules/Skill/Skill.R")
+source("modules/Saving/Saving.R")
+source("modules/Visualization/Visualization.R")
+source("tools/prepare_outputs.R")
+
+# Read recipe
+#args = commandArgs(trailingOnly = TRUE)
+#recipe_file <- args[1]
+#recipe <- read_atomic_recipe(recipe_file)
+## to test a single recipe:
+ #  recipe_file <- "recipe_tas_seasonal_units.yml"
+ #  recipe_file <- "recipe_prlr_seasonal_units.yml"
+
+recipe <- prepare_outputs(recipe_file)
+
+# Load datasets
+data <- load_datasets(recipe)
+# Units transformation
+source("modules/Units/Units.R")
+test <- Units(recipe, data)
+# Calibrate datasets
+data <- calibrate_datasets(recipe, test)
+# Compute skill metrics 
+skill_metrics <- compute_skill_metrics(recipe, data)
+# Compute percentiles and probability bins
+probabilities <- compute_probabilities(recipe, data)
+# Export all data to netCDF
+## TODO: Fix plotting
+# save_data(recipe, data, skill_metrics, probabilities)
+# Plot data
+plot_data(recipe, data, skill_metrics, probabilities, significance = T)
+
+
+
+
diff --git a/modules/Loading/Loading.R b/modules/Loading/Loading.R
index d945d5af3547013c1ec55eac55a869569249f20c..8cc46dc7f5e3fd25b82d5846eca2254b997962b5 100644
--- a/modules/Loading/Loading.R
+++ b/modules/Loading/Loading.R
@@ -365,26 +365,26 @@ load_datasets <- function(recipe) {
     
     # Convert prlr from m/s to mm/day
     ## TODO: Make a unit conversion function
-    if (variable[[var_idx]] == "prlr") {
-      # Verify that the units are m/s and the same in obs and hcst
-      if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
-           (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) &&
-          ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
-           (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) {
-        info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.")
-        obs$data[, var_idx, , , , , , , ] <-
-          obs$data[, var_idx, , , , , , , ]*86400*1000
-        obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
-        hcst$data[, var_idx, , , , , , , ] <-
-          hcst$data[, var_idx, , , , , , , ]*86400*1000
-        hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
-        if (!is.null(fcst)) {
-          fcst$data[, var_idx, , , , , , , ] <- 
-            fcst$data[, var_idx, , , , , , , ]*86400*1000
-          fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
-        }
-      }
-    }
+#    if (variable[[var_idx]] == "prlr") {
+#      # Verify that the units are m/s and the same in obs and hcst
+#      if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
+#           (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) &&
+#          ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
+#           (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) {
+#        info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.")
+#        obs$data[, var_idx, , , , , , , ] <-
+#          obs$data[, var_idx, , , , , , , ]*86400*1000
+#        obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
+#        hcst$data[, var_idx, , , , , , , ] <-
+#          hcst$data[, var_idx, , , , , , , ]*86400*1000
+#        hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
+#        if (!is.null(fcst)) {
+#          fcst$data[, var_idx, , , , , , , ] <- 
+#            fcst$data[, var_idx, , , , , , , ]*86400*1000
+#          fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
+#        }
+#      }
+#    }
   }
   # Compute anomalies if requested
   # Print a summary of the loaded data for the user, for each object
diff --git a/modules/Units/R/deaccumulate.R b/modules/Units/R/deaccumulate.R
new file mode 100644
index 0000000000000000000000000000000000000000..e1fcd6f59cc246b4c99267aba35291e622d6e173
--- /dev/null
+++ b/modules/Units/R/deaccumulate.R
@@ -0,0 +1,6 @@
+deaccumulate <- function(data, ncores = NULL) {
+  data[[1]]$data <- Apply(list(data[[1]]$data), target_dim = 'time',
+                          fun = function(x) {
+                            c(x[1], diff(x))
+                          }, ncores = ncores)$output1
+}
diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R
new file mode 100644
index 0000000000000000000000000000000000000000..a24ef851559928b49537de62402375ac0c829ee0
--- /dev/null
+++ b/modules/Units/R/transform_units_precipitation.R
@@ -0,0 +1,120 @@
+transform_units_precipitation <- function(data, original_units, new_units,
+                                          var_name, freq, flux = FALSE, ncores = NULL,
+                                          var_index = 1) {
+  ## TODO: Hard-coded subsetting
+## TODO consider higher frequencies (e.g. 6hourly) 
+## could create a constant hours <- 24 or hours <- 6 and use the same code
+  if (original_units == "ms-1") {
+    if (new_units == "mm") {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] * 3600 * 24 * 1000
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm'
+    } else if (new_units == "m") {
+     data[[1]]$data[ , var_index, , , , , , , ] <- 
+       data[[1]]$data[ , var_index, , , , , , , ] * 3600 * 24
+     data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm'
+    } else if (new_units == "kgm-2s-1") {
+     data[[1]]$data[ , var_index, , , , , , , ] <- 
+       data[[1]]$data[ , var_index, , , , , , , ] * 1000 
+     data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'kg m-2 s-1'
+    } else {
+      stop(paste("Unknown transformation from", original_units, "to",
+                 new_units))
+    }
+  } else if (original_units == "mm") {
+    if (new_units == "ms-1") {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] / (3600 * 24 * 1000)
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm s-1'
+    } else if (new_units == "m") {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] / 1000
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm'
+    } else if (new_units == "kgm-2s-1") {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] / (3600 * 24 )
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'kg m-2 s-1'
+    } else {
+      stop(paste("Unknown transformation from", original_units, "to",
+                 new_units))
+    }
+  } else if (original_units == "m") {
+    if (new_units == "ms-1") {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] / (3600 * 24)
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm s-1'
+    } else if (new_units == "mm") {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] * 1000
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm'
+    } else if (new_units == "kgm-2s-1") {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] * 1000 / (3600 * 24) 
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'kg m-2 s-1'
+    } else {
+      stop(paste("Unknown transformation from", original_units, "to",
+                 new_units))
+    }
+  } else if (original_units == "kgm-2s-1") {
+    if (new_units == "ms-1") {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] / 1000
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm s-1'
+    } else if (new_units == "mm") {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] * 3600 * 24
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm'
+    } else if (new_units == "m") {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] * 3600 * 24 / 1000
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm'
+    } else {
+      stop(paste("Unknown transformation from", original_units, "to",
+                 new_units))
+    }
+  } else {
+    stop("Unknown precipitation units transformation")
+  }
+  if (flux) {
+    data[[1]]$attrs$Variable$metadata[[var_name]]$units <- paste0(
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units, "/day")
+  } else {
+    ## TODO: Shouldn't use time dimension, need it for Compute()
+    if (freq == "monthly_mean") { # could it not be mean?
+      time_pos <- which(lapply(data[[1]]$attrs$Variable$metadata[[var_name]]$dim,
+                               function(x) {x$name}) == 'time')
+      cal <- tolower(data[[1]]$attrs$Variable$metadata[[var_name]]$dim[[time_pos]]$calendar)
+      data_subset <- Subset(data[[1]]$data, along = "var", indices = var_index, drop = 'selected')
+      data[[1]]$data[ , var_index, , , , , , , ] <-
+        Apply(list(data_subset, data[[1]]$attrs$Dates),
+              target_dim = list(c('syear'), c('syear')),
+              extra_info = list(cal = cal, days_in_month = .days_in_month),
+              fun = function(x, y) {
+                date <- as.Date(y, "%Y-%m-%d")
+                num_days <- .days_in_month(date, cal = .cal)
+                res <- x * num_days
+              }, ncores = ncores)$output1
+    }
+  }
+  return(data)
+}
+
+.days_in_month <- function(x, cal) {
+  if (cal %in% c('gregorian', 'standard', 'proleptic_gregorian')) {
+    N_DAYS_IN_MONTHS <- lubridate:::N_DAYS_IN_MONTHS
+    if (leap_year(year(x))) {
+      N_DAYS_IN_MONTHS[2] <- N_DAYS_IN_MONTHS[2] + 1
+    }
+  } else if (cal %in% c('360', '360_day')) {
+    N_DAYS_IN_MONTHS <- rep(30, 12)
+    names(N_DAYS_IN_MONTHS) <- month.abb
+  } else if (cal %in% c('365', '365_day')) {
+    N_DAYS_IN_MONTHS <- lubridate:::N_DAYS_IN_MONTHS
+  } else {
+    stop("Unknown calendar")
+  } 
+  month_x <- month(x, label = TRUE, locale = "C")
+  n_days <- N_DAYS_IN_MONTHS[month_x]
+  n_days[month_x == "Feb" & leap_year(x)] <- 29L
+  return(n_days)
+}
diff --git a/modules/Units/R/transform_units_pressure.R b/modules/Units/R/transform_units_pressure.R
new file mode 100644
index 0000000000000000000000000000000000000000..9712e9fe54b82638fb201bf276cbdc3f871553c9
--- /dev/null
+++ b/modules/Units/R/transform_units_pressure.R
@@ -0,0 +1,33 @@
+transform_units_pressure <- function(data, original_units, new_units, var_name,
+                                     var_index = 1) {
+
+  ## TODO: Hard-coded subsetting
+  if (original_units == 'pa') {
+    if (new_units == 'hpa') {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] /100
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'hPa'
+    } else if (new_units == 'mb') {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] /100
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mb'
+    }
+  } else if (original_units == 'hpa') {
+    if (new_units == 'pa') {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] * 100
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "Pa"
+    } else if (new_units == "mb") {
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "mb"
+    }
+  } else if (original_units == "mb") {
+    if (new_units == 'pa') {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] * 100
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "Pa"
+    } else if (new_units == "hPa") {
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "hPa"
+    }  
+  }
+  return(data)
+}
diff --git a/modules/Units/R/transform_units_temperature.R b/modules/Units/R/transform_units_temperature.R
new file mode 100644
index 0000000000000000000000000000000000000000..366f0d34cee4a523845351d0b8ffea12c239e037
--- /dev/null
+++ b/modules/Units/R/transform_units_temperature.R
@@ -0,0 +1,17 @@
+transform_units_temperature <- function(data, original_units, new_units,
+                                        var_name, var_index = 1,
+                                        var_dim = "var") {
+  ## TODO: Hard-coded subsetting
+  if (original_units == 'c' & new_units == 'k') {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] + 273.15
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'K'
+  }
+  if (original_units == 'k' & new_units == 'c') {
+      data[[1]]$data[ , var_index, , , , , , , ] <- 
+        data[[1]]$data[ , var_index, , , , , , , ] - 273.15
+      data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "C"
+
+  }
+  return(data)
+}
diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R
new file mode 100644
index 0000000000000000000000000000000000000000..11e05f3116ca2870e742606db99276064592aa85
--- /dev/null
+++ b/modules/Units/R/units_transform.R
@@ -0,0 +1,37 @@
+# data is a s2dv_cube object
+# units as character
+units_transform <- function(data, orig_units, user_units, var_name, freq, flux = FALSE,
+                            ncores = NULL) {
+  ## TODO: Change how argument 'flux' works
+  for (i in 1:length(var_name)) {
+    if (!(orig_units[i] == user_units[i])) {
+      if (orig_units[i] %in% c("c", "k")) {
+        if (user_units[i] %in% c("c", "k")) {
+          data <- transform_units_temperature(data, orig_units[i], user_units[i],
+                                              var_name[i], var_index = i)
+        } else {
+          stop("Transformation temperature units not available.")
+        }
+      } else if (orig_units[i] %in% c("ms-1", "kgm-2s-1", "mm", "m")) {
+        if (user_units[i] %in% c("ms-1", "kgm-2s-1", "mm", "m")) {
+          data <- transform_units_precipitation(data, orig_units[i], user_units[i],
+                                                var_name[i], freq, flux, 
+                                                ncores = ncores, var_index = i)
+        } else {
+          stop("Transformation precipitation units not available.")
+        }
+      } else if (orig_units[i] %in% c("pa", "hpa", "mb")) {
+         if (user_units[i] %in% c("pa", "hpa", "mb")) {
+          data <- transform_units_pressure(data, orig_units[i], user_units[i],
+                                           var_name[i])
+        } else {
+          stop("Transformation precipitation units not available.")
+        }
+      } else {
+        stop(paste("Transformation unknown from", orig_units[i], "to", user_units[i]))
+      }
+    }
+  }
+  return(data)
+}
+
diff --git a/modules/Units/Units.R b/modules/Units/Units.R
new file mode 100644
index 0000000000000000000000000000000000000000..07e175436a66563fb81ec2df9cd2ffd0447c8b75
--- /dev/null
+++ b/modules/Units/Units.R
@@ -0,0 +1,86 @@
+# This function aims to convert units
+## 
+source("modules/Units/R/units_transform.R")
+source("modules/Units/R/transform_units_temperature.R")
+source("modules/Units/R/transform_units_precipitation.R")
+source("modules/Units/R/transform_units_pressure.R")
+# Requires recipe to include:
+# recipe$Analysis$Variables:
+# $freq
+# $flux (if precipitation)
+# $units (optional)
+Units <- function(recipe, data) {
+  # from recipe read the user defined units
+  # from data the original units
+  # deaccumulate option for CDS accumulated variables?
+  ncores <- recipe$Analysis$ncores
+  ## Do we need to convert other than ECVs?
+  var_names <- lapply(data, function(x) {x$attrs$Variable$varName})
+  freq <- recipe$Analysis$Variable$freq
+  if (is.null(recipe$Analysis$Variable$flux)) {
+    flux <- FALSE
+  } else {
+    flux <- recipe$Analysis$Variable$flux
+  }
+  orig_units <- list()
+  for (element in names(var_names)) {
+    orig_units[[element]] <- c()
+    for (x in var_names[[element]]) {
+      orig_units[[element]] <- c(orig_units[[element]], 
+                                 data[[element]]$attrs$Variable$metadata[[x]]$units)
+    }
+  }
+  if (is.null(recipe$Analysis$Variables$units)) {
+    user_units <- orig_units[[which(!is.null(orig_units))[1]]]
+  } else {
+    user_units <- recipe$Analysis$Variables$units
+    if (!is.null(names(user_units))) {
+      # Change to vector and ensure that the units are in the correct order
+      user_units <- unlist(user_units[var_names[[1]]])
+    } else if (length(var_names[[1]] > 1)) {
+      ## TODO: How to handle spaces in multi-var case?
+      user_units <- strsplit(recipe$Analysis$Variables$units, ", |,")[[1]]
+    }
+  }
+  # remove spaces, "**", "*" and "per" from units
+  user_units <- tolower(user_units) 
+  user_units <- gsub(" ", "", user_units)
+  user_units <- gsub("\\**", "", user_units)
+  user_units <- gsub("\\*", "", user_units)
+  orig_units <- lapply(orig_units, function(x) {
+                       x <- tolower(x)
+                       x <- gsub(" ", "", x)
+                       x <- gsub("\\**", "", x)
+                       x <- gsub("\\*", "", x)
+                       x <- unlist(ifelse(rep(length(x)==0, length(user_units)), user_units, x)) #when fcst is NULL
+                       }) 
+  ## TODO:
+  ## if "/" appears substitute by -1 in at the end of next unit. How to know?
+  # Check if all units are equal (if so, no conversion needed; else convert)
+  if (length(unique(c(unlist(orig_units), user_units))) == length(user_units)) {
+    info(recipe$Run$logger, "##### NO UNIT CONVERSION NEEDED #####")
+    res <- data
+  } else {
+    if (recipe$Run$filesystem == 'esarchive' &&
+        !(length(unique(orig_units)) == 1)) {
+      warn(recipe$Run$logger, 
+           paste("The units in", names(orig_units), "are not all equal.",
+                 "If this is not expected, please contact the ES data team."))
+    }
+    res <- sapply(1:length(data), function(x) {
+                  if (!all(orig_units[x] == user_units)) {
+                    result <- units_transform(data[x], 
+                                              orig_units = orig_units[[x]],
+                                              user_units = user_units, 
+                                              var_names[[x]], freq = freq,
+                                              flux = flux,
+                                              ncores = ncores)
+                  } else {
+                    result <- data[x]
+                  }
+                  return(result)
+                  }, simplify = TRUE) # instead of lapply to get the named list directly
+    info(recipe$Run$logger, "##### UNIT CONVERSION COMPLETE #####")
+  }
+  return(res)
+}
diff --git a/recipe_prlr_seasonal_units.yml b/recipe_prlr_seasonal_units.yml
new file mode 100644
index 0000000000000000000000000000000000000000..e03428ac44147c68af42d3bdfa0c120fe07d8c75
--- /dev/null
+++ b/recipe_prlr_seasonal_units.yml
@@ -0,0 +1,63 @@
+Description:
+  Author: nperez
+  Info: ECVs Oper ESS ECMWF SEAS5 Seasonal Forecast recipe (monthly mean, tas)
+
+Analysis:
+  Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal
+  Variables:
+    name: prlr
+    freq: monthly_mean
+    units: mm
+  Datasets:
+    System:
+      name: ECMWF-SEAS5.1 # Mandatory, str: system5c3s system21_m1 system35c3s
+    Multimodel: no # Mandatory, bool: Either yes/true or no/false
+    Reference:
+      name: ERA5 # Mandatory, str: Reference codename. See docu.
+  Time:
+    sdate: '0601' ## MMDD
+    fcst_year:  # Optional, int: Forecast year 'YYYY'
+    hcst_start: '2014' # Mandatory, int: Hindcast start year 'YYYY'
+    hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY'
+    ftime_min: 1 # Mandatory, int: First leadtime time step in months 
+    ftime_max: 6 # Mandatory, int: Last leadtime time step in months
+  Region:
+    latmin: 30 # Mandatory, int: minimum latitude
+    latmax: 50 # Mandatory, int: maximum latitude
+    lonmin: -10 # Mandatory, int: minimum longitude
+    lonmax: 10 # Mandatory, int: maximum longitude
+  Regrid:
+    method: bilinear # Mandatory, str: Interpolation method. See docu.
+    type: "to_system"
+    #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
+  Workflow:
+    Anomalies:
+      compute: no
+      cross_validation: no
+      save: none
+    Calibration:
+      method: evmos # Mandatory, str: Calibration method. See docu.
+      cross_validation: yes
+      save: none
+    Skill:
+      metric: mean_bias EnsCorr rpss crpss bss10 bss90
+      save: 'all'
+      cross_validation: yes
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds.
+      save: 'all'
+    Indicators:
+      index: no
+    Visualization:
+      plots: skill_metrics forecast_ensemble_mean most_likely_terciles
+      multi_panel: no
+      mask_terciles: both
+  ncores: 4 # Optional, int: number of cores, defaults to 1
+  remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+  logo: yes
+Run:
+  Loglevel: INFO
+  Terminal: TRUE
+  output_dir: /esarchive/scratch/nperez/cs_oper/
+  code_dir: /esarchive/scratch/nperez/git/s2s-suite/
diff --git a/recipe_tas_seasonal_units.yml b/recipe_tas_seasonal_units.yml
new file mode 100644
index 0000000000000000000000000000000000000000..d2b25321b8b6c43a2842c577d265469ba5131281
--- /dev/null
+++ b/recipe_tas_seasonal_units.yml
@@ -0,0 +1,63 @@
+Description:
+  Author: nperez
+  Info: ECVs Oper ESS ECMWF SEAS5 Seasonal Forecast recipe (monthly mean, tas)
+
+Analysis:
+  Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal
+  Variables:
+    name: tas
+    freq: monthly_mean
+    units: K
+  Datasets:
+    System:
+      name: ECMWF-SEAS5.1 # Mandatory, str: system5c3s system21_m1 system35c3s
+    Multimodel: no # Mandatory, bool: Either yes/true or no/false
+    Reference:
+      name: ERA5 # Mandatory, str: Reference codename. See docu.
+  Time:
+    sdate: '0601' ## MMDD
+    fcst_year: '2023' # Optional, int: Forecast year 'YYYY'
+    hcst_start: '2009' # Mandatory, int: Hindcast start year 'YYYY'
+    hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY'
+    ftime_min: 1 # Mandatory, int: First leadtime time step in months 
+    ftime_max: 6 # Mandatory, int: Last leadtime time step in months
+  Region:
+    latmin: 30 # Mandatory, int: minimum latitude
+    latmax: 50 # Mandatory, int: maximum latitude
+    lonmin: -10 # Mandatory, int: minimum longitude
+    lonmax: 10 # Mandatory, int: maximum longitude
+  Regrid:
+    method: bilinear # Mandatory, str: Interpolation method. See docu.
+    type: "to_system"
+    #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
+  Workflow:
+    Anomalies:
+      compute: no
+      cross_validation: no
+      save: none
+    Calibration:
+      method: evmos # Mandatory, str: Calibration method. See docu.
+      cross_validation: yes
+      save: none
+    Skill:
+      metric: mean_bias EnsCorr rpss crpss bss10 bss90
+      save: 'all'
+      cross_validation: yes
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds.
+      save: 'all'
+    Indicators:
+      index: no
+    Visualization:
+      plots: skill_metrics forecast_ensemble_mean most_likely_terciles
+      multi_panel: no
+      mask_terciles: both
+  ncores: 4 # Optional, int: number of cores, defaults to 1
+  remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+  logo: yes
+Run:
+  Loglevel: INFO
+  Terminal: TRUE
+  output_dir: /esarchive/scratch/nperez/cs_oper/
+  code_dir: /esarchive/scratch/nperez/git/s2s-suite/
diff --git a/recipes/atomic_recipes/recipe_test_multivar.yml b/recipes/atomic_recipes/recipe_test_multivar.yml
index 8eb3e962e0ae5f05ead34ffa8dd23aa88d9b3293..a8c2e71bced0d04b51f18dac5a2d17d1033b3b32 100644
--- a/recipes/atomic_recipes/recipe_test_multivar.yml
+++ b/recipes/atomic_recipes/recipe_test_multivar.yml
@@ -4,8 +4,9 @@ Description:
 Analysis:
   Horizon: Seasonal
   Variables:
-    name: tas prlr
+    name: tas, prlr
     freq: monthly_mean
+    units: {tas: C, prlr: mm}
   Datasets:
     System:
       name: ECMWF-SEAS5
@@ -18,7 +19,7 @@ Analysis:
     hcst_start: '1999'
     hcst_end: '2010'
     ftime_min: 1
-    ftime_max: 2
+    ftime_max: 1
   Region:
     latmin: -10
     latmax: 10
diff --git a/recipes/examples/recipe_ecvs_seasonal_oper.yml b/recipes/examples/recipe_ecvs_seasonal_oper.yml
new file mode 100644
index 0000000000000000000000000000000000000000..d47fd1593749cb07c04dd8166a73075ac2058d76
--- /dev/null
+++ b/recipes/examples/recipe_ecvs_seasonal_oper.yml
@@ -0,0 +1,73 @@
+Description:
+  Author: nperez
+  Info: ECVs Oper ESS ECMWF SEAS5 Seasonal Forecast recipe (monthly mean, tas)
+
+Analysis:
+  Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal
+  Variables:
+    - {name: tas, freq: monthly_mean}
+    - {name: prlr, freq: monthly_mean}
+  Datasets:
+    System:
+      - {name: ECMWF-SEAS5.1} # system21_m1 system35c3s
+    Multimodel: no # Mandatory, bool: Either yes/true or no/false
+    Reference:
+      - {name: ERA5} # Mandatory, str: Reference codename. See docu.
+  Time:
+    sdate: '0701' ## MMDD
+    fcst_year: '2023' # Optional, int: Forecast year 'YYYY'
+    hcst_start: '1993' # Mandatory, int: Hindcast start year 'YYYY'
+    hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY'
+    ftime_min: 1 # Mandatory, int: First leadtime time step in months 
+    ftime_max: 6 # Mandatory, int: Last leadtime time step in months
+  Region:
+   - {name: "UE", latmin: 20, latmax: 80, lonmin: -20, lonmax: 40}
+  Regrid:
+    method: bilinear # Mandatory, str: Interpolation method. See docu.
+    type: "to_system"
+    #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
+  Workflow:
+    Anomalies:
+      compute: no
+      cross_validation: no
+      save: none
+    Calibration:
+      method: evmos # Mandatory, str: Calibration method. See docu.
+      cross_validation: yes
+      save: none
+    Skill:
+      metric: mean_bias EnsCorr rpss crpss bss10 bss90
+      save: 'all'
+      cross_validation: yes
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds.
+      save: 'all'
+    Indicators:
+      index: no
+    Visualization:
+      plots: skill_metrics forecast_ensemble_mean most_likely_terciles
+      multi_panel: no
+      projection: lambert_europe
+  ncores: 4 # Optional, int: number of cores, defaults to 1
+  remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+  logo: yes
+Run:
+  Loglevel: INFO
+  Terminal: yes
+  filesystem: esarchive
+  output_dir: /esarchive/scratch/nperez/cs_oper/seasonal/ # replace with the directory where you want to save the outputs
+  code_dir: /esarchive/scratch/nperez/git/auto-s2s/ # replace with the directory where your code is
+  autosubmit: yes
+  # fill only if using autosubmit
+  auto_conf:
+    script: /esarchive/scratch/nperez/git/auto-s2s/exec_ecvs_seasonal_oper.R # replace with the path to your script
+    expid: a68v # replace with your EXPID
+    hpc_user: bsc32339 # replace with your hpc username
+    wallclock: 02:00 # hh:mm
+    processors_per_job: 4 
+    platform: nord3v2 
+    email_notifications: yes # enable/disable email notifications. Change it if you want to.
+    email_address: nuria.perez@bsc.es # replace with your email address
+    notify_completed: yes # notify me by email when a job finishes
+    notify_failed: yes # notify me by email when a job fails
diff --git a/recipes/examples/recipe_prlr_seasonal_oper.yml b/recipes/examples/recipe_prlr_seasonal_oper.yml
new file mode 100644
index 0000000000000000000000000000000000000000..fb6d1a1cd039df755abc45da9aae0d4486abdc5c
--- /dev/null
+++ b/recipes/examples/recipe_prlr_seasonal_oper.yml
@@ -0,0 +1,62 @@
+Description:
+  Author: nperez
+  Info: ECVs Oper ESS ECMWF SEAS5 Seasonal Forecast recipe (monthly mean, tas)
+
+Analysis:
+  Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal
+  Variables:
+    name: prlr
+    freq: monthly_mean
+  Datasets:
+    System:
+      name: ECMWF-SEAS5.1 # Mandatory, str: system5c3s system21_m1 system35c3s
+    Multimodel: no # Mandatory, bool: Either yes/true or no/false
+    Reference:
+      name: ERA5 # Mandatory, str: Reference codename. See docu.
+  Time:
+    sdate: '0701' ## MMDD
+    fcst_year: '2023' # Optional, int: Forecast year 'YYYY'
+    hcst_start: '1993' # Mandatory, int: Hindcast start year 'YYYY'
+    hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY'
+    ftime_min: 1 # Mandatory, int: First leadtime time step in months 
+    ftime_max: 6 # Mandatory, int: Last leadtime time step in months
+  Region:
+    latmin: 20 # Mandatory, int: minimum latitude
+    latmax: 80 # Mandatory, int: maximum latitude
+    lonmin: -20 # Mandatory, int: minimum longitude
+    lonmax: 40 # Mandatory, int: maximum longitude
+  Regrid:
+    method: bilinear # Mandatory, str: Interpolation method. See docu.
+    type: "to_system"
+    #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
+  Workflow:
+    Anomalies:
+      compute: no
+      cross_validation: no
+      save: none
+    Calibration:
+      method: evmos # Mandatory, str: Calibration method. See docu.
+      cross_validation: yes
+      save: none
+    Skill:
+      metric: mean_bias EnsCorr rpss crpss bss10 bss90
+      save: 'all'
+      cross_validation: yes
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds.
+      save: 'all'
+    Indicators:
+      index: no
+    Visualization:
+      plots: skill_metrics forecast_ensemble_mean most_likely_terciles
+      multi_panel: no
+      dots: both
+  ncores: 4 # Optional, int: number of cores, defaults to 1
+  remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+  logo: yes
+Run:
+  Loglevel: INFO
+  Terminal: TRUE
+  output_dir: /esarchive/scratch/nperez/cs_oper/
+  code_dir: /esarchive/scratch/nperez/git/s2s-suite/
diff --git a/recipes/examples/recipe_prlr_seasonal_units.yml b/recipes/examples/recipe_prlr_seasonal_units.yml
new file mode 100644
index 0000000000000000000000000000000000000000..e03428ac44147c68af42d3bdfa0c120fe07d8c75
--- /dev/null
+++ b/recipes/examples/recipe_prlr_seasonal_units.yml
@@ -0,0 +1,63 @@
+Description:
+  Author: nperez
+  Info: ECVs Oper ESS ECMWF SEAS5 Seasonal Forecast recipe (monthly mean, tas)
+
+Analysis:
+  Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal
+  Variables:
+    name: prlr
+    freq: monthly_mean
+    units: mm
+  Datasets:
+    System:
+      name: ECMWF-SEAS5.1 # Mandatory, str: system5c3s system21_m1 system35c3s
+    Multimodel: no # Mandatory, bool: Either yes/true or no/false
+    Reference:
+      name: ERA5 # Mandatory, str: Reference codename. See docu.
+  Time:
+    sdate: '0601' ## MMDD
+    fcst_year:  # Optional, int: Forecast year 'YYYY'
+    hcst_start: '2014' # Mandatory, int: Hindcast start year 'YYYY'
+    hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY'
+    ftime_min: 1 # Mandatory, int: First leadtime time step in months 
+    ftime_max: 6 # Mandatory, int: Last leadtime time step in months
+  Region:
+    latmin: 30 # Mandatory, int: minimum latitude
+    latmax: 50 # Mandatory, int: maximum latitude
+    lonmin: -10 # Mandatory, int: minimum longitude
+    lonmax: 10 # Mandatory, int: maximum longitude
+  Regrid:
+    method: bilinear # Mandatory, str: Interpolation method. See docu.
+    type: "to_system"
+    #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
+  Workflow:
+    Anomalies:
+      compute: no
+      cross_validation: no
+      save: none
+    Calibration:
+      method: evmos # Mandatory, str: Calibration method. See docu.
+      cross_validation: yes
+      save: none
+    Skill:
+      metric: mean_bias EnsCorr rpss crpss bss10 bss90
+      save: 'all'
+      cross_validation: yes
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds.
+      save: 'all'
+    Indicators:
+      index: no
+    Visualization:
+      plots: skill_metrics forecast_ensemble_mean most_likely_terciles
+      multi_panel: no
+      mask_terciles: both
+  ncores: 4 # Optional, int: number of cores, defaults to 1
+  remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+  logo: yes
+Run:
+  Loglevel: INFO
+  Terminal: TRUE
+  output_dir: /esarchive/scratch/nperez/cs_oper/
+  code_dir: /esarchive/scratch/nperez/git/s2s-suite/
diff --git a/recipes/examples/recipe_tas_seasonal_oper.yml b/recipes/examples/recipe_tas_seasonal_oper.yml
new file mode 100644
index 0000000000000000000000000000000000000000..75918265ba4d12a20226bd69511ef53eb1de1d4e
--- /dev/null
+++ b/recipes/examples/recipe_tas_seasonal_oper.yml
@@ -0,0 +1,62 @@
+Description:
+  Author: nperez
+  Info: ECVs Oper ESS ECMWF SEAS5 Seasonal Forecast recipe (monthly mean, tas)
+
+Analysis:
+  Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal
+  Variables:
+    name: tas
+    freq: monthly_mean
+  Datasets:
+    System:
+      name: ECMWF-SEAS5.1 # Mandatory, str: system5c3s system21_m1 system35c3s
+    Multimodel: no # Mandatory, bool: Either yes/true or no/false
+    Reference:
+      name: ERA5 # Mandatory, str: Reference codename. See docu.
+  Time:
+    sdate: '0601' ## MMDD
+    fcst_year: '2023' # Optional, int: Forecast year 'YYYY'
+    hcst_start: '1993' # Mandatory, int: Hindcast start year 'YYYY'
+    hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY'
+    ftime_min: 1 # Mandatory, int: First leadtime time step in months 
+    ftime_max: 6 # Mandatory, int: Last leadtime time step in months
+  Region:
+    latmin: 20 # Mandatory, int: minimum latitude
+    latmax: 80 # Mandatory, int: maximum latitude
+    lonmin: -20 # Mandatory, int: minimum longitude
+    lonmax: 40 # Mandatory, int: maximum longitude
+  Regrid:
+    method: bilinear # Mandatory, str: Interpolation method. See docu.
+    type: "to_system"
+    #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
+  Workflow:
+    Anomalies:
+      compute: no
+      cross_validation: no
+      save: none
+    Calibration:
+      method: evmos # Mandatory, str: Calibration method. See docu.
+      cross_validation: yes
+      save: none
+    Skill:
+      metric: mean_bias EnsCorr rpss crpss bss10 bss90
+      save: 'all'
+      cross_validation: yes
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds.
+      save: 'all'
+    Indicators:
+      index: no
+    Visualization:
+      plots: skill_metrics forecast_ensemble_mean most_likely_terciles
+      multi_panel: no
+      mask_terciles: both
+  ncores: 4 # Optional, int: number of cores, defaults to 1
+  remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+  logo: yes
+Run:
+  Loglevel: INFO
+  Terminal: TRUE
+  output_dir: /esarchive/scratch/nperez/cs_oper/
+  code_dir: /esarchive/scratch/nperez/git/s2s-suite/
diff --git a/recipes/examples/recipe_tas_seasonal_units.yml b/recipes/examples/recipe_tas_seasonal_units.yml
new file mode 100644
index 0000000000000000000000000000000000000000..d2b25321b8b6c43a2842c577d265469ba5131281
--- /dev/null
+++ b/recipes/examples/recipe_tas_seasonal_units.yml
@@ -0,0 +1,63 @@
+Description:
+  Author: nperez
+  Info: ECVs Oper ESS ECMWF SEAS5 Seasonal Forecast recipe (monthly mean, tas)
+
+Analysis:
+  Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal
+  Variables:
+    name: tas
+    freq: monthly_mean
+    units: K
+  Datasets:
+    System:
+      name: ECMWF-SEAS5.1 # Mandatory, str: system5c3s system21_m1 system35c3s
+    Multimodel: no # Mandatory, bool: Either yes/true or no/false
+    Reference:
+      name: ERA5 # Mandatory, str: Reference codename. See docu.
+  Time:
+    sdate: '0601' ## MMDD
+    fcst_year: '2023' # Optional, int: Forecast year 'YYYY'
+    hcst_start: '2009' # Mandatory, int: Hindcast start year 'YYYY'
+    hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY'
+    ftime_min: 1 # Mandatory, int: First leadtime time step in months 
+    ftime_max: 6 # Mandatory, int: Last leadtime time step in months
+  Region:
+    latmin: 30 # Mandatory, int: minimum latitude
+    latmax: 50 # Mandatory, int: maximum latitude
+    lonmin: -10 # Mandatory, int: minimum longitude
+    lonmax: 10 # Mandatory, int: maximum longitude
+  Regrid:
+    method: bilinear # Mandatory, str: Interpolation method. See docu.
+    type: "to_system"
+    #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
+  Workflow:
+    Anomalies:
+      compute: no
+      cross_validation: no
+      save: none
+    Calibration:
+      method: evmos # Mandatory, str: Calibration method. See docu.
+      cross_validation: yes
+      save: none
+    Skill:
+      metric: mean_bias EnsCorr rpss crpss bss10 bss90
+      save: 'all'
+      cross_validation: yes
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds.
+      save: 'all'
+    Indicators:
+      index: no
+    Visualization:
+      plots: skill_metrics forecast_ensemble_mean most_likely_terciles
+      multi_panel: no
+      mask_terciles: both
+  ncores: 4 # Optional, int: number of cores, defaults to 1
+  remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+  logo: yes
+Run:
+  Loglevel: INFO
+  Terminal: TRUE
+  output_dir: /esarchive/scratch/nperez/cs_oper/
+  code_dir: /esarchive/scratch/nperez/git/s2s-suite/
diff --git a/tools/data_summary.R b/tools/data_summary.R
index 714c5a0b39b61631fc0f0bab1abd31e0c661aec1..11f365cdda1d566e1a644a5aa3098537a077156f 100644
--- a/tools/data_summary.R
+++ b/tools/data_summary.R
@@ -34,8 +34,11 @@ data_summary <- function(data_cube, recipe) {
   info(recipe$Run$logger, paste0("Statistical summary of the data in ",
                                  object_name, ":"))
   for (var_index in 1:data_cube$dims[['var']]) {
+    variable_name <- data_cube$attrs$Variable$varName[var_index]
+    variable_units <- data_cube$attrs$Variable$metadata[[variable_name]]$units
     info(recipe$Run$logger,
-         paste("Variable:", data_cube$attrs$Variable$varName[var_index]))
+         paste0("Variable: ", variable_name,
+                " (units: ", variable_units, ")"))
     output_string <- capture.output(summary(Subset(data_cube$data,
                                                    along = "var",
                                                    indices = var_index)))