diff --git a/R/PlotVsLTime.R b/R/PlotVsLTime.R
index dd068adc1dfbcad141d517e7f50a2c30f86493c1..ff7a37c1c7c85de5613616b10731ef88f7965246 100644
--- a/R/PlotVsLTime.R
+++ b/R/PlotVsLTime.R
@@ -31,8 +31,12 @@ PlotVsLTime <- function(var, toptitle = '', ytitle = '', monini = 1, freq = 12,
   nexp <- dim(var)[1]
   nobs <- dim(var)[2]
   if (is.null(limits) == TRUE) {
-    ll <- min(var, na.rm = TRUE)
-    ul <- max(var, na.rm = TRUE)
+	if (all(is.na(var > 0))) {
+      ll <- ul <- 0
+	} else {
+      ll <- min(var, na.rm = TRUE)
+      ul <- max(var, na.rm = TRUE)
+    }
     if (biglab) {
       ul <- ul + 0.4 * (ul - ll)
     } else {
diff --git a/inst/config/BSC.conf b/inst/config/BSC.conf
index 1769e689bf18fd87966a5cc461964726cfe5a484..d9d617eb1a4f505013fd9af4bf7b998e28176015 100644
--- a/inst/config/BSC.conf
+++ b/inst/config/BSC.conf
@@ -22,7 +22,7 @@ DEFAULT_DIM_NAME_MEMBERS = ensemble
 # Helper variables
 EXP_FILE = $VAR_NAME$_*$START_DATE$*.nc
 OBS_FILE = $VAR_NAME$_$YEAR$$MONTH$*.nc
-RECON_LIST = (20thcr_v2|copernicus012|ds083.2|ecco|era40|era40-erainterim|eraclim|erainterim|erainterim-lores|eraland|gecco_v2|gfs|glorys2_v1|glorys2_v3|glosea5|ishii-kimoto|jra55|merra|merra_v2|ncep-reanalysis|oaflux|oi_v2|orap5|piomas|seaice-lim2|sst|tropflux)
+RECON_LIST = (20thcr_v2|copernicus012|ds083.2|ecco|era40|era40-erainterim|eraclim|erainterim|erainterim-lores|eraland|gecco_v2|gfs|glorys2_v1|glorys2_v3|glosea5|ishii-kimoto|jra55|merra|merra_v2|ncep-reanalysis|oaflux|oi_v2|orap5|piomas|seaice-lim2|sst|tropflux|nemovar_system4)
 
 # Defaults for extra variables from Load
 DEFAULT_FILE_GRID = 
diff --git a/inst/doc/plot_maps.R b/inst/doc/plot_maps.R
index c5970008b3e3d69e4a73af77d91ed5c4008c2ed1..0393f833e843a66cc13e2654c248dce663f04f40 100755
--- a/inst/doc/plot_maps.R
+++ b/inst/doc/plot_maps.R
@@ -1,26 +1,39 @@
 #!/usr/bin/env Rscript
 
 library(s2dverification)
+library(ncdf4)
+library(abind)
 args <- commandArgs(TRUE)
 
-comptrend <- T    # Trend as a function of the start date for each leadtime
-compcor <- T      # Correlation Coefficient
-comprms <- T      # Root Mean Square Error
-comprmsss <- T    # Root Mean Square Skill Score
-compratrms <- T   # Ratio RMSE expid1 / expid2
+comptrend <- TRUE    # Trend as a function of the start date for each leadtime
+compcor <- TRUE      # Correlation Coefficient
+comprms <- TRUE      # Root Mean Square Error
+comprmsss <- TRUE    # Root Mean Square Skill Score
+compratrms <- TRUE   # Ratio RMSE expid1 / expid2
 
 var <- args[1]    # tos/tas/pr
-season <- args[2] # Year/DJF/MAM/JJA/SON
-ly <- args[3]     # ly1/ly2-5/ly6-9 for Forecast year 1 / years 2 to 5 / years 
+#season <- args[2] # Year/DJF/MAM/JJA/SON
+season <- 'all'
+#ly <- args[3]     # ly1/ly2-5/ly6-9 for Forecast year 1 / years 2 to 5 / years 
                   # 6 to 9
+ly <- 'all'
 nltimemax <- 124  # number of leadtimes max in the experiments (in months)
+nltimechunk <- 4  # number of leadtimes per chunk (in months)
 lstexpid <- c('i00k','b02p') # list of ids
+grid <- '320x160' # atmospheric grid for tos/pr (ocean/land only)
 mon0 <- 11        # initial month
 year0 <- 1960     # first start date
 yearf <- 2005     # last start date
 intsdate <- 5     # interval between start dates
+runmeanlen <- 12  # length of the window for running mean (in months)
+members <- list(19900101 = c('r1i1p1')) 
+#PUT IN ORDER. NONE CAN BE MISSING
+chunks <- list(19900101 = c('199001-199001', '199002-199002', '199003-199003'))
+
+
 
 obs <- switch(var, 'tas' = 'ghcnersstgiss', 'tos' = 'ersst_v3b', 'pr' = 'gpcc1x1_v6')
+
 syears <- seq(year0, yearf, intsdate)
 imon2 <- paste("0", as.character(mon0), sep = "")
 sdates <- paste(as.character(syears), substr(imon2, nchar(imon2) - 1, 
@@ -30,57 +43,167 @@ savename <- paste(var, '_', season, '_', ly, sep = '')
 for (expid in lstexpid ) {
   savename <- paste(savename, '_', expid, sep = '')
 }
-savename <- paste(savename, '.sav', sep = '')
 
-if (file.exists(savename)) {
-  load(savename)
+if (file.exists(paste(savename, '.sav', sep = ''))) {
+  cat(paste0("Loading existing data from backup file ", savename, "...\n"))
+  load(paste(savename, '.sav', sep = ''))
+  if (is.null(toto1$mod) || is.null(toto1$obs)) {
+    stop("Missing experimental or observational data in backup file. Remove the backup to try retrieving again.")
+  }
 } else {
-  if (is.na(match('b02p', lstexpid)) == TRUE) {
-    lstload <- lstexpid
-  } else {
-    lstload <- lstexpid[-match('b02p', lstexpid)]
+  latbnd <- switch(var, 'tos' = c(-60, 65), 'pr' = c(-60, 65), c(-90, 90))
+  varname <- switch(var, 'siarea' = paste(var, tolower(pole), sep = ''), 
+                         'siextent' = paste(var, tolower(pole), sep = ''), 
+                         'sivol' = paste(var, tolower(pole), sep = ''), 
+                          var)
+  #if (is.na(match('b02p', lstexpid)) == TRUE) {
+    #lstload <- lstexpid
+  #} else {
+    #lstload <- lstexpid[-match('b02p', lstexpid)]
+  #}
+  exp <- list(path = paste0('/es*/exp/ecearth/', lstexpid[1], '/monthly_mean/$VAR_NAME$*/$VAR_NAME$_*_S$START_DATE$_$MEMBER$_$CHUNK$.nc'))
+  b02p <- list(name = 'b02p')
+  lstobs <- obs
+
+
+  toto1 <- list()
+  leading_empty_sdates <- 0
+  greatest_n_of_members <- max(sapply(members, length))
+  greatest_n_of_ltimes <- max(sapply(chunks, length)) * nltimechunk
+  if (any(sapply(chunks, length) == 1)) {
+	greatest_n_of_ltimes <- max(nltimeout, greatest_n_of_ltimes)
   }
-  toto <- Load(var, lstload, obs,sdates, nleadtime = nltimemax,
-               leadtimemin = switch(ly, 'ly1' = 1, 'ly2-5' = 13, 'ly6-9' = 61),
-               leadtimemax = switch(ly, 'ly1' = switch(season, 'SON' = 13, 12),
-               'ly2-5' = switch(season, 'SON' = 61, 60), 
-               'ly6-9' = switch(season, 'SON' = 109, 108)), output = 'lonlat')
-  if (is.na(match('b02p', lstexpid)) == FALSE) {
-    toto1bis <- Load(var, 'b02p', obs = NULL, '19501101', output = 'lonlat')
-    toto1ter <- Histo2Hindcast(toto1bis$mod, '19501101', paste(as.character(
-                syears + switch(ly, 'ly1' = 0, 'ly2-5' = 1, 'ly6-9' = 5)),
-                substr(imon2, nchar(imon2) - 1, nchar(imon2)), '01', sep = ""),
-                nleadtimesout = switch(ly, 'ly1' = switch(season, 'SON' = 13,
-                12), switch(season, 'SON' = 49, 48)))
-    toto1beta <- array(dim = c(dim(toto$mod)[1] + dim(toto1ter)[1], 
-                 max(dim(toto$mod)[2], dim(toto1ter)[2]), dim(toto$mod)[3:6]))
-    toto1beta[1:dim(toto$mod)[1], 1:dim(toto$mod)[2], , , , ] <- toto$mod
-    toto1beta[(dim(toto$mod)[1] + 1):(dim(toto$mod)[1] + dim(toto1ter)[1]),
+  for (sdate in sdates) {
+    sdate_number <- which(sdates == sdate)
+	members_exist <- FALSE
+	chunks_exist <- FALSE
+	if (sdate %in% names(members)) {
+      if (length(members[[sdate]]) > 0) {
+		members_exist <- TRUE
+	  }
+	}
+	if (sdate %in% names(chunks)) {
+	  if (length(chunks[[sdate]]) > 0) {
+		chunks_exist <- TRUE
+	  }
+	}
+	load_exp_data <- FALSE
+    toto1_exp_sdate <- NULL
+	if (members_exist && chunks_exist) {
+      load_exp_data <- TRUE
+	} else {
+      if (is.null(toto1$mod)) {
+		leading_empty_sdates <- leading_empty_sdates + 1
+	  } else {
+		sdate_dims <- dim(toto1$mod)
+		sdate_dims[3] <- 1
+	    toto1_exp_sdate <- array(dim = sdate_dims)
+      }
+	}
+	load_obs_data <- FALSE
+	if (!is.null(lstobs)) {
+	  load_obs_data <- TRUE
+	}
+    ltimes_to_take <- ifelse(length(chunks[[sdate]]) == 1, nltimeout, nltimechunk)
+	toto1_obs_sdate <- NULL
+    for (member in members[[sdate]]) {
+      member_number <- which(members[[sdate]] == member)
+      if (load_obs_data && member_number > 1) {
+        load_obs_data <- FALSE
+      }
+      toto1_exp_member <- NULL
+      for (chunk in chunks[[sdate]]) {
+        chunk_number <- which(chunks[[sdate]] == chunk)
+		if (load_exp_data) {
+          lstload <- list(exp)
+          lstload[[1]][['path']] <- gsub('$MEMBER$', member, lstload[[1]][['path']], fixed = TRUE)
+          lstload[[1]][['path']] <- gsub('$CHUNK$', chunk, lstload[[1]][['path']], fixed = TRUE)
+          toto1_exp_chunk <- Load(varname, lstload, NULL, sdate, leadtimemax = ltimes_to_take, dimnames = list(member = 'region'), output = 'lonlat', grid = grid)
+          if (is.null(toto1_exp_chunk$mod)) {
+            stop("Failed retrieving data for the experiment for the start date ", sdate, ", member ", member, " and chunk ", chunk, ".")
+          }
+          toto1_exp_member <- abind(toto1_exp_member, toto1_exp_chunk$mod, along = 4)
+          if (!('lon' %in% names(toto1))) {
+			toto1[['lon']] <- toto1_exp_chunk$lon
+			toto1[['lat']] <- toto1_exp_chunk$lat
+	      }
+        }
+        if (load_obs_data) {
+          toto1_obs_chunk <- Load(varname, NULL, lstobs, sdate, leadtimemin = ltimes_to_take * (chunk_number - 1) + 1, leadtimemax = ltimes_to_take * chunk_number, output = 'lonlat', grid = grid)
+          if (is.null(toto1_obs_chunk$obs)) {
+            stop("Failed retrieving data for the observation for the start date ", sdate, ", member ", member, " and chunk ", chunk, ".")
+          }
+          toto1_obs_sdate <- abind(toto1_obs_sdate, toto1_obs_chunk$obs, along = 4)
+        }
+      }
+	  if (dim(toto1_exp_member)[4] < greatest_n_of_ltimes) {
+		padding_array_dims <- dim(toto1_exp_member)
+		padding_array_dims[4] <- greatest_n_of_ltimes - padding_array_dims[4]
+		toto1_exp_member <- abind(toto1_exp_member, array(dim = padding_array_dims), along = 4)
+      }
+      toto1_exp_sdate <- abind(toto1_exp_sdate, toto1_exp_member, along = 2)
+	  if (dim(toto1_obs_sdate)[4] < greatest_n_of_ltimes) {
+		padding_array_dims <- dim(toto1_obs_sdate)
+		padding_array_dims[4] <- greatest_n_of_ltimes - padding_array_dims[4]
+		toto1_obs_sdate <- abind(toto1_obs_sdate, array(dim = padding_array_dims), along = 4)
+      }
+	}
+	if (!is.null(toto1_exp_sdate)) {
+      if (dim(toto1_exp_sdate)[2] < greatest_n_of_members) {
+        padding_array_dims <- dim(toto1_exp_sdate)
+	    padding_array_dims[4] <- greatest_n_of_ltimes - padding_array_dims[4]
+	    toto1_exp_sdate <- abind(toto1_exp_sdate, array(dim = padding_array_dims), along = 4)
+      }
+      toto1$mod <- abind(toto1$mod, toto1_exp_sdate, along = 3)
+	}
+	if (!is.null(toto1_obs_sdate)) {
+	  toto1$obs <- abind(toto1$obs, toto1_obs_sdate, along = 3)
+    }
+  }
+  if (leading_empty_sdates > 0) {
+	padding_array_dims <- dim(toto1$mod)
+	padding_array_dims[3] <- leading_empty_sdates
+	toto1$mod <- abind(array(dim = leading_empty_sdates), toto1$mod, along = 3)
+  }
+  if (('b02p' %in% lstexpid)) {
+    toto1bis <- Load(varname, 'b02p', obs = NULL, '19501101', output = 'lonlat', grid = grid)
+    if (is.null(toto1bis$mod)) {
+      stop("Failed retrieving data for b02p.")
+    }
+    toto1ter <- Histo2Hindcast(toto1bis$mod, '19501101', sdates, 
+                               nleadtimesout = ifelse(max(sapply(chunks, length)) == 1, nltimeout, max(sapply(chunks, length)) * nltimechunk))
+    toto1beta <- array(dim = c(dim(toto1$mod)[1] + dim(toto1ter)[1], 
+                       max(dim(toto1$mod)[2], dim(toto1ter)[2]),
+                       dim(toto1$mod)[3:6]))
+    toto1beta[1:dim(toto1$mod)[1], 1:dim(toto1$mod)[2], , , , ] <- toto1$mod
+    toto1beta[(dim(toto1$mod)[1] + 1):(dim(toto1$mod)[1] + dim(toto1ter)[1]), 
               1:dim(toto1ter)[2], , , , ] <- toto1ter
-    toto$mod <- toto1beta
-    lstexpid <- c(lstload, 'b02p')
+    toto1$mod <- toto1beta
   }
-  toto_exp <- InsertDim(Mean1Dim(Season(toto$mod, 4, mon0, switch(season,
-                        'Year' = mon0, 'DJF' = 12, 'MAM' = 3, 'JJA' = 6,
-                        'SON' = 9), switch(season, 
-                        'Year' = (mon0 + 12 - 2) %% 12 + 1, 'DJF' = 2, 
-                        'MAM' = 5, 'JJA' = 8, 'SON' = 11)), 4), 4, 1)
-  toto_obs <- InsertDim(Mean1Dim(Season(toto$obs, 4, mon0, switch(season,
-                        'Year' = mon0, 'DJF' = 12, 'MAM' = 3, 'JJA' = 6,
-                        'SON' = 9), switch(season, 
-                        'Year' = (mon0 + 12 - 2) %% 12 + 1, 'DJF' = 2, 
-                        'MAM' = 5, 'JJA' = 8, 'SON' = 11)), 4), 4, 1)
+
+  toto1$mod <- InsertDim(Mean1Dim(Season(toto1$mod, 4, mon0, mon0, (mon0 + min(11, greatest_n_of_ltimes - 1) - 1) %% 12 + 1), 4), 4, 1)
+  toto1$obs <- InsertDim(Mean1Dim(Season(toto1$obs, 4, mon0, mon0, (mon0 + min(11, greatest_n_of_ltimes - 1) - 1) %% 12 + 1), 4), 4, 1)
+
   if (var == 'pr') {
-    toto_exp <- toto_exp * 1000 * 3600 * 24
-    toto_obs <- toto_obs * 1000 * 3600 * 24
+    toto1$mod <- toto1$mod * 3600 #1000 * 3600 * 24
+    toto1$obs <- toto1$obs * 1000 * 3600 * 24
+  }
+  if (var == 'heatcsum' | var == 'heatcsum1-17' | var == 'heatcsum18-22' | var == 'heatcsum23-46') {
+    toto1$mod <- toto1$mod / 1e22
+    toto1$obs <- toto1$obs / 1e22
+  }
+  if (var == 'siarea' | var=='siextent' | var=='sivol') {
+    toto1$mod <- toto1$mod/1000
+    if (var == 'sivol') {
+      toto1$obs <- toto1$obs/1000
+    }
   }
-  toto=list(mod=toto_exp,obs=toto_obs,lat=toto$lat,lon=toto$lon)
-  save(toto,file=savename)
+  save(toto1, file = paste(savename, '.sav', sep = ''))
 }
 
-clims <- Clim(toto$mod, toto$obs)
-ano_exp <- Ano(toto$mod, clims$clim_exp)
-ano_obs <- Ano(toto$obs, clims$clim_obs)
+clims <- Clim(toto1$mod, toto1$obs)
+ano_exp <- Ano(toto1$mod, clims$clim_exp)
+ano_obs <- Ano(toto1$obs, clims$clim_obs)
 
 if (compcor) {
   cor <- Corr(Mean1Dim(ano_exp, 2), Mean1Dim(ano_obs, 2), 1, 2)
@@ -92,7 +215,7 @@ if (compcor) {
     flag[which(cor[jexp, 1, 2, 1, , ] > cor[jexp, 1, 4, 1, , ])] <- T
     postscript(paste('CorCoef2d_', var, '_', lstexpid[jexp], '_', season, '_',
                ly, '.eps', sep = ''))
-    PlotEquiMap(cor[jexp, 1, 2, 1, , ], toto$lon, toto$lat, 
+    PlotEquiMap(cor[jexp, 1, 2, 1, , ], toto1$lon, toto1$lat, 
                 toptitle = paste('Correlation Coefficient', lstexpid[jexp],
                 switch(season, 'Year' = 'annual', season), switch(var, 
                 'tas' = 'near surface temperature', 
@@ -102,7 +225,7 @@ if (compcor) {
                 'ly1' = 'Year1', 'ly2-5' = 'Year2-5', 'ly6-9' = 'Year6-9')),
                 sizetit = 0.8, brks = lims, cols = cols, colNA = switch(var,
                 'pr' = 'white', grey(0.4)), filled.continents = switch(var,
-                'tas' = F, 'tos' = T, 'pr' = F), dots = t(flag), intylat = 45)
+                'tas' = F, 'tos' = T, 'pr' = F), dots = flag, intylat = 45)
     dev.off()
   }
 }
@@ -117,7 +240,7 @@ if (comprms) {
   for (jexp in 1:length(lstexpid)) {
     postscript(paste('RMSE2d_', var, '_', lstexpid[jexp], '_', season, '_', ly,
                '.eps', sep = ''))
-    PlotEquiMap(rmse[jexp, 1, 2, 1, , ], toto$lon, toto$lat, 
+    PlotEquiMap(rmse[jexp, 1, 2, 1, , ], toto1$lon, toto1$lat, 
                 toptitle = paste('RMSE', lstexpid[jexp], switch(season, 
                 'Year' = 'annual', season), switch(var, 
                 'tas' = 'near surface temperature', 
@@ -143,7 +266,7 @@ if (comprmsss) {
     rmsss[which(-1 > rmsss)] = -1
     postscript(paste('RMSSS2d_', var, '_', lstexpid[jexp], '_', season, '_', ly,
                '.eps', sep = ''))
-    PlotEquiMap(rmsss[jexp, 1, 1, 1, , ], toto$lon, toto$lat, 
+    PlotEquiMap(rmsss[jexp, 1, 1, 1, , ], toto1$lon, toto1$lat, 
                 toptitle = paste('RMSSS', lstexpid[jexp], switch(season,
                 'Year' = 'annual', season), switch(var, 
                 'tas' = 'near surface temperature', 
@@ -153,7 +276,7 @@ if (comprmsss) {
                 'ly1' = 'Year1', 'ly2-5' = 'Year2-5', 'ly6-9' = 'Year6-9')),
                 sizetit = 0.8, brks = lims, cols = cols, colNA = switch(var,
                 'pr' = 'white', grey(0.4)), filled.continents = switch(var,
-                'tas' = F, 'tos' = T, 'pr' = F), dots = t(flag), intylat = 45)
+                'tas' = F, 'tos' = T, 'pr' = F), dots = flag, intylat = 45)
     dev.off()
   }
 }
@@ -169,7 +292,7 @@ if (compratrms) {
   flag[which(ratrms[2, , ] < 0.05)] <- T
   postscript(paste('Rati_RMSE2d_', var, '_', lstexpid[1], '_', lstexpid[2], 
                    '_', season, '_', ly, '.eps', sep = ''))
-  PlotEquiMap(ratrms[1, , ], toto$lon, toto$lat, toptitle = paste('RMSE',
+  PlotEquiMap(ratrms[1, , ], toto1$lon, toto1$lat, toptitle = paste('RMSE',
               lstexpid[1], '/ RMSE', lstexpid[2], switch(season, 
               'Year' = 'annual', season), switch(var, 
               'tas' = 'near surface temperature', 
@@ -179,6 +302,6 @@ if (compratrms) {
               'ly6-9' = 'Year6-9')), sizetit = 0.8, brks = lims, cols = cols,
               colNA = switch(var, 'pr' = 'white', grey(0.4)), 
               filled.continents = switch(var, 'tas' = F, 'tos' = T, 'pr' = F),
-              dots = t(flag), intylat = 45)
+              dots = flag, intylat = 45)
   dev.off()
 }
diff --git a/inst/doc/plot_timeseries.R b/inst/doc/plot_timeseries.R
index 2d36397f26315323a7f9af122bb17630c6807b30..cf92ec0489fd6b663cc1fe24bf8a1c03ef092ca3 100755
--- a/inst/doc/plot_timeseries.R
+++ b/inst/doc/plot_timeseries.R
@@ -2,6 +2,7 @@
 
 library(s2dverification)
 library(ncdf4)
+library(abind)
 args <- commandArgs(TRUE)
 
 comptrend <- TRUE
@@ -10,10 +11,11 @@ comprms <- TRUE
 compspread <- TRUE
 plotano <- TRUE
 
-var <- args[1]    # siextent/siarea/sivol/tos/tas/pr/ohcsum/ohcsum1-17/ohcsum18-22/ohcsum23-46/vsftmyz
+var <- args[1]    # siextent/siarea/sivol/tos/tas/pr/heatcsum/heatcsum1-17/heatcsum18-22/heatcsum23-46/vsftmyz
 pole <- args[2]   # N/S only for siarea/siextent
 nltimemax <- 124  # number of leadtimes max in the experiments (in months)
 nltimeout <- 60   # number of leadtimes to postprocess(in months)
+nltimechunk <- 4  # number of leadtimes per chunk (in months)
 lstexpid <- c('i00k', 'b02p') # list of ids
 grid <- '320x160' # atmospheric grid for tos/pr (ocean/land only)
 vertnem <- 'L42'  # Number of vertical levels in nemo
@@ -22,29 +24,32 @@ year0 <- 1960     # first start date
 yearf <- 2005     # last start date
 intsdate <- 5     # interval between start dates
 runmeanlen <- 12  # length of the window for running mean (in months)
+members <- list(19900101 = c('r1i1p1')) 
+#PUT IN ORDER. NONE CAN BE MISSING
+chunks <- list(19900101 = c('199001-199001', '199002-199002', '199003-199003'))
 
 obs <- switch(var, 'siarea' = c('hadisst_v1.1'), 'siextent' = c('hadisst_v1.1'), 
               'tas' = c('ncep-reanalysis', 'era40'), 'tos' = c('ersst_v3b', 
                                                                'hadisst_v1.1'),
-              'pr' = c('cru_v3.0', 'gpcc1x1_v6'), 'ohcsum' = c('nemovar_system4'),
-              'ohcsum18-22' = c('nemovar_system4'), 'ohcsum23-46' = c('nemovar_system4'),
-              'ohcsum1-17' = c('nemovar_system4'), 'vsftmyz' = c('nemovar_system4'),
+              'pr' = c('cru_v3.0', 'gpcc1x1_v6'), 'heatcsum' = c('nemovar_system4'),
+              'heatcsum18-22' = c('nemovar_system4'), 'heatcsum23-46' = c('nemovar_system4'),
+              'heatcsum1-17' = c('nemovar_system4'), 'vsftmyz' = c('nemovar_system4'),
               'sivol' = 'piomas')
 toptitle2 <- switch(var, 'siarea' = "sea ice area", 'siextent' = "sea ice extent",
                     'sivol' = "sea ice volume", 'tas' = "global T2m", 
                     'tos' = "global SST (60S-65N)", 
                     'pr' = 'land precipitation (60S-65N)', 
-                    'ohcsum' = "global ocean heat content", 
-                    'ohcsum1-17' = 'global 800m-bottom ocean heat content',
-                    'ohcsum18-22' = 'global 350m-800m ocean heat content',
-                    'ohcsum23-46' = 'global 0-350m ocean heat content',
+                    'heatcsum' = "global ocean heat content", 
+                    'heatcsum1-17' = 'global 800m-bottom ocean heat content',
+                    'heatcsum18-22' = 'global 350m-800m ocean heat content',
+                    'heatcsum23-46' = 'global 0-350m ocean heat content',
                     'vsftmyz' = 'Atlantic Overturning Streamfunction (40-55N, 1-2km)'
                     )
 ytitle1 <- switch(var, 'siarea' = "Millions km2", 'siextent' = "Millions km2", 
                   'sivol' = 'Thousands km3', 'tas' = 'K', 'tos' = 'K', 
-                  'pr' = 'mm/day', 'ohcsum' = '10e22 Joules', 
-                  'ohcsum1-17' = '10e22 Joules', 'ohcsum18-22' = '10e22 Joules',
-                  'ohcsum23-46' = '10e22 Joules', 'vsftmyz' = 'Sv')
+                  'pr' = 'mm/day', 'heatcsum' = '10e22 Joules', 
+                  'heatcsum1-17' = '10e22 Joules', 'heatcsum18-22' = '10e22 Joules',
+                  'heatcsum23-46' = '10e22 Joules', 'vsftmyz' = 'Sv')
 
 syears <- seq(year0, yearf, intsdate)
 imon2 <- paste("0", as.character(mon0), sep = "")
@@ -104,19 +109,117 @@ if (file.exists(paste(savename, '.sav', sep = ''))) {
                          'siextent' = paste(var, tolower(pole), sep = ''), 
                          'sivol' = paste(var, tolower(pole), sep = ''), 
                           var)
-  if (is.na(match('b02p', lstexpid)) == TRUE) { 
-    lstload <- lstexpid
-  } else {
-    lstload <- lstexpid[-match('b02p', lstexpid)]
-  }       
-  toto1 <- Load(varname, lstload, obs, sdates, latmin = latbnd[1],
-                latmax = latbnd[2], nleadtime = nltimemax, 
-                leadtimemax = nltimeout, maskmod = lstmask, 
-                maskobs = lstmaskobs)
-  if (is.null(toto1$mod) || is.null(toto1$obs)) {
-    stop("Failed retrieving data for all the experiments or observations.")
+#  if (is.na(match('b02p', lstexpid)) == TRUE) { 
+#    lstload <- lstexpid
+#  } else {
+#    lstload <- lstexpid[-match('b02p', lstexpid)]
+#  }
+  exp <- list(path = paste0('/es*/exp/ecearth/', lstexpid[1], '/monthly_mean/$VAR_NAME$*/$VAR_NAME$_*_S$START_DATE$_$MEMBER$_$CHUNK$.nc'))
+  b02p <- list(name = 'b02p')
+  lstobs <- obs
+
+  toto1 <- list()
+  leading_empty_sdates <- 0
+  greatest_n_of_members <- max(sapply(members, length))
+  greatest_n_of_ltimes <- max(sapply(chunks, length)) * nltimechunk
+  if (any(sapply(chunks, length) == 1)) {
+	greatest_n_of_ltimes <- max(nltimeout, greatest_n_of_ltimes)
+  }
+  for (sdate in sdates) {
+    sdate_number <- which(sdates == sdate)
+	members_exist <- FALSE
+	chunks_exist <- FALSE
+	if (sdate %in% names(members)) {
+      if (length(members[[sdate]]) > 0) {
+		members_exist <- TRUE
+	  }
+	}
+	if (sdate %in% names(chunks)) {
+	  if (length(chunks[[sdate]]) > 0) {
+		chunks_exist <- TRUE
+	  }
+	}
+	load_exp_data <- FALSE
+    toto1_exp_sdate <- NULL
+	if (members_exist && chunks_exist) {
+      load_exp_data <- TRUE
+	} else {
+      if (is.null(toto1$mod)) {
+		leading_empty_sdates <- leading_empty_sdates + 1
+	  } else {
+		sdate_dims <- dim(toto1$mod)
+		sdate_dims[3] <- 1
+	    toto1_exp_sdate <- array(dim = sdate_dims)
+      }
+	}
+	load_obs_data <- FALSE
+	if (!is.null(lstobs)) {
+	  load_obs_data <- TRUE
+	}
+    ltimes_to_take <- ifelse(length(chunks[[sdate]]) == 1, nltimeout, nltimechunk)
+	toto1_obs_sdate <- NULL
+    for (member in members[[sdate]]) {
+      member_number <- which(members[[sdate]] == member)
+      if (load_obs_data && member_number > 1) {
+        load_obs_data <- FALSE
+      }
+      toto1_exp_member <- NULL
+      for (chunk in chunks[[sdate]]) {
+        chunk_number <- which(chunks[[sdate]] == chunk)
+		if (load_exp_data) {
+          lstload <- list(exp)
+          lstload[[1]][['path']] <- gsub('$MEMBER$', member, lstload[[1]][['path']], fixed = TRUE)
+          lstload[[1]][['path']] <- gsub('$CHUNK$', chunk, lstload[[1]][['path']], fixed = TRUE)
+          toto1_exp_chunk <- Load(varname, lstload, NULL, sdate, latmin = latbnd[1],
+                                  latmax = latbnd[2], leadtimemax = ltimes_to_take, 
+                                  maskmod = lstmask, maskobs = lstmaskobs, dimnames = list(member = 'region'))
+          if (is.null(toto1_exp_chunk$mod)) {
+            stop("Failed retrieving data for the experiment for the start date ", sdate, ", member ", member, " and chunk ", chunk, ".")
+          }
+          toto1_exp_member <- abind(toto1_exp_member, toto1_exp_chunk$mod, along = 4)
+        }
+        if (load_obs_data) {
+          toto1_obs_chunk <- Load(varname, NULL, lstobs, sdate, latmin = latbnd[1],
+                                  latmax = latbnd[2], 
+                                  leadtimemin = ltimes_to_take * (chunk_number - 1) + 1,
+                                  leadtimemax = ltimes_to_take * chunk_number, 
+                                  maskmod = lstmask, maskobs = lstmaskobs)
+          if (is.null(toto1_obs_chunk$obs)) {
+            stop("Failed retrieving data for the observation for the start date ", sdate, ", member ", member, " and chunk ", chunk, ".")
+          }
+          toto1_obs_sdate <- abind(toto1_obs_sdate, toto1_obs_chunk$obs, along = 4)
+        }
+      }
+	  if (dim(toto1_exp_member)[4] < greatest_n_of_ltimes) {
+		padding_array_dims <- dim(toto1_exp_member)
+		padding_array_dims[4] <- greatest_n_of_ltimes - padding_array_dims[4]
+		toto1_exp_member <- abind(toto1_exp_member, array(dim = padding_array_dims), along = 4)
+      }
+      toto1_exp_sdate <- abind(toto1_exp_sdate, toto1_exp_member, along = 2)
+	  if (dim(toto1_obs_sdate)[4] < greatest_n_of_ltimes) {
+		padding_array_dims <- dim(toto1_obs_sdate)
+		padding_array_dims[4] <- greatest_n_of_ltimes - padding_array_dims[4]
+		toto1_obs_sdate <- abind(toto1_obs_sdate, array(dim = padding_array_dims), along = 4)
+      }
+	}
+	if (!is.null(toto1_exp_sdate)) {
+      if (dim(toto1_exp_sdate)[2] < greatest_n_of_members) {
+        padding_array_dims <- dim(toto1_exp_sdate)
+	    padding_array_dims[4] <- greatest_n_of_ltimes - padding_array_dims[4]
+	    toto1_exp_sdate <- abind(toto1_exp_sdate, array(dim = padding_array_dims), along = 4)
+      }
+      toto1$mod <- abind(toto1$mod, toto1_exp_sdate, along = 3)
+	}
+	if (!is.null(toto1_obs_sdate)) {
+	  toto1$obs <- abind(toto1$obs, toto1_obs_sdate, along = 3)
+    }
+  }
+  if (leading_empty_sdates > 0) {
+	padding_array_dims <- dim(toto1$mod)
+	padding_array_dims[3] <- leading_empty_sdates
+	toto1$mod <- abind(array(dim = leading_empty_sdates), toto1$mod, along = 3)
   }
-  if (is.na(match('b02p', lstexpid)) == FALSE) {
+  if (('b02p' %in% lstexpid)) {
     toto1bis <- Load(varname, 'b02p', obs = NULL, '19501101', 
                      latmin = latbnd[1], latmax = latbnd[2], maskmod = lstmask,
                      maskobs = lstmaskobs)
@@ -124,7 +227,7 @@ if (file.exists(paste(savename, '.sav', sep = ''))) {
       stop("Failed retrieving data for b02p.")
     }
     toto1ter <- Histo2Hindcast(toto1bis$mod, '19501101', sdates, 
-                               nleadtimesout = nltimeout)
+                               nleadtimesout = ifelse(max(sapply(chunks, length)) == 1, nltimeout, max(sapply(chunks, length)) * nltimechunk))
     toto1beta <- array(dim = c(dim(toto1$mod)[1] + dim(toto1ter)[1], 
                        max(dim(toto1$mod)[2], dim(toto1ter)[2]),
                        dim(toto1$mod)[3:4]))
@@ -132,18 +235,17 @@ if (file.exists(paste(savename, '.sav', sep = ''))) {
     toto1beta[(dim(toto1$mod)[1] + 1):(dim(toto1$mod)[1] + dim(toto1ter)[1]), 
               1:dim(toto1ter)[2], , ] <- toto1ter
     toto1$mod <- toto1beta
-    lstexpid <- c(lstload, 'b02p')
   }
   if (var == 'pr') {
-    toto1$mod <- toto1$mod * 1000 * 3600 * 24
+    toto1$mod <- toto1$mod * 3600 * 24 #1000 * 3600 * 24
     toto1$obs <- toto1$obs * 1000 * 3600 * 24
   }
-  if (var == 'ohcsum' | var == 'ohcsum1-17' | var == 'ohcsum18-22' | var == 'ohcsum23-46') {
+  if (var == 'heatcsum' | var == 'heatcsum1-17' | var == 'heatcsum18-22' | var == 'heatcsum23-46') {
     toto1$mod <- toto1$mod / 1e22
     toto1$obs <- toto1$obs / 1e22
   }
   if (var == 'siarea' | var=='siextent' | var=='sivol') {
-    toto1$mod <- toto1$mod/1000
+    #toto1$mod <- toto1$mod/1000
     if (var == 'sivol') {
       toto1$obs <- toto1$obs/1000
     }