diff --git a/NAMESPACE b/NAMESPACE
index 2e56a9c2..5f6a0f09 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -61,6 +61,7 @@ export(readUSGS)
 export(readindustry_subsectors_specific)
 export(readvanRuijven2016)
 export(readworldsteel)
+export(toolAggregateTimeSteps)
 export(toolBiomassSupplyAggregate)
 export(toolCubicFunctionAggregate)
 export(toolCubicFunctionDisaggregate)
diff --git a/R/calcFEdemand.R b/R/calcFEdemand.R
index 3515fe82..86d28880 100644
--- a/R/calcFEdemand.R
+++ b/R/calcFEdemand.R
@@ -2,7 +2,7 @@
 #'
 #' Returns the Edge data at the Remind level
 #'
-#' @param subtype Final energy (FE) or Energy service (ES) or Useful/Final
+#' @param subtype Final energy (FE) or Useful/Final
 #'   Energy items from EDGEv3 corresponding to REMIND FE items (UE_for_Eff,
 #'   FE_for_Eff)
 #' @param use_ODYM_RECC per-capita pathways for `SDP_xx` scenarios?  (Defaults
@@ -29,7 +29,7 @@
 #' @author Michaja Pehl, Robin Hasse, Falk Benke
 calcFEdemand <- function(subtype, use_ODYM_RECC = FALSE) {
 
-    if (!subtype %in% c("FE", "ES", "FE_for_Eff", "UE_for_Eff")) {
+    if (!subtype %in% c("FE", "FE_for_Eff", "UE_for_Eff")) {
       stop(paste0("Unsupported subtype: ", subtype ))
     }
 
@@ -38,9 +38,9 @@ calcFEdemand <- function(subtype, use_ODYM_RECC = FALSE) {
       getNames(mag, dim = "scenario")
     }
 
-    addDim <- function(x, addnm, dim, dimCode = 3.2) {
-      do.call("mbind", lapply(addnm, function(item) {
-        add_dimension(x, dim = dimCode, add = dim, nm = item)
+    addDimensions <- function(x, dimVals, dimName, dimCode = 3.2) {
+      do.call("mbind", lapply(dimVals, function(item) {
+        add_dimension(x, dim = dimCode, add = dimName, nm = item)
       }))
     }
 
@@ -479,16 +479,9 @@ calcFEdemand <- function(subtype, use_ODYM_RECC = FALSE) {
         }
       }
 
-    } else if (subtype == "ES"){
-      Unit2Million = 1e-6
-
-      services <- readSource("EDGE",subtype="ES_buildings")
-      getSets(services) <- gsub("data", "item", getSets(services))
-      data <- services*Unit2Million
-      unit_out = "million square meters times degree [1e6.m2.C]"
-      description_out = "demand pathways for energy service in buildings"
     }
 
+    # SAME FOR ALL ----
 
     mapping = toolGetMapping(type = "sectoral", name = "structuremappingIO_outputs.csv", where = "mappingfolder")
 
@@ -503,9 +496,6 @@ calcFEdemand <- function(subtype, use_ODYM_RECC = FALSE) {
         mutate(REMINDitems_out = "feelb")
     )
 
-
-    #----- PROCESS DATA ------------------
-
     regions  <- getRegions(data)
     years    <- getYears(data)
     scenarios <- getScens(data)
@@ -567,6 +557,8 @@ calcFEdemand <- function(subtype, use_ODYM_RECC = FALSE) {
     getNames(reminditems) <- gsub("^SSP","gdp_SSP",getNames(reminditems))
     getNames(reminditems) <- gsub("SDP","gdp_SDP",getNames(reminditems))
 
+    # FE only ----
+
     if ('FE' == subtype) {
 
       # ---- _modify SSP1/SSP2 data of CHN/IND further ----
@@ -1881,19 +1873,21 @@ calcFEdemand <- function(subtype, use_ODYM_RECC = FALSE) {
       reminditems <- mbind(
         mselect(reminditems, scenario = nonDuplScenarios),
         mselect(reminditems, scenario = duplScenarios, item = nonIndustryItems),
-        addDim(mselect(reminditems, scenario = "gdp_SSP2EU", item = industryItems,
-                       collapseNames = TRUE),
-               c(paste0("gdp_SSP2EU_NAV_", c("act", "tec", "ele", "lce", "all")),
-                 paste0("gdp_SSP2EU_CAMP_", c("weak", "strong"))),
-               "scenario", 3.1)
+        addDimensions(x = mselect(reminditems, scenario = "gdp_SSP2EU", item = industryItems,
+                                  collapseNames = TRUE),
+                      dimVals = c(paste0("gdp_SSP2EU_NAV_", c("act", "tec", "ele", "lce", "all")),
+                                  paste0("gdp_SSP2EU_CAMP_", c("weak", "strong"))),
+                      dimName = "scenario",
+                      dimCode = 3.1)
       )
     }
 
+    # SAME FOR ALL ----
+
     structure_data <- switch(subtype,
                              FE = "^gdp_(SSP[1-5].*|SDP.*)\\.(fe|ue)",
                              FE_for_Eff = "^gdp_(SSP[1-5]|SDP).*\\.fe.*(b|s)$",
                              UE_for_Eff = "^gdp_(SSP[1-5]|SDP).*\\.fe.*(b|s)$",
-                             ES = "^gdp_(SSP[1-5]|SDP).*\\.esswb$",
                              "^gdp_(SSP[1-5].*|SDP.*)\\.fe..s\\.ue.*b\\.te_ue.*b$")
 
     return(list(x=reminditems,weight=NULL,
diff --git a/R/calcFeDemandBuildings.R b/R/calcFeDemandBuildings.R
index a5b13c55..ca0ff146 100644
--- a/R/calcFeDemandBuildings.R
+++ b/R/calcFeDemandBuildings.R
@@ -116,11 +116,13 @@ calcFeDemandBuildings <- function(subtype) {
   # change item names back from UE to FE
   if (subtype == "UE") {
     getItems(remind, "item") <- gsub("^ue", "fe", getItems(remind, "item"))
-    description <- "useful energy demand in buildings"
-  } else {
-    description <- "demand pathways for final energy in buildings and industry in the original file"
   }
 
+  description <- switch(subtype,
+    FE = "demand pathways for final energy in buildings and industry in the original file",
+    UE = "useful energy demand in buildings"
+  )
+
   outputStructure <- switch(subtype,
     FE = "^gdp_(SSP[1-5]|SDP).*\\..*\\.fe.*b$",
     UE = "^gdp_(SSP[1-5]|SDP).*\\..*\\.fe.*b$"
diff --git a/R/convertEDGE.R b/R/convertEDGE.R
index edba3369..cdc6a3bc 100644
--- a/R/convertEDGE.R
+++ b/R/convertEDGE.R
@@ -11,7 +11,6 @@
 #' mbind
 
 convertEDGE <- function(x, subtype = "FE_stationary") {
-
   #---- Functions -------------
   noYearDim <- function(x) setYears(x, NULL)
 
@@ -65,19 +64,18 @@ convertEDGE <- function(x, subtype = "FE_stationary") {
 
 
   #---- Parameters and Mappings ------
-  rem_years_hist <- seq(1990,2150,5)
-  keep_years <- getYears(x)
+  rem_years_hist <- seq(1990, 2150, 5)
 
-  struct_mapping_path = toolGetMapping(type = "sectoral", name = "structuremappingIO_outputs.csv",
+  struct_mapping_path <- toolGetMapping(type = "sectoral", name = "structuremappingIO_outputs.csv",
                                       returnPathOnly = TRUE, where = "mappingfolder")
-  struct_mapping = read.csv2(struct_mapping_path, na.strings = "")
+  struct_mapping <- read.csv2(struct_mapping_path, na.strings = "")
 
-  #Select the relevant part of the mapping
-  struct_mapping = struct_mapping[!is.na(struct_mapping$weight_convertEDGE),]
-  struct_mapping = unique(struct_mapping[c( "weight_convertEDGE", "EDGEitems")])
+  # Select the relevant part of the mapping
+  struct_mapping <- struct_mapping[!is.na(struct_mapping$weight_convertEDGE), ]
+  struct_mapping <- unique(struct_mapping[c("weight_convertEDGE", "EDGEitems")])
 
 
-  if (subtype %in% c("FE_stationary","FE_buildings")){
+  if (subtype %in% c("FE_stationary", "FE_buildings")) {
     #---- Explanations
     # For the historical data, weights are directly taken from the IEA
     # to ensure the consistency at the country level
@@ -85,21 +83,21 @@ convertEDGE <- function(x, subtype = "FE_stationary") {
     # multiplied by the growth rate of the country
 
     # Load the regional mapping which depends upon the model used
-    if (subtype == "FE_stationary"){
+    if (subtype == "FE_stationary") {
 
       mappingfile <- toolGetMapping(type = "regional", name = "regionmappingREMIND.csv",
                                     returnPathOnly = TRUE, where = "mappingfolder")
       mapping <- read.csv2(mappingfile)
-      region_col = which(names(mapping) == "RegionCode")
-      iso_col = which(names(mapping) == "CountryCode")
+      region_col <- which(names(mapping) == "RegionCode")
+      iso_col <- which(names(mapping) == "CountryCode")
 
-    }else if (subtype %in% c("FE_buildings")){
+    } else if (subtype %in% c("FE_buildings")) {
 
       mappingfile <- toolGetMapping(type = "regional", name = "regionmappingEDGE.csv",
                                     returnPathOnly = TRUE, where = "mappingfolder")
       mapping <- read.csv2(mappingfile)
-      region_col = which(names(mapping) == "RegionCodeEUR_ETP")
-      iso_col = which(names(mapping) == "CountryCode")
+      region_col <- which(names(mapping) == "RegionCodeEUR_ETP")
+      iso_col <- which(names(mapping) == "CountryCode")
 
     }
 
@@ -112,9 +110,9 @@ convertEDGE <- function(x, subtype = "FE_stationary") {
 
 
     #--- Then load the final energy data
-    hist_fe_stationary = calcOutput("IOEdgeBuildings", subtype = "output_EDGE", aggregate = F)
-    hist_fe_buildings = calcOutput("IOEdgeBuildings", subtype = "output_EDGE_buildings", aggregate = F)
-    hist_fe_transport = calcOutput("IO", subtype="output", aggregate = F)
+    hist_fe_stationary <- calcOutput("IOEdgeBuildings", subtype = "output_EDGE", aggregate = FALSE)
+    hist_fe_buildings <- calcOutput("IOEdgeBuildings", subtype = "output_EDGE_buildings", aggregate = FALSE)
+    hist_fe_transport <- calcOutput("IO", subtype = "output", aggregate = FALSE)
 
     wfe <- mbind(hist_fe_stationary, hist_fe_buildings)
 
@@ -124,58 +122,58 @@ convertEDGE <- function(x, subtype = "FE_stationary") {
     # Replace NAs
     x[is.na(x)] <- 0
 
-    if (any(wfe < 0 )){
+    if (any(wfe < 0)) {
       warning("calcOutput('IOEdgeBuildings', subtype = X), with X in (output_EDGE, output_EDGE_buildings) produces negative values, set to 0")
-      wfe[wfe < 0 ] = 0
+      wfe[wfe < 0] <- 0
     }
 
-    #Select last year of X available in the historical data set
-    maxYear_X_in_FE = max(getYears(x, as.integer = T)[getYears(x, as.integer = T) %in% getYears(wfe, as.integer = T)])
+    # Select last year of X available in the historical data set
+    maxYear_X_in_FE <- max(getYears(x, as.integer = TRUE)[getYears(x, as.integer = TRUE) %in% getYears(wfe, as.integer = TRUE)])
     # Deduce the scenario periods
-    exceeding_years = getYears(x, as.integer = T)[getYears(x, as.integer = T) > maxYear_X_in_FE]
+    exceeding_years <- getYears(x, as.integer = TRUE)[getYears(x, as.integer = TRUE) > maxYear_X_in_FE]
 
 
-    #FE_stationary projections are not updated. Therefore, we correct here for the newly published past data
-    #For historical years, the data is substituted. For projections years, there is first a transition period,
-    #before the FE_stationary projections are fully taken up
-    if(subtype == "FE_stationary"){
-      fe_stationary = time_interpolate(hist_fe_stationary[,getYears(hist_fe_stationary)[getYears(hist_fe_stationary,T) <= maxYear_X_in_FE],], # The years exceeding maxYear might not be meaningful. Therefore we exclude them
-                              interpolated_year = c(maxYear_X_in_FE,exceeding_years),
-                              integrate_interpolated_years = T,
+    # FE_stationary projections are not updated. Therefore, we correct here for the newly published past data
+    # For historical years, the data is substituted. For projections years, there is first a transition period,
+    # before the FE_stationary projections are fully taken up
+    if (subtype == "FE_stationary") {
+      fe_stationary <- time_interpolate(hist_fe_stationary[, getYears(hist_fe_stationary)[getYears(hist_fe_stationary, TRUE) <= maxYear_X_in_FE], ], # The years exceeding maxYear might not be meaningful. Therefore we exclude them
+                              interpolated_year = c(maxYear_X_in_FE, exceeding_years),
+                              integrate_interpolated_years = TRUE,
                               extrapolation_type = "constant")
-      fe_stationary = addSSPnames(fe_stationary)
+      fe_stationary <- addSSPnames(fe_stationary)
 
-      fe_transport = time_interpolate(hist_fe_transport[,getYears(hist_fe_transport)[getYears(hist_fe_transport,T) <= maxYear_X_in_FE],], # The years exceeding maxYear might not be meaningful. Therefore we exclude them
-                                       interpolated_year = c(maxYear_X_in_FE,exceeding_years),
-                                       integrate_interpolated_years = T,
+      fe_transport <- time_interpolate(hist_fe_transport[, getYears(hist_fe_transport)[getYears(hist_fe_transport, TRUE) <= maxYear_X_in_FE], ], # The years exceeding maxYear might not be meaningful. Therefore we exclude them
+                                       interpolated_year = c(maxYear_X_in_FE, exceeding_years),
+                                       integrate_interpolated_years = TRUE,
                                        extrapolation_type = "constant")
-      fe_transport = addSSPnames(fe_transport)
-      getSets(fe_transport) = c("region","period","scenario","input","output","tech")
-
-      #change the regional resolution of fe_stationary to match the EDGE_stationary resolution
-      #iso_col and region_col are originally designed for the weights, that is why names are confusing here
-      fe_stationary = toolAggregate(fe_stationary,mappingfile, from = iso_col, to = region_col)
-      fe_transport = toolAggregate(fe_transport,mappingfile, from = iso_col, to = region_col)
-
-      #Item names differ slightly for the input of EDGE_stationary (fe_stationary) and the output
-      #The main issue concerns transport. We therefore restrict to the variables of interest in each data set of
-      #historical data
-      stationary_items = grep("^(fenon|feagr|feind|feoth)", getNames(x,T)[[2]], value = T) # Stationary, non-buildings names
-      transport_items = grep("^(fepet|fedie|feelt)", getNames(x,T)[[2]], value = T) #Transport names
-
-      #simplify transport
-      fe_transport = dimSums(fe_transport,dim = c("input","tech"))
-
-      #create lambda vector that gives 0 to the historical data and 1 after 2030
-      lambda =  calcLambda(exceeding_years,2030,getYears(x)[getYears(x,T) <= maxYear_X_in_FE ])
-      #Replace
-      x[,,stationary_items] = fe_stationary[,getYears(x),stationary_items] * (1 - lambda) + x[,,stationary_items] * lambda
-      x[,,transport_items] = fe_transport[,getYears(x),transport_items] * (1 - lambda) + x[,,transport_items] * lambda
+      fe_transport <- addSSPnames(fe_transport)
+      getSets(fe_transport) <- c("region", "period", "scenario", "input", "output", "tech")
+
+      # change the regional resolution of fe_stationary to match the EDGE_stationary resolution
+      # iso_col and region_col are originally designed for the weights, that is why names are confusing here
+      fe_stationary <- toolAggregate(fe_stationary, mappingfile, from = iso_col, to = region_col)
+      fe_transport <- toolAggregate(fe_transport, mappingfile, from = iso_col, to = region_col)
+
+      # Item names differ slightly for the input of EDGE_stationary (fe_stationary) and the output
+      # The main issue concerns transport. We therefore restrict to the variables of interest in each data set of
+      # historical data
+      stationary_items <- grep("^(fenon|feagr|feind|feoth)", getNames(x, TRUE)[[2]], value = TRUE) # Stationary, non-buildings names
+      transport_items <- grep("^(fepet|fedie|feelt)", getNames(x, TRUE)[[2]], value = TRUE) # Transport names
+
+      # simplify transport
+      fe_transport <- dimSums(fe_transport, dim = c("input", "tech"))
+
+      # create lambda vector that gives 0 to the historical data and 1 after 2030
+      lambda <-  calcLambda(exceeding_years, 2030, getYears(x)[getYears(x, TRUE) <= maxYear_X_in_FE])
+      # Replace
+      x[, , stationary_items] <- fe_stationary[, getYears(x), stationary_items] * (1 - lambda) + x[, , stationary_items] * lambda
+      x[, , transport_items] <- fe_transport[, getYears(x), transport_items] * (1 - lambda) + x[, , transport_items] * lambda
     }
 
     # Scale GDP and FE weights so that they can be added
-    wg = wg/dimSums(wg, dim = 1, na.rm = T)
-    wfe = wfe/dimSums(wfe, dim =1, na.rm = T)
+    wg <- wg / dimSums(wg, dim = 1, na.rm = TRUE)
+    wfe <- wfe / dimSums(wfe, dim = 1, na.rm = TRUE)
 
     # Add some corrections
     wg[is.na(wg)] <- 0
@@ -183,53 +181,53 @@ convertEDGE <- function(x, subtype = "FE_stationary") {
 
     # Add some corrections to the FE data set + add the scenario dimension
     wfe[is.na(wfe)] <- 0
-    wfe = addSSPnames(wfe)
+    wfe <- addSSPnames(wfe)
 
     # Compute lambda
-    lambda = calcLambda(exceeding_years, 2060)
+    lambda <- calcLambda(exceeding_years, 2060)
     # For the future periods, the weight will be a linear combination of last FE weight and of the GDP size.
     # until maxYear_X_in_FE this will be exclusively FE, in 2060 (depending on the threshold value above), exclusively GDP
 
     wfe <- mbind(wfe,
-      lambda[,exceeding_years, ] * wg[,exceeding_years, ] +
+      lambda[, exceeding_years, ] * wg[, exceeding_years, ] +
         (1 - lambda[, exceeding_years, ]) * (noYearDim(wfe[, maxYear_X_in_FE, ]))
     )
 
     # In cases where the variables in EDGE do not exist in the mapping for computing the final energy,
-    #e.g. when EDGE produces further disaggregations, or when it gives REMIND items without computing them
-    wfe = mbind(wfe, renameExtraWeights(x,wfe, struct_mapping))
+    # e.g. when EDGE produces further disaggregations, or when it gives REMIND items without computing them
+    wfe <- mbind(wfe, renameExtraWeights(x, wfe, struct_mapping))
 
-    #Reduce the dimensions of the weights
-    wfe = wfe[,getYears(x), getNames(x, dim = "item")]
+    # Reduce the dimensions of the weights
+    wfe <- wfe[, getYears(x), getNames(x, dim = "item")]
 
-    #Disaggregate and fill the gaps
-    xadd <- toolAggregate(x,mappingfile,weight=wfe,
+    # Disaggregate and fill the gaps
+    xadd <- toolAggregate(x, mappingfile, weight = wfe,
                           from = region_col,
                           to = iso_col)
-    result <- toolCountryFill(xadd,0)
+    result <- toolCountryFill(xadd, 0)
 
-    if(subtype == "FE_stationary"){
+    if (subtype == "FE_stationary") {
       # re-calculating fepet and fedie final energy based on updated EDGE shares
-      share <- readSource(type="EDGETransport", subtype = "shares_LDV_transport")
+      share <- readSource(type = "EDGETransport", subtype = "shares_LDV_transport")
       # for EU regions use JRC data instead
-      JRC_reg <- c("MLT","EST","CYP","LVA","LTU","LUX","SVK","SVN","HRV","BGR","HUN","ROU","FIN","DNK","IRL","CZE","GRC","AUT","PRT","SWE","BEL","NLD","POL","ESP","ITA","GBR","FRA","DEU")
-      JRC <- calcOutput("JRC_IDEES", subtype="Transport", aggregate = FALSE)
-      JRC_share <- new.magpie(JRC_reg,getYears(share),getNames(share),fill=0)
+      JRC_reg <- c("MLT", "EST", "CYP", "LVA", "LTU", "LUX", "SVK", "SVN", "HRV", "BGR", "HUN", "ROU", "FIN", "DNK", "IRL", "CZE", "GRC", "AUT", "PRT", "SWE", "BEL", "NLD", "POL", "ESP", "ITA", "GBR", "FRA", "DEU")
+      JRC <- calcOutput("JRC_IDEES", subtype = "Transport", aggregate = FALSE)
+      JRC_share <- new.magpie(JRC_reg, getYears(share), getNames(share), fill = 0)
       # for years lower or equal to 2015 assume bunkers equal to JRC historical values
-      y1 <- getYears(JRC)[getYears(JRC, as.integer = TRUE)<=2015]
-      JRC_share[JRC_reg,y1,] <- JRC[JRC_reg,y1,"FE|Transport|LDV|Liquids (EJ/yr)"]/(JRC[JRC_reg,y1,"FE|Transport|non-LDV|Liquids (EJ/yr)"]+JRC[JRC_reg,y1,"FE|Transport|LDV|Liquids (EJ/yr)"])
+      y1 <- getYears(JRC)[getYears(JRC, as.integer = TRUE) <= 2015]
+      JRC_share[JRC_reg, y1, ] <- JRC[JRC_reg, y1, "FE|Transport|LDV|Liquids (EJ/yr)"] / (JRC[JRC_reg, y1, "FE|Transport|non-LDV|Liquids (EJ/yr)"] + JRC[JRC_reg, y1, "FE|Transport|LDV|Liquids (EJ/yr)"])
       # for years after 2015 assume LDV share constant and eqaul to JRC 2015 values
-      y2 <- getYears(share)[getYears(share, as.integer = TRUE)>2015]
-      JRC_share[,y2,] <- JRC_share[,2015,]
+      y2 <- getYears(share)[getYears(share, as.integer = TRUE) > 2015]
+      JRC_share[, y2, ] <- JRC_share[, 2015, ]
       ## setting EU shares equal to JRC values
-      varname_SSP2 <- getNames(share[,, "gdp_SSP2"])[1]
-      share[JRC_reg,getYears(JRC_share), varname_SSP2] <- JRC_share[JRC_reg,getYears(JRC_share),]
+      varname_SSP2 <- getNames(share[, , "gdp_SSP2"])[1]
+      share[JRC_reg, getYears(JRC_share), varname_SSP2] <- JRC_share[JRC_reg, getYears(JRC_share), ]
       # redefining LDV and non-LDV liquids
-      feTotal <- dimSums(result[,,c("fepet","fedie")],dim=3.2)
-      feShares <- new.magpie(cells_and_regions = getRegions(share), years = intersect(getYears(share),getYears(result)), names = getNames(result[,,c("fepet","fedie")]))
-      feShares[,,"fepet"] <- setNames(setNames(share[getRegions(share),getYears(feShares),"share_LDV_totliq"],"fepet"),NULL)
-      feShares[,,"fedie"] <- (1-setNames(setNames(share[getRegions(share),getYears(feShares),"share_LDV_totliq"],"fepet"),NULL))
-      feTransp <- new.magpie(cells_and_regions = getRegions(share), years = getYears(feShares), names = getNames(result[,,c("fepet","fedie")]))
+      feTotal <- dimSums(result[, , c("fepet", "fedie")], dim = 3.2)
+      feShares <- new.magpie(cells_and_regions = getRegions(share), years = intersect(getYears(share), getYears(result)), names = getNames(result[, , c("fepet", "fedie")]))
+      feShares[, , "fepet"] <- setNames(setNames(share[getRegions(share), getYears(feShares), "share_LDV_totliq"], "fepet"), NULL)
+      feShares[, , "fedie"] <- (1 - setNames(setNames(share[getRegions(share), getYears(feShares), "share_LDV_totliq"], "fepet"), NULL))
+      feTransp <- new.magpie(cells_and_regions = getRegions(share), years = getYears(feShares), names = getNames(result[, , c("fepet", "fedie")]))
 
       for (i in c("SSP1", "SSP2", "SSP3", "SSP4", "SSP5", "SDP", "SDP_EI", "SDP_RC", "SDP_MC", "SSP2EU")) {
         i1 <- paste0(i, ".fepet")
@@ -239,25 +237,25 @@ convertEDGE <- function(x, subtype = "FE_stationary") {
       }
 
       # extrapolating missing historical years
-      result[,getYears(feTransp),getNames(feTransp)] <- feTransp[,getYears(feTransp),getNames(feTransp)]
+      result[, getYears(feTransp), getNames(feTransp)] <- feTransp[, getYears(feTransp), getNames(feTransp)]
     }
 
-    if(subtype == "FE_buildings"){
+    if (subtype == "FE_buildings") {
       # Attribute the growth in water heating demand of the EDGE Region OCD to TUR,
       # and retrieve it from AUS, CAN, CHE (Swiss), NOR, NZL
       # For SSP1, SSP2 and SDP
-      names_2_change = grep("(SSP1|SSP2|SDP|SDP_EI|SDP_RC|SDP_MC|SSP2EU).*water_heating", getNames(result), value = TRUE)
-      names_2_change_elec = grep("elec",names_2_change,value = T)
-      names_2_change_nonelec = grep("elec",names_2_change,value = T, invert = T)
-      regs_OCD = c("AUS","CAN","CHE","NOR","NZL")
-      reg_TUR = "TUR"
-      end_of_history = 2015
-      scenario_time = getYears(result, T)[getYears(result, T) > end_of_history]
-
-      WH_growth = result[regs_OCD,scenario_time,names_2_change] -  dimReduce(result[regs_OCD,end_of_history,names_2_change])
-      WH_growth[,,names_2_change_elec] = WH_growth[,,names_2_change_elec] * 0.5
+      names_2_change <- grep("(SSP1|SSP2|SDP|SDP_EI|SDP_RC|SDP_MC|SSP2EU).*water_heating", getNames(result), value = TRUE)
+      names_2_change_elec <- grep("elec", names_2_change, value = TRUE)
+      names_2_change_nonelec <- grep("elec", names_2_change, value = TRUE, invert = TRUE)
+      regs_OCD <- c("AUS", "CAN", "CHE", "NOR", "NZL")
+      reg_TUR <- "TUR"
+      end_of_history <- 2015
+      scenario_time <- getYears(result, TRUE)[getYears(result, TRUE) > end_of_history]
+
+      WH_growth <- result[regs_OCD, scenario_time, names_2_change] -  dimReduce(result[regs_OCD, end_of_history, names_2_change])
+      WH_growth[, , names_2_change_elec] <- WH_growth[, , names_2_change_elec] * 0.5
       WH_growth[WH_growth < 0] <- 0
-      WH_growth_agg = dimSums(WH_growth, dim = 1)
+      WH_growth_agg <- dimSums(WH_growth, dim = 1)
 
       result[getRegions(WH_growth), getYears(WH_growth), getNames(WH_growth)] <- result[getRegions(WH_growth), getYears(WH_growth), getNames(WH_growth)] - WH_growth
       result[reg_TUR, getYears(WH_growth), getNames(WH_growth)] <- result[reg_TUR, getYears(WH_growth), getNames(WH_growth)] + WH_growth_agg
@@ -266,86 +264,68 @@ convertEDGE <- function(x, subtype = "FE_stationary") {
 
 
 
-  } else if(subtype %in% c("Capital")){
+  } else if (subtype %in% c("Capital")) {
 
     mappingfile <- toolGetMapping(type = "regional", name = "regionmappingEDGE.csv",
                                   returnPathOnly = TRUE, where = "mappingfolder")
     mapping <- read.csv2(mappingfile)
-    region_col = which(names(mapping) == "RegionCodeEUR_ETP")
-    iso_col = which(names(mapping) == "CountryCode")
+    region_col <- which(names(mapping) == "RegionCodeEUR_ETP")
+    iso_col <- which(names(mapping) == "CountryCode")
 
-    x <- x[,getYears(x,T)[which(getYears(x,T) <= 2100)],]
+    x <- x[, getYears(x, TRUE)[which(getYears(x, TRUE) <= 2100)], ]
     getItems(x, 3.1) <- paste0("gdp_", getItems(x, 3.1))
 
-    wg     <- calcOutput("GDP", aggregate=F)
-    wfe    <- calcOutput("FEdemand", subtype = "FE", aggregate = F)
+    wg     <- calcOutput("GDP", aggregate = FALSE)
+    wfe    <- calcOutput("FEdemand", subtype = "FE", aggregate = FALSE)
 
-    getSets(wg) = gsub("variable","scenario",getSets(wg))
-    getSets(wfe) = gsub("item","data", getSets(wfe))
-    wg = add_dimension(wg,dim = 3.2, add = "data",nm ="kap")
+    getSets(wg) <- gsub("variable", "scenario", getSets(wg))
+    getSets(wfe) <- gsub("item", "data", getSets(wfe))
+    wg <- add_dimension(wg, dim = 3.2, add = "data", nm = "kap")
 
-    #***Reproduce this in the aggregation of CapitalUnit in calcCapital
-    corres_ener_cap = c(kapal = "fealelb",
+    #*** Reproduce this in the aggregation of CapitalUnit in calcCapital
+    corres_ener_cap <- c(kapal = "fealelb",
                         kapsc = "fescelb",
                         kaphc = "ueswb")
-    wfe = do.call(mbind,
-                  lapply(names(corres_ener_cap), function(kap_nm){
-                    ener_nm = corres_ener_cap[kap_nm]
-                    tmp = wfe[,,ener_nm]
-                    getNames(tmp) = gsub(ener_nm,kap_nm,getNames(tmp))
+    wfe <- do.call(mbind,
+                  lapply(names(corres_ener_cap), function(kap_nm) {
+                    ener_nm <- corres_ener_cap[kap_nm]
+                    tmp <- wfe[, , ener_nm]
+                    getNames(tmp) <- gsub(ener_nm, kap_nm, getNames(tmp))
                     return(tmp)
                   })
     )
 
-    years_select = intersect(  intersect(getYears(x),getYears(wg))
-                               ,getYears(wfe))
+    years_select <- intersect(intersect(getYears(x), getYears(wg)),
+getYears(wfe))
 
-    wfe = wfe[,years_select,]
-    wg = wg[,years_select,]
+    wfe <- wfe[, years_select, ]
+    wg <- wg[, years_select, ]
 
-    weights = mbind(wfe,wg)
+    weights <- mbind(wfe, wg)
 
-    x = toolAggregate(x[,years_select,],mappingfile, weight = weights[,,getNames(x)], from = region_col, to = iso_col )
-    result = x
+    x <- toolAggregate(x[, years_select, ], mappingfile, weight = weights[, , getNames(x)], from = region_col, to = iso_col)
+    result <- x
 
-  } else if(subtype %in% c("CapitalUnit")){
+  } else if (subtype %in% c("CapitalUnit")) {
 
     mappingfile <- toolGetMapping(type = "regional", name = "regionmappingEDGE.csv",
                                   returnPathOnly = TRUE, where = "mappingfolder")
     mapping <- read.csv2(mappingfile)
-    region_col = which(names(mapping) == "RegionCodeEUR_ETP")
-    iso_col = which(names(mapping) == "CountryCode")
+    region_col <- which(names(mapping) == "RegionCodeEUR_ETP")
+    iso_col <- which(names(mapping) == "CountryCode")
 
     wg     <- NULL
 
 
-    x = toolAggregate(x[,,],mappingfile, weight = wg, from = region_col, to = iso_col )
-    result = x
-
-  } else if(subtype %in% c("ES_buildings")){
-
-    mappingfile <- toolGetMapping(type = "regional", name = "regionmappingEDGE.csv",
-                                  returnPathOnly = TRUE, where = "mappingfolder")
-    mapping <- read.csv2(mappingfile)
-    region_col = which(names(mapping) == "RegionCodeEUR_ETP")
-    iso_col = which(names(mapping) == "CountryCode")
-
-    select_years = intersect(getYears(x,as.integer = T),rem_years_hist)
-    wg <- calcOutput("GDP", years = select_years, aggregate = FALSE)
-
-    # duplicate SSP2 for SSP2_lowEn an SSP2EU for Navigate and Campaigners scenarios
-    wg <- duplScens(wg)
-    getNames(wg) = gsub("gdp_","", getNames(wg))
-
-    x = toolAggregate(x[,select_years,],mappingfile, weight = wg[,,getNames(x,dim=1)], from = region_col, to = iso_col )
-    result = x
+    x <- toolAggregate(x[, , ], mappingfile, weight = wg, from = region_col, to = iso_col)
+    result <- x
 
   } else if (subtype == "Floorspace") {
     mappingfile <- toolGetMapping(type = "regional", name = "regionmappingEDGE.csv",
                                   returnPathOnly = TRUE, where = "mappingfolder")
     mapping <- read.csv2(mappingfile)
-    region_col = which(names(mapping) == "RegionCodeEUR_ETP")
-    iso_col = which(names(mapping) == "CountryCode")
+    region_col <- which(names(mapping) == "RegionCodeEUR_ETP")
+    iso_col <- which(names(mapping) == "CountryCode")
 
     getNames(x) <- paste0("gdp_", getNames(x))
     wp <- calcOutput("Population", years = rem_years_hist, aggregate = FALSE)
@@ -355,7 +335,7 @@ convertEDGE <- function(x, subtype = "FE_stationary") {
     wp <- duplScens(wp)
 
     x <- toolAggregate(x[, rem_years_hist, ], mappingfile, weight = wp,
-                       from = region_col, to = iso_col )
+                       from = region_col, to = iso_col)
     result <- x
   }
   return(result)
diff --git a/R/fullREMIND.R b/R/fullREMIND.R
index efa8d961..32c242de 100644
--- a/R/fullREMIND.R
+++ b/R/fullREMIND.R
@@ -33,8 +33,8 @@ fullREMIND <- function() {
   calcOutput("TaxLimits", subtype = "propFeSubsidy",  round = 2,  file = "f21_prop_fe_sub.cs4r")
   calcOutput("PETaxes", subtype = "subsidies",        round = 2,  file = "f21_tau_pe_sub.cs4r")
   calcOutput("TaxXport",                              round = 2,  file = "p21_tau_xpres_tax.cs4r")   # not default, overwritten with 0
-  calcOutput("Capital",                               signif = 4,  file = "f29_capitalQuantity.cs4r")
-  calcOutput("Capital",   subtype = "CapitalUnit",    round = 6,  file = "f29_capitalUnitProjections.cs4r")
+  calcOutput("Capital",   subtype = "Capital",        signif = 4, file = "f29_capitalQuantity.cs4r") # TODO
+  calcOutput("Capital",   subtype = "CapitalUnit",    round = 6,  file = "f29_capitalUnitProjections.cs4r") # TODO
   calcOutput("FEdemand",  subtype = "FE",             signif = 4,  file = "f_fedemand.cs4r")
   calcOutput("ExogDemScen",                           round = 8,  file = "p47_exogDemScen.cs4r") # exogenous demand scenarios activated by cm_exogDem_scen
   calcOutput(
@@ -47,7 +47,6 @@ fullREMIND <- function() {
       madrat_mule())
   calcOutput("FeDemandBuildings",  subtype = "FE",   round = 8,  file = "f_fedemand_build.cs4r")
   calcOutput("FeDemandBuildings",  subtype = "UE",   round = 8,  file = "f36_uedemand_build.cs4r")
-  calcOutput("FEdemand",  subtype = "ES",             round = 6,  file = "f29_esdemand.cs4r")
   calcOutput("ChemicalFeedstocksShare", round = 2, file = "p37_chemicals_feedstock_share.cs4r")
   calcOutput("EnergyEffPaths",                        round = 6,  file = "f29_efficiency_growth.cs4r")
   calcOutput("Floorspace", onlyTotal = TRUE,          round = 1,  file = "p36_floorspace_scen.cs4r")
diff --git a/R/readEDGE.R b/R/readEDGE.R
index 1c17f9ae..296ed0c3 100644
--- a/R/readEDGE.R
+++ b/R/readEDGE.R
@@ -12,8 +12,7 @@
 #' @author Antoine Levesque, Robin Hasse
 #' @seealso \code{\link{readSource}}
 #' @importFrom magclass read.magpie mselect as.magpie mbind add_dimension
-readEDGE <- function(subtype = c("FE_stationary", "FE_buildings", "Capital", "CapitalUnit", "Floorspace",
-                                 "ES_buildings")) {
+readEDGE <- function(subtype = c("FE_stationary", "FE_buildings", "Capital", "CapitalUnit", "Floorspace")) {
   subtype <- match.arg(subtype)
 
   # input data version
@@ -41,19 +40,22 @@ readEDGE <- function(subtype = c("FE_stationary", "FE_buildings", "Capital", "Ca
         addDim(mselect(data, scenario = "SSP2", collapseNames = TRUE),
                scenarios$SSP2s, "scenario", 3.1),
         addDim(mselect(data, scenario = "SSP1", collapseNames = TRUE),
-               scenarios$SDPs, "scenario", 3.1))},
+               scenarios$SDPs, "scenario", 3.1))
+    },
     FE_buildings = {
       data <- read.csv(file.path(ver, "EDGE_buildings_energy.csv"))
       data <- as.magpie(data)
       getNames(data) <- gsub("rcp", "", getNames(data))
       getNames(data) <- gsub("NoC", "fixed", getNames(data))
-      getSets(data) <- c("region", "year", "scenario", "rcp", "item")},
+      getSets(data) <- c("region", "year", "scenario", "rcp", "item")
+    },
     Capital = {
       data <- read.csv(file.path(ver, "capitalProjections.csv"))
       data <- as.magpie(data)
       data <- collapseNames(data)
       getItems(data, 3.1) <- sub("gdp_", "", getItems(data, 3.1))
-      getSets(data) <- c("region", "year", "scenario", "variable")},
+      getSets(data) <- c("region", "year", "scenario", "variable")
+    },
     CapitalUnit = {
       mcapitalunitCap <- read.csv(file.path(ver, "capitalUnitCost_cap.csv"))
       mcapitalunitCap$type <- "cap"
@@ -62,18 +64,14 @@ readEDGE <- function(subtype = c("FE_stationary", "FE_buildings", "Capital", "Ca
       data <- rbind(mcapitalunitCap, mcapitalunitInv)
       data <- data[c(setdiff(colnames(data), "value"), "value")]
       data <- as.magpie(data, tidy = TRUE)
-      data <- collapseNames(data)},
+      data <- collapseNames(data)
+    },
     Floorspace = {
       data <- read.csv(file.path(ver, "EDGE_buildings_floorspace.csv"))
       data <- as.magpie(data)
       data <- collapseNames(data)
-      getSets(data) <- c("region", "year", "scenario", "variable")},
-    ES_buildings = {
-      data <- read.csv(file.path(ver, "EDGE_buildings_service.csv"))
-      data <- as.magpie(data)
-      # Only consider trajectories with fixed climate for services
-      data <- mselect(data, rcp = "rcpNoC", collapseNames = TRUE)
-      getSets(data) <- c("region", "year", "scenario", "item")}
+      getSets(data) <- c("region", "year", "scenario", "variable")
+    }
   )
 
   if ("scenario" %in% getSets(data)) {