diff --git a/autosubmit/auto-multimodel.sh b/autosubmit/auto-multimodel.sh
new file mode 100644
index 0000000000000000000000000000000000000000..a9912666046bf3d8fc33643f764d864fa390963d
--- /dev/null
+++ b/autosubmit/auto-multimodel.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+############ AUTOSUBMIT INPUTS ############
+proj_dir=%PROJDIR%
+outdir=%common.OUTDIR%
+script=%common.SCRIPT%
+SPLIT=%SPLIT%
+###############################
+
+cd $proj_dir
+
+source split_to_recipe
+# atomic_recipe_number=$(printf "%02d" $CHUNK)
+atomic_recipe=${outdir}/logs/recipes/multimodel/atomic_recipe_sys-Multimodel${recipe}.yml
+
+source MODULES
+
+Rscript ${script} ${atomic_recipe}
diff --git a/autosubmit/auto-verification-CERISE.sh b/autosubmit/auto-verification-CERISE.sh
index caf2dd0ec8194b2868187eb26697a197003b5d1a..675b41d474064d53ca435d5681105a85f16d880a 100644
--- a/autosubmit/auto-verification-CERISE.sh
+++ b/autosubmit/auto-verification-CERISE.sh
@@ -9,8 +9,10 @@ CHUNK=%CHUNK%
 
 cd $proj_dir
 
-atomic_recipe_number=$(printf "%02d" $CHUNK)
-atomic_recipe=${outdir}/logs/recipes/atomic_recipe_${atomic_recipe_number}.yml
+source chunk_to_recipe
+
+# atomic_recipe_number=$(printf "%02d" $CHUNK)
+atomic_recipe=${outdir}/logs/recipes/atomic_recipe_${recipe}.yml
 
 ## Workaround to avoid bug in conda activate/source activate when running
 ## inside bash script
diff --git a/autosubmit/auto-verification.sh b/autosubmit/auto-verification.sh
index 0089e322d63e89a5578c98a6a64a1369b5b9b108..e909dbfb433ff6d59816734973547e918513bf76 100644
--- a/autosubmit/auto-verification.sh
+++ b/autosubmit/auto-verification.sh
@@ -9,8 +9,10 @@ CHUNK=%CHUNK%
 
 cd $proj_dir
 
-atomic_recipe_number=$(printf "%02d" $CHUNK)
-atomic_recipe=${outdir}/logs/recipes/atomic_recipe_${atomic_recipe_number}.yml
+source chunk_to_recipe
+
+# atomic_recipe_number=$(printf "%02d" $CHUNK)
+atomic_recipe=${outdir}/logs/recipes/atomic_recipe_${recipe}.yml
 
 source MODULES
 
diff --git a/autosubmit/conf_esarchive/jobs.yml b/autosubmit/conf_esarchive/jobs.yml
index a3c8934bf70f92d15a8e644a6f34947afcc28847..7e2a19480289b5da5a36820290b94caefa1938c2 100644
--- a/autosubmit/conf_esarchive/jobs.yml
+++ b/autosubmit/conf_esarchive/jobs.yml
@@ -5,7 +5,19 @@ JOBS:
     WALLCLOCK: 
     NOTIFY_ON: 
     PLATFORM: nord3v2
-    PROCESSORS: 
+    PROCESSORS:
+      # SPLITS: # n_atomic_recipes, number of atomic recipes
+  multimodel:
+    FILE: autosubmit/auto-multimodel.sh
+    RUNNING: once
+    WALLCLOCK:
+    NOTIFY_ON:
+    PLATFORM: nord3v2
+    PROCESSORS:
+    DEPENDENCIES:
+      verification:
+        SPLITS_FROM:
+    SPLITS: # n_atomic_recipes/n_models = n_multimodels
   scorecards:
     FILE: autosubmit/auto-scorecards.sh
     WALLCLOCK: 00:10
@@ -13,4 +25,5 @@ JOBS:
     NOTIFY_ON:
     PROCESSORS: 1
     DEPENDENCIES: verification
+  ## TODO: Add scorecards-multimodel with multimodel dependency? 
 
diff --git a/recipes/recipe_multimodel_seasonal.yml b/recipes/recipe_multimodel_seasonal.yml
index 4addb764a80d3c62dace6326939bc0610d3cc82c..b2442d1f27ec0da278392109a594cfab59893fba 100644
--- a/recipes/recipe_multimodel_seasonal.yml
+++ b/recipes/recipe_multimodel_seasonal.yml
@@ -64,16 +64,17 @@ Run:
   filesystem: esarchive
   output_dir: /esarchive/scratch/vagudets/auto-s2s-outputs/ # replace with the directory where you want to save the outputs
   code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/ # replace with the directory where your code is
-  autosubmit: no
+  script: 
+  autosubmit: yes
   # fill only if using autosubmit
   auto_conf:
-    script: /esarchive/scratch/cdelgado/gitlat/SUNSET/main_multimodel_seasonal.R # replace with the path to your script
-    expid: XXXX # replace with your EXPID
-    hpc_user: bsc32924 # replace with your hpc username
+    script: ./example_scripts/multimodel_seasonal.R # replace with the path to your script
+    expid: a6wq # replace with your EXPID
+    hpc_user: bsc32762 # 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: carlos.delgado@bsc.es # replace with your email address
+    email_address: victoria.agudetse@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/split.R b/split.R
index 23607aad63d02fe2178acc1aa4e327c02eadd910..faadafd68051b934cafcc2b8f2c2617dee3f849b 100755
--- a/split.R
+++ b/split.R
@@ -34,7 +34,10 @@ recipe <- prepare_outputs(recipe_file = arguments$recipe,
 run_parameters <- divide_recipe(recipe)
 
 if (!is.null(recipe$Run$autosubmit) && (recipe$Run$autosubmit)) {
-  write_autosubmit_conf(recipe, run_parameters$n_atomic_recipes)
+  write_autosubmit_conf(recipe = recipe, 
+                        nchunks = run_parameters$n_atomic_recipes,
+                        chunk_to_recipe = run_parameters$chunk_to_recipe,
+                        split_to_recipe = run_parameters$split_to_recipe)
   sink(arguments$tmpfile, append = FALSE)
   # Run with...
   cat("autosubmit")
diff --git a/tools/check_recipe.R b/tools/check_recipe.R
index 5061f14642cfef62eac541a485c79c1bb5760063..4bea7dab35108baff06d966b935c6cda7e200133 100644
--- a/tools/check_recipe.R
+++ b/tools/check_recipe.R
@@ -23,21 +23,21 @@ check_recipe <- function(recipe) {
   if (!("Analysis" %in% names(recipe))) {
     error(recipe$Run$logger,
           "The recipe must contain an element called 'Analysis'.")
-    error_status <- T
+    error_status <- TRUE
   }
 
   if (!all(PARAMS %in% names(recipe$Analysis))) {
     error(recipe$Run$logger, 
           paste0("The element 'Analysis' in the recipe must contain all of ",
                  "the following: ", paste(PARAMS, collapse = ", "), "."))
-    error_status <- T
+    error_status <- TRUE
   }
 
   if (!any(HORIZONS %in% tolower(recipe$Analysis$Horizon))) {
     error(recipe$Run$logger,
           paste0("The element 'Horizon' in the recipe must be one of the ",
                 "following: ", paste(HORIZONS, collapse = ", "), "."))
-    error_status <- T
+    error_status <- TRUE
   }
   # Check time settings
   if (tolower(recipe$Analysis$Horizon) == "seasonal") {
@@ -48,7 +48,7 @@ check_recipe <- function(recipe) {
             paste0("The element 'Time' in the recipe must contain all of the ",
                    "following: ", paste(TIME_SETTINGS_SEASONAL,
                                         collapse = ", "), "."))
-      error_status <- T
+      error_status <- TRUE
     }
   } else if (tolower(recipe$Analysis$Horizon) == "decadal") {
     archive <- read_yaml(ARCHIVE_DECADAL)[[recipe$Run$filesystem]]
@@ -57,7 +57,7 @@ check_recipe <- function(recipe) {
             paste0("The element 'Time' in the recipe must contain all of the ",
                    "following: ", paste(TIME_SETTINGS_DECADAL,
                                         collapse = ", "), "."))
-      error_status <- T
+      error_status <- TRUE
     }
   } else {
     archive <- NULL
@@ -82,81 +82,94 @@ check_recipe <- function(recipe) {
   # Check system names
   if (!is.null(archive)) {
     if (!all(recipe$Analysis$Datasets$System$name %in% 
-             c(names(archive$System),'Multimodel'))) {
+             c(names(archive$System), 'Multimodel'))) {
       error(recipe$Run$logger,
             "The specified System name was not found in the archive.")
-      error_status <- T
+      error_status <- TRUE
     }
     # Check reference names
     if (!all(recipe$Analysis$Datasets$Reference$name %in% 
         names(archive$Reference))) {
       error(recipe$Run$logger,
             "The specified Reference name was not found in the archive.")
-      error_status <- T
+      error_status <- TRUE
     }
   }
   # Check multimodel
-  if (!is.null(recipe$Analysis$Datasets$Multimodel) && 
-       !tolower(recipe$Analysis$Datasets$Multimodel) %in% c('no','false')){
-    if (!tolower(recipe$Analysis$Datasets$Multimodel$execute) %in% 
-        c('yes','true','no','false','both')){
-      error(recipe$Run$logger,
-            paste("The specified execution for the multimodel is not valid.",
-                  "Please specify yes/true, no/false or both."))
-      error_status <- T
-    }
-    if (!tolower(recipe$Analysis$Datasets$Multimodel$approach) %in% 
-        c('pooled')){ #,'mean','median')){
-      error(recipe$Run$logger,
-            paste("The specified approach for the multimodel is not valid.",
-                  "Please specify pooled.")) #, mean or median."))
-      error_status <- T
-    }
-    if (!tolower(recipe$Analysis$Datasets$Multimodel$createFrom) %in% 
-        c('calibration','anomalies','indicators')){
-      error(recipe$Run$logger,
-            paste("The specified 'createFrom' for the multimodel is not valid.",
-                  "Please specify Calibration, Anomalies, Indicators."))
-      error_status <- T
+  if (is.null(recipe$Analysis$Datasets$Multimodel) || 
+      (is.logical(recipe$Analysis$Datasets$Multimodel) &&
+       !(recipe$Analysis$Datasets$Multimodel))) {
+    recipe$Analysis$Datasets$Multimodel <- list(execute = FALSE)
+  }
+  if (tolower(recipe$Analysis$Datasets$Multimodel$execute) == 'false') {
+    multimodel <- FALSE
+  } else {
+    multimodel <- TRUE
+  }
+  MULTIMODEL_METHODS <- c("pooled") ## to be added: mean, median...
+  MULTIMODEL_CREATEFROM <- c("calibration", "anomalies", "indicators")
+  if (multimodel) {
+    if (!is.null(recipe$Analysis$Datasets$Multimodel)) {
+      if (!tolower(recipe$Analysis$Datasets$Multimodel$execute) %in% 
+          c('true', 'false', 'both')) {
+        error(recipe$Run$logger,
+              paste("The specified execution for the multimodel is not valid.",
+                    "Please specify yes/true, no/false or 'both'."))
+        error_status <- TRUE
+      }
+      if (!tolower(recipe$Analysis$Datasets$Multimodel$approach) %in% 
+          MULTIMODEL_METHODS) {
+        error(recipe$Run$logger,
+              paste("The specified approach for the multimodel is not valid.",
+                    "Please specify pooled.")) #, mean or median."))
+        error_status <- TRUE
+      }
+      if (!tolower(recipe$Analysis$Datasets$Multimodel$createFrom) %in% 
+          MULTIMODEL_CREATEFROM) {
+        error(recipe$Run$logger,
+              paste("The specified 'createFrom' for the multimodel is not valid.",
+                    "Please specify Calibration, Anomalies, Indicators."))
+        error_status <- TRUE
+      }
     }
   } else {
-    recipe$Analysis$Datasets$Multimodel <- 'no'
+      recipe$Analysis$Datasets$Multimodel <- FALSE
   }
   # Check ftime_min and ftime_max
   if ((!(recipe$Analysis$Time$ftime_min > 0)) ||
       (!is.integer(recipe$Analysis$Time$ftime_min))) {
     error(recipe$Run$logger,
           "The element 'ftime_min' must be an integer larger than 0.")
-    error_status <- T
+    error_status <- TRUE
   }
   if ((!(recipe$Analysis$Time$ftime_max > 0)) ||
       (!is.integer(recipe$Analysis$Time$ftime_max))) {
     error(recipe$Run$logger,
           "The element 'ftime_max' must be an integer larger than 0.")
-    error_status <- T
+    error_status <- TRUE
   }
   if (recipe$Analysis$Time$ftime_max < recipe$Analysis$Time$ftime_min) {
     error(recipe$Run$logger,
           "'ftime_max' cannot be smaller than 'ftime_min'.")
-    error_status <- T
+    error_status <- TRUE
   }
   # Check consistency of hindcast years
   if (!(as.numeric(recipe$Analysis$Time$hcst_start) %% 1 == 0) ||
       (!(recipe$Analysis$Time$hcst_start > 0))) {
     error(recipe$Run$logger,
           "The element 'hcst_start' must be a valid year.")
-    error_status <- T
+    error_status <- TRUE
   }
   if (!(as.numeric(recipe$Analysis$Time$hcst_end) %% 1 == 0) ||
       (!(recipe$Analysis$Time$hcst_end > 0))) {
     error(recipe$Run$logger,
           "The element 'hcst_end' must be a valid year.")
-    error_status <- T
+    error_status <- TRUE
   }
   if (recipe$Analysis$Time$hcst_end < recipe$Analysis$Time$hcst_start) {
     error(recipe$Run$logger,
           "'hcst_end' cannot be smaller than 'hcst_start'.")
-    error_status <- T
+    error_status <- TRUE
   }
   ## TODO: Is this needed?
   if (is.null(recipe$Analysis$Time$fcst_year) || 
@@ -201,18 +214,14 @@ check_recipe <- function(recipe) {
   if (length(recipe$Analysis$Regrid) != 2) {
     error(recipe$Run$logger,
           "The 'Regrid' element must specify the 'method' and 'type'.")
-    error_status <- T
+    error_status <- TRUE
   }
-  if (!is.null(recipe$Analysis$Datasets$Multimodel) && 
-      !tolower(recipe$Analysis$Datasets$Multimodel) %in% c('no','false')){
-    if (recipe$Analysis$Regrid$type == 'to_system' && 
-        tolower(recipe$Analysis$Datasets$Multimodel$execute) 
-        %in% c('both','yes','true')) {
+
+  if (recipe$Analysis$Regrid$type == 'to_system' && multimodel) {
       error(recipe$Run$logger,
             paste0("The 'Regrid$type' cannot be 'to_system' if ",
                    "'Multimodel$execute' is yes/true or both."))
-      error_status <- T
-    }
+      error_status <- TRUE
   }
   # TODO: Add Workflow checks?
   # ...
@@ -220,7 +229,7 @@ check_recipe <- function(recipe) {
   if (length(recipe$Analysis$Horizon) > 1) {
     error(recipe$Run$logger,
           "Only one single Horizon can be specified in the recipe")
-    error_status <- T
+    error_status <- TRUE
   }
 
   ## TODO: Refine this
@@ -252,7 +261,7 @@ check_recipe <- function(recipe) {
         error(recipe$Run$logger,
               paste0("There must be 4 elements in 'Region': ",
                    paste(LIMITS, collapse = ", "), "."))
-        error_status <- T
+        error_status <- TRUE
       }
     }
     if (length(recipe$Analysis$Region) > 1) {
@@ -269,7 +278,7 @@ check_recipe <- function(recipe) {
     error(recipe$Run$logger,
           paste0("There must be 4 elements in 'Region': ",
                  paste(LIMITS, collapse = ", "), "."))
-    error_status <- T
+    error_status <- TRUE
   }
     ## TODO: Implement multiple regions
     # nregions <- length(recipe$Analysis$Region)
@@ -280,7 +289,7 @@ check_recipe <- function(recipe) {
     #           paste0("Each region defined in element 'Region' ",
     #                  "should have 4 elements: ",
     #                  paste(limits, collapse = ", "), "."))
-    #     error_status <- T
+    #     error_status <- TRUE
     #   }
     #   if (length(recipe$Analysis$Region) > 1) {
     #     if (!("name" %in% names(recipe$Analysis$Region[[i]]))) {
@@ -307,7 +316,7 @@ check_recipe <- function(recipe) {
     if (is.null(recipe$Analysis$Workflow$Calibration$method)) {
       error(recipe$Run$logger,
             "The 'Calibration' element 'method' must be specified.")
-      error_status <- T
+      error_status <- TRUE
     }
     SAVING_OPTIONS_CALIB <- c("all", "none", "exp_only", "fcst_only")
     if ((is.null(recipe$Analysis$Workflow$Calibration$save)) ||
@@ -316,7 +325,7 @@ check_recipe <- function(recipe) {
             paste0("Please specify which Calibration module outputs you want ",
                    "to save with the 'save' parameter. The options are: ",
                    paste(SAVING_OPTIONS_CALIB, collapse = ", "), "."))
-      error_status <- T
+      error_status <- TRUE
     }
   }
   # Anomalies
@@ -325,12 +334,12 @@ check_recipe <- function(recipe) {
     if (is.null(recipe$Analysis$Workflow$Anomalies$compute)) {
       error(recipe$Run$logger,
             "Parameter 'compute' must be defined under 'Anomalies'.")
-      error_status <- T
+      error_status <- TRUE
     } else if (!(is.logical(recipe$Analysis$Workflow$Anomalies$compute))) {
       error(recipe$Run$logger,
             paste("Parameter 'Anomalies:compute' must be a logical value",
                   "(True/False or yes/no)."))
-      error_status <- T
+      error_status <- TRUE
     } else if ((recipe$Analysis$Workflow$Anomalies$compute)) {
       # Cross-validation check
       if (!is.logical(recipe$Analysis$Workflow$Anomalies$cross_validation)) {
@@ -338,7 +347,7 @@ check_recipe <- function(recipe) {
               paste("If anomaly computation is requested, parameter",
                     "'cross_validation' must be defined under 'Anomalies',
                     and it must be a logical value (True/False or yes/no)."))
-        error_status <- T
+        error_status <- TRUE
       }
       # Saving checks
       SAVING_OPTIONS_ANOM <- c("all", "none", "exp_only", "fcst_only")
@@ -348,7 +357,7 @@ check_recipe <- function(recipe) {
               paste0("Please specify which Anomalies module outputs you want ",
                      "to save with the 'save' parameter. The options are: ",
                      paste(SAVING_OPTIONS_ANOM, collapse = ", "), "."))
-        error_status <- T
+        error_status <- TRUE
       }
     }
   }
@@ -375,21 +384,21 @@ check_recipe <- function(recipe) {
               paste0("The type of Downscaling request in the recipe is not ",
             	   "available. It must be one of the following: ",
             	   paste(DOWNSCAL_TYPES, collapse = ", "), "."))
-        error_status <- T
+        error_status <- TRUE
       }
       if ((downscal_params$type %in% c("int", "intbc", "intlr", "logreg")) &&
           (is.null(downscal_params$target_grid))) {
         error(recipe$Run$logger,
               paste("A target grid is required for the downscaling method",
             	  "requested in the recipe."))
-        error_status <- T
+        error_status <- TRUE
       }
       if (downscal_params$type == "int") {
         if (is.null(downscal_params$int_method)) {
           error(recipe$Run$logger,
                 paste("Downscaling type 'int' was requested, but no",
             	    "interpolation method is provided in the recipe."))
-          error_status <- T
+          error_status <- TRUE
         }
       } else if (downscal_params$type %in% 
                  c("int", "intbc", "intlr", "logreg")) {
@@ -398,32 +407,32 @@ check_recipe <- function(recipe) {
                 paste("Downscaling type", downscal_params$type, 
                       "was requested in the recipe, but no",
             	      "interpolation method is provided."))
-          error_status <- T
+          error_status <- TRUE
         }
       } else if (downscal_params$type == "intbc") {
         if (is.null(downscal_params$bc_method)) {
           error(recipe$Run$logger,
                 paste("Downscaling type 'intbc' was requested in the recipe, but",
             	    "no bias correction method is provided."))
-          error_status <- T
+          error_status <- TRUE
         } else if (!(downscal_params$bc_method %in% BC_METHODS)) {
           error(recipe$Run$logger,
                 paste0("The accepted Bias Correction methods for the downscaling",
             	     " module are: ", paste(BC_METHODS, collapse = ", "), "."))
-          error_status <- T
+          error_status <- TRUE
         }
       } else if (downscal_params$type == "intlr") {
         if (length(downscal_params$lr_method) == 0) {
           error(recipe$Run$logger,
                 paste("Downscaling type 'intlr' was requested in the recipe, but",
             	    "no linear regression method was provided."))
-          error_status <- T
+          error_status <- TRUE
         } else if (!(downscal_params$lr_method %in% LR_METHODS)) {
           error(recipe$Run$logger,
                 paste0("The accepted linear regression methods for the",
             	     " downscaling module are: ",
             	     paste(LR_METHODS, collapse = ", "), "."))
-          error_status <- T
+          error_status <- TRUE
         }
       } else if (downscal_params$type == "analogs") {
         if (is.null(downscal_params$nanalogs)) {
@@ -436,19 +445,19 @@ check_recipe <- function(recipe) {
           error(recipe$Run$logger,
                 paste("Downscaling type 'logreg' was requested in the recipe, but",
             	    "no interpolation method was provided."))
-          error_status <- T
+          error_status <- TRUE
         }
         if (is.null(downscal_params$log_reg_method)) {
           error(recipe$Run$logger,
                 paste("Downscaling type 'logreg' was requested in the recipe,",
             	    "but no logistic regression method is provided."))
-          error_status <- T
+          error_status <- TRUE
         } else if (!(downscal_params$log_reg_method %in% LOGREG_METHODS)) {
           error(recipe$Run$logger,
                 paste0("The accepted logistic regression methods for the ",
             	     "downscaling module are: ",
             	     paste(LOGREG_METHODS, collapse = ", "), "."))
-          error_status <- T
+          error_status <- TRUE
         }
       }
     }
@@ -462,12 +471,12 @@ check_recipe <- function(recipe) {
       error(recipe$Run$logger,
             paste0("Indices uses Anomalies as input, but Anomalies are missing",
                    "in the recipe."))
-      error_status <- T
+      error_status <- TRUE
     } else if (!(recipe$Analysis$Workflow$Anomalies$compute)) {
         error(recipe$Run$logger,
               paste0("Indices uses Anomalies as input, but the parameter",
                      "'Anomalies:compute' is set as no/False."))
-      error_status <- T
+      error_status <- TRUE
     }
     recipe_indices <- tolower(names(recipe$Analysis$Workflow$Indices))
     if (!all(recipe_indices %in% indices)) {
@@ -475,7 +484,7 @@ check_recipe <- function(recipe) {
             paste0("Some of the indices under 'Indices' are not available.",
                    "The available Indices are: 'NAO', 'Nino1+2', 'Nino3', ",
                    "'Nino3.4' and 'Nino4'."))
-      error_status <- T
+      error_status <- TRUE
     }
     # Check that variables correspond with indices requested
     if (("nao" %in% recipe_indices) &&
@@ -484,7 +493,7 @@ check_recipe <- function(recipe) {
             paste0("It is not possible to compute the NAO with some of the ",
                    "variables requested. To compute the NAO, please make sure",
                    "your recipe requests only psl and/or z500."))
-      error_status <- T
+      error_status <- TRUE
     }
     if ((any(nino_indices %in% recipe_indices)) &&
         (!all(recipe_variables %in% c("tos", "sst")))) {
@@ -492,7 +501,7 @@ check_recipe <- function(recipe) {
             paste0("It is not possible to compute El Nino indices with some ",
                    "of the variables requested. To compute El Nino, please ",
                    "make sure your recipe requests only tos."))
-      error_status <- T
+      error_status <- TRUE
     }
   }
 
@@ -506,7 +515,7 @@ check_recipe <- function(recipe) {
     if (is.null(recipe$Analysis$Workflow$Skill$metric)) {
     error(recipe$Run$logger,
           "Parameter 'metric' must be defined under 'Skill'.")
-    error_status <- T
+    error_status <- TRUE
     } else {
       requested_metrics <- strsplit(recipe$Analysis$Workflow$Skill$metric, 
                                     ", | |,")[[1]]
@@ -515,7 +524,7 @@ check_recipe <- function(recipe) {
               paste0("Some of the metrics requested under 'Skill' are not ",
                      "available in SUNSET. Check the documentation to see the ",
                      "full list of accepted skill metrics."))
-        error_status <- T
+        error_status <- TRUE
       }
     }
     # Saving checks
@@ -526,7 +535,7 @@ check_recipe <- function(recipe) {
             paste0("Please specify whether you want to save the Skill metrics ",
                    "with the 'save' parameter. The options are: ",
                    paste(SAVING_OPTIONS_SKILL, collapse = ", "), "."))
-      error_status <- T
+      error_status <- TRUE
     }
   }
 
@@ -535,12 +544,12 @@ check_recipe <- function(recipe) {
     if (is.null(recipe$Analysis$Workflow$Probabilities$percentiles)) {
       error(recipe$Run$logger,
             "Parameter 'percentiles' must be defined under 'Probabilities'.")
-      error_status <- T
+      error_status <- TRUE
     } else if (!is.list(recipe$Analysis$Workflow$Probabilities$percentiles)) {
       error(recipe$Run$logger,
             paste("Parameter 'Probabilities:percentiles' expects a list.",
                   "See documentation in the wiki for examples."))
-      error_status <- T
+      error_status <- TRUE
     }
     # Saving checks
     SAVING_OPTIONS_PROBS <- c("all", "none", "bins_only", "percentiles_only")
@@ -551,7 +560,7 @@ check_recipe <- function(recipe) {
                    "and probability bins with the 'save' parameter. The ",
                    "options are: ",
                    paste(SAVING_OPTIONS_PROBS, collapse = ", "), "."))
-      error_status <- T
+      error_status <- TRUE
     }
   }
 
@@ -563,7 +572,7 @@ check_recipe <- function(recipe) {
     if (is.null(recipe$Analysis$Workflow$Visualization$plots)) {
       error(recipe$Run$logger,
             "The 'plots' element must be defined under 'Visualization'.")
-      error_status <- T
+      error_status <- TRUE
     } else {
       plots <- strsplit(recipe$Analysis$Workflow$Visualization$plots,
                         ", | |,")[[1]]
@@ -571,7 +580,7 @@ check_recipe <- function(recipe) {
         error(recipe$Run$logger,
               paste0("The options available for the plots are: ",
                      paste(PLOT_OPTIONS, collapse = ", "), "."))
-        error_status <- T
+        error_status <- TRUE
       }
     }
     # Check multi_panel option
@@ -584,7 +593,7 @@ check_recipe <- function(recipe) {
       error(recipe$Run$logger,
             paste0("Parameter 'Visualization:multi_panel' must be a logical ",
                    "value: either 'yes/True' or 'no/False'"))
-      error_status <- T
+      error_status <- TRUE
     }
     # Check projection
     if (is.null(recipe$Analysis$Workflow$Visualization$projection)) {
@@ -603,7 +612,7 @@ check_recipe <- function(recipe) {
         error(recipe$Run$logger,
               paste0("Parameter Visualization:mask_terciles must be one of: ",
                      "yes/True, no/False, 'both'"))
-        error_status <- T
+        error_status <- TRUE
       }
       if (is.null(recipe$Analysis$Workflow$Visualization$dots)) {
         warn(recipe$Run$logger,
@@ -622,7 +631,7 @@ check_recipe <- function(recipe) {
     if (is.null(recipe$Analysis$Workflow$Scorecards$metric)) {
     error(recipe$Run$logger,
           "Parameter 'metric' must be defined under 'Scorecards'.")
-    error_status <- T
+    error_status <- TRUE
     } else {
       sc_metrics <- strsplit(recipe$Analysis$Workflow$Scorecards$metric, 
                              ", | |,")[[1]]
@@ -630,7 +639,7 @@ check_recipe <- function(recipe) {
         error(recipe$Run$logger,
               paste0("All of the metrics requested under 'Scorecards' must ",
                      "be requested in the 'Skill' section."))
-        error_status <- T
+        error_status <- TRUE
       }
     }
   }
@@ -649,28 +658,28 @@ check_recipe <- function(recipe) {
     error(recipe$Run$logger, paste("Recipe element 'Run' must contain",
                                    "all of the following fields:",
                                    paste(RUN_FIELDS, collapse=", "), "."))
-    error_status <- T
+    error_status <- TRUE
   }
   if (!is.character(recipe$Run$output_dir)) {
     error(recipe$Run$logger,
           paste("The Run element 'output_dir' in", recipe$name, "file",
                 "should be a character string indicating the path where",
                 "the outputs should be saved."))
-    error_status <- T
+    error_status <- TRUE
   }
   if (!is.character(recipe$Run$code_dir)) {
     error(recipe$Run$logger,
           paste("The Run element 'code_dir' in", recipe$name, "file ",
                 "should be a character string indicating the path",
                 "where the code is."))
-    error_status <- T
+    error_status <- TRUE
   }
   if (!is.logical(recipe$Run$Terminal)) {
     error(recipe$Run$logger,
           paste("The Run element 'Terminal' in", recipe$name, "file ",
                 "should be a boolean value indicating whether or not to",
                 "print the logs in the terminal."))
-    error_status <- T
+    error_status <- TRUE
   }
   ## TODO: Review this case, since default value is allowed
   if (!is.character(recipe$Run$Loglevel) || 
@@ -679,7 +688,7 @@ check_recipe <- function(recipe) {
           paste("The Run element 'Loglevel' in", recipe$name, "file",
                 "should be a character string specifying one of the levels available:",
                 paste0(LOG_LEVELS, collapse='/')))
-    error_status <- T
+    error_status <- TRUE
   }
 
   # ---------------------------------------------------------------------
@@ -701,22 +710,22 @@ check_recipe <- function(recipe) {
     if (!("auto_conf" %in% names(recipe$Run))) {
       error(recipe$Run$logger,
             "The 'auto_conf' is missing from the 'Run' section of the recipe.")
-      error_status <- T
+      error_status <- TRUE
     } else if (!all(AUTO_PARAMS %in% names(recipe$Run$auto_conf))) {
       error(recipe$Run$logger,
             paste0("The element 'Run:auto_conf' must contain all of the ",
                    "following: ", paste(AUTO_PARAMS, collapse = ", "), "."))
-      error_status <- T
+      error_status <- TRUE
     }
     # Check that the script is not NULL and exists
     if (is.null(recipe$Run$auto_conf$script)) {
       error(recipe$Run$logger,
             "A script must be provided to run the recipe with autosubmit.")
-      error_status <- T
+      error_status <- TRUE
     } else if (!file.exists(recipe$Run$auto_conf$script)) {
       error(recipe$Run$logger,
             "Could not find the file for the script in 'auto_conf'.")
-      error_status <- T
+      error_status <- TRUE
     }
     # Check that the experiment ID exists
     if (is.null(recipe$Run$auto_conf$expid)) {
@@ -728,30 +737,30 @@ check_recipe <- function(recipe) {
       error(recipe$Run$logger,
             paste("autosubmit expid -H", auto_specs$platform,
                   "-d <Description>"))
-      error_status <- T
+      error_status <- TRUE
     } else if (!dir.exists(paste0(auto_specs$experiment_dir,
                                   recipe$Run$auto_conf$expid))) {
       error(recipe$Run$logger,
             paste0("No folder in ", auto_specs$experiment_dir,
                    " for the EXPID", recipe$Run$auto_conf$expid,
                    ". Please make sure it is correct."))
-      error_status <- T
+      error_status <- TRUE
     }
     if ((recipe$Run$auto_conf$email_notifications) &&
         (is.null(recipe$Run$auto_conf$email_address))) {
       error(recipe$Run$logger,
             "Autosubmit notifications are enabled but email address is empty!")
-      error_status <- T
+      error_status <- TRUE
     }
     if (is.null(recipe$Run$auto_conf$hpc_user)) {
       error(recipe$Run$logger,
             "The 'Run:auto_conf:hpc_user' field can not be empty.")
-      error_status <- T
+      error_status <- TRUE
     } else if ((recipe$Run$filesystem == "esarchive") &&
                (!substr(recipe$Run$auto_conf$hpc_user, 1, 5) == "bsc32")) {
       error(recipe$Run$logger,
             "Please check your hpc_user ID. It should look like: 'bsc32xxx'")
-      error_status <- T
+      error_status <- TRUE
     }
   }
 
diff --git a/tools/divide_recipe.R b/tools/divide_recipe.R
index 214362d24f30a4cf3a70ad527b6ad6c9a0c6f4e3..3b2e6eee3399fec5be8daebe097222951f65123e 100644
--- a/tools/divide_recipe.R
+++ b/tools/divide_recipe.R
@@ -171,33 +171,48 @@ divide_recipe <- function(recipe) {
     }
   } # Rest of horizons
   # Save all recipes in separate YAML files
+  chunk <- 1
+  split <- 1
+  chunk_to_recipe <- list()
+  split_to_recipe <- list()
+  total_models <- length(recipe$Analysis$Datasets$System)
   for (reci in 1:length(all_recipes)) {
-    ## TODO: Sort dependencies
-    # if (reci < 10) {
-    #   recipe_number <- paste0("0", reci)
-    # } else {
-    #   recipe_number <- reci
-    # }
-    recipe_number <- paste0('var-',all_recipes[[reci]]$Analysis$Variables$name,
-                            '_sys-',gsub('\\.', '', all_recipes[[reci]]$Analysis$Datasets$System$name),
-                            '_ref-',all_recipes[[reci]]$Analysis$Datasets$Reference$name,
-                            '_reg-',all_recipes[[reci]]$Analysis$Region$name,
-                            '_sdate-',all_recipes[[reci]]$Analysis$Time$sdate)
+    ## TODO: Document
+    recipe_model <- paste0("sys-",
+                           gsub('\\.', '', all_recipes[[reci]]$Analysis$Datasets$System$name))
+    # 
+    recipe_split <- paste0("_ref-", all_recipes[[reci]]$Analysis$Datasets$Reference$name,
+                           "_var-", all_recipes[[reci]]$Analysis$Variables$name,
+                           "_reg-", all_recipes[[reci]]$Analysis$Region$name,
+                           "_sdate-", all_recipes[[reci]]$Analysis$Time$sdate)
+    recipe_name <- paste0(recipe_model, recipe_split)
+    
+
     if (all_recipes[[reci]]$Analysis$Datasets$System$name == 'Multimodel') {
       recipe_dir <- paste0(recipe$Run$output_dir, "/logs/recipes/multimodel/")
+      split_to_recipe[split] <- recipe_split
+      split <- split + 1
     } else {
       recipe_dir <- paste0(recipe$Run$output_dir, "/logs/recipes/")
+      chunk_to_recipe[chunk] <- recipe_name
+      chunk <- chunk + 1
     }
     write_yaml(all_recipes[[reci]], 
-               paste0(recipe_dir, "atomic_recipe_", recipe_number, ".yml"))
+               paste0(recipe_dir, "atomic_recipe_", recipe_name, ".yml"))
   }
+
+  # Print information for user
   info(recipe$Run$logger,
-       paste("The main recipe has been divided into", length(all_recipes),
-             "atomic recipes."))
+       paste("The main recipe has been divided into", length(chunk_to_recipe),
+             "single model atomic recipes, plus", length(split_to_recipe),
+             "multi-model atomic recipes."))
   text <- paste0("Check output directory ", recipe$Run$output_dir,
                  "/logs/recipes/ to see all the individual atomic recipes.")
   info(recipe$Run$logger, text)
   ## TODO: Change returns?
-  return(list(n_atomic_recipes = length(all_recipes),
-              outdir = recipe$Run$output_dir))
+  return(list(n_atomic_recipes = length(chunk_to_recipe), # length(all_recipes)
+              outdir = recipe$Run$output_dir,
+              chunk_to_recipe = chunk_to_recipe,
+              split_to_recipe = split_to_recipe))
 }
+
diff --git a/tools/write_autosubmit_conf.R b/tools/write_autosubmit_conf.R
index 95ca93f0f09f10f8b8b67a5076f24918c3d6db7e..0bffaa4387b8b1c4205e59c2390e20c5a5f74a8a 100644
--- a/tools/write_autosubmit_conf.R
+++ b/tools/write_autosubmit_conf.R
@@ -1,5 +1,19 @@
-# Function to write autosubmit configuration from an Auto-S2S recipe
-write_autosubmit_conf <- function(recipe, nchunks) {
+# Function to write autosubmit configuration from an SUNSET recipe. The function
+# reads the corresponding AS configuration file templates and fills them with
+# the information needed to run the experiment. The modified configuration
+# files are saved in the `conf/` folder of the Autosubmit experimet.
+# 
+# recipe: the SUNSET recipe
+# nchunks: the number of 'chunks' to be processed by Autosubmit, as returned
+#   by divide_recipe().
+# chunk_to_recipe: list with the correspondence between the chunk number and 
+#   the name of the atomic recipe, as returned by divide_recipe().
+# split_to_recipe: list with the correspondence between the split number and
+#   the name of the multi-model atomic recipe, as returned by divide_recipe().
+
+write_autosubmit_conf <- function(recipe, nchunks,
+                                  chunk_to_recipe,
+                                  split_to_recipe) {
   # Experiment ID
   expid <- recipe$Run$auto_conf$expid
   # Directory with the experiment templates
@@ -8,6 +22,7 @@ write_autosubmit_conf <- function(recipe, nchunks) {
   auto_specs <- read_yaml("conf/autosubmit.yml")[[recipe$Run$filesystem]]
   # Output directory
   dest_dir <- paste0(auto_specs$experiment_dir, expid, "/conf/")
+  proj_dir <- paste0(auto_specs$experiment_dir, expid, "/proj/auto-s2s/")
   # Modify the configuration files according to the info in the recipe
   for (file in list.files(template_dir)) {
     conf_type <- strsplit(file, split = "[.]")[[1]][1]
@@ -35,6 +50,12 @@ write_autosubmit_conf <- function(recipe, nchunks) {
     } else if (conf_type == "jobs") {
       # Section 3: jobs
       ## wallclock, notify_on, platform?, processors, 
+      # Create bash file to associate chunk number to recipe name
+      chunk_file <- paste0(proj_dir, "chunk_to_recipe")
+      .create_bash_file(fileout = chunk_file,
+                        dictionary = chunk_to_recipe,
+                        variable = "CHUNK")
+      # Define job parameters
       conf$JOBS$verification$WALLCLOCK <- recipe$Run$auto_conf$wallclock
       if (recipe$Run$auto_conf$notify_completed) {
         conf$JOBS$verification$NOTIFY_ON <- paste(conf$JOBS$verification$NOTIFY_ON,
@@ -52,6 +73,7 @@ write_autosubmit_conf <- function(recipe, nchunks) {
 	  (!recipe$Analysis$Workflow$Scorecards$execute)) {
         conf$JOBS$scorecards <- NULL
       } else {
+        ## TODO: Add multimodel dependency
         if (recipe$Run$auto_conf$notify_completed) {
           conf$JOBS$scorecards$NOTIFY_ON <- paste(conf$JOBS$scorecards$NOTIFY_ON,
                                                   "COMPLETED")
@@ -61,6 +83,41 @@ write_autosubmit_conf <- function(recipe, nchunks) {
                                                   "FAILED")
         }
       }
+      # Only include Multimodel job if sections exists in the recipe
+      # is set to execute = 'True' or 'both'
+      if (!is.null(recipe$Analysis$Datasets$Multimodel) && 
+          tolower(recipe$Analysis$Datasets$Multimodel$execute) == "false") {
+        conf$JOBS$multimodel <- NULL
+      } else {
+        # Create bash file to associate split number to recipe name
+        split_file <- paste0(proj_dir, "split_to_recipe")
+        .create_bash_file(fileout = split_file,
+                          dictionary = split_to_recipe,
+                          variable = "SPLIT")
+        # Define multimodel dependencies in the format required by AS config
+        mm_dependencies <- lapply(split_to_recipe, grep, chunk_to_recipe)
+        mm_dependencies <- lapply(mm_dependencies, paste, collapse = ",")
+        names(mm_dependencies) <- paste(1:length(mm_dependencies))
+        for (split in names(mm_dependencies)) {
+          conf$JOBS$multimodel$DEPENDENCIES$verification$SPLITS_FROM[[split]]$CHUNKS_TO <- 
+            mm_dependencies[[split]]
+        }
+        # 'Splits' parameter should be the number of mulimodel jobs
+        conf$JOBS$multimodel$SPLITS <- length(mm_dependencies)
+        # Define the rest of the parameters
+        if (recipe$Run$auto_conf$notify_completed) {
+          conf$JOBS$multimodel$NOTIFY_ON <- paste(conf$JOBS$multimodel$NOTIFY_ON,
+                                                  "COMPLETED")
+        }
+        if (recipe$Run$auto_conf$notify_failed) {
+          conf$JOBS$multimodel$NOTIFY_ON <- paste(conf$JOBS$multimodel$NOTIFY_ON,
+                                                  "FAILED")
+        }
+
+        conf$JOBS$multimodel$PROCESSORS <- recipe$Run$auto_conf$processors_per_job
+        conf$JOBS$multimodel$CUSTOM_DIRECTIVES <- recipe$Run$auto_conf$custom_directives
+        conf$JOBS$multimodel$WALLCLOCK <- recipe$Run$auto_conf$wallclock
+      }
     } else if (conf_type == "platforms") {
       # Section 4: platform configuration
       ## nord3v2 configuration... platform name? user, processors_per_node
@@ -107,3 +164,15 @@ write_autosubmit_conf <- function(recipe, nchunks) {
     print(paste("nohup autosubmit run", expid, "& disown"))
   }
 }
+
+.create_bash_file <- function(fileout, dictionary, variable) {
+  file_connection <- file(fileout)
+  script_lines <- paste0("case $", variable, " in")
+  for (item in 1:length(dictionary)) {
+    script_command <- paste0("  ", item, ") recipe='", dictionary[[item]], "' ;;")
+    script_lines <- c(script_lines, script_command)
+  }
+  script_lines <- c(script_lines, "esac")
+  writeLines(script_lines, file_connection)
+  close(file_connection)
+}