diff --git a/.buildlibrary b/.buildlibrary
index b751bb4..8ad7c82 100644
--- a/.buildlibrary
+++ b/.buildlibrary
@@ -1,4 +1,4 @@
-ValidationKey: '2742180'
+ValidationKey: '2955900'
 AutocreateReadme: yes
 AcceptedWarnings:
 - 'Warning: package ''.*'' was built under R version'
diff --git a/.github/workflows/check.yaml b/.github/workflows/check.yaml
index 7d564a1..870f216 100644
--- a/.github/workflows/check.yaml
+++ b/.github/workflows/check.yaml
@@ -23,6 +23,7 @@ jobs:
       - uses: r-lib/actions/setup-r-dependencies@v2
         with:
           extra-packages: |
+            gamstransfer=?ignore
             any::lucode2
             any::covr
             any::madrat
@@ -56,6 +57,8 @@ jobs:
 
       - name: Test coverage
         shell: Rscript {0}
-        run: covr::codecov(quiet = FALSE)
+        run: |
+          nonDummyTests <- setdiff(list.files("./tests/testthat/"), c("test-dummy.R", "_snaps"))
+          if(length(nonDummyTests) > 0) covr::codecov(quiet = FALSE)
         env:
           NOT_CRAN: "true"
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 5d2e4ca..2f13466 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -3,7 +3,7 @@
 exclude: '^tests/testthat/_snaps/.*$'
 repos:
 -   repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.4.0
+    rev: v4.5.0
     hooks:
     -   id: check-case-conflict
     -   id: check-json
@@ -15,7 +15,7 @@ repos:
     -   id: mixed-line-ending
 
 -   repo: https://github.com/lorenzwalthert/precommit
-    rev: v0.3.2.9019
+    rev: v0.3.2.9025
     hooks:
     -   id: parsable-R
     -   id: deps-in-desc
diff --git a/CITATION.cff b/CITATION.cff
index 64a9feb..795b8fd 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -2,8 +2,8 @@ cff-version: 1.2.0
 message: If you use this software, please cite it using the metadata from this file.
 type: software
 title: 'mrdrivers: Create GDP and Population Scenarios'
-version: 1.4.0
-date-released: '2023-08-18'
+version: 1.5.0
+date-released: '2023-12-15'
 abstract: Create GDP and population scenarios This package constructs the GDP and
   population scenarios used as drivers in both the REMIND and MAgPIE models.
 authors:
diff --git a/DESCRIPTION b/DESCRIPTION
index ed7f300..4cd8afa 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: mrdrivers
 Type: Package
 Title: Create GDP and Population Scenarios
-Version: 1.4.0
+Version: 1.5.0
 Authors@R: c(person(given = "Johannes", 
                     family = "Koch", 
                     email = "jokoch@pik-potsdam.de", 
@@ -19,7 +19,6 @@ Depends:
     madrat (>= 2.5.1),
     magclass (>= 6.0.3)
 Imports: 
-    bezier,
     countrycode,
     dplyr,
     lifecycle,
@@ -49,6 +48,6 @@ Suggests:
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
-Date: 2023-08-18
+Date: 2023-12-15
 Config/testthat/edition: 3
 VignetteBuilder: knitr
diff --git a/R/calcDriver.R b/R/calcDriver.R
index 400894d..3c03058 100644
--- a/R/calcDriver.R
+++ b/R/calcDriver.R
@@ -121,21 +121,34 @@ calcScenarioConstructor <- function(driver,
   harmonizedData$x <- toolFinishingTouches(x = harmonizedData$x, extension2150 = extension2150, naming = naming)
 
   weight <- NULL
+  description <- harmonizedData$description
   if (popAsWeight) {
-    weight <- calcOutput("Population", scenario = scenario, extension2150 = extension2150, aggregate = FALSE)
+    weight <- calcOutput("Population",
+                         scenario = scenario,
+                         extension2150 = extension2150,
+                         aggregate = FALSE,
+                         supplementary = TRUE)
     # Give weight same names as data, so that aggregate doesn't mess up data dim
-    getNames(weight) <- getNames(harmonizedData$x)
+    getNames(weight$x) <- getNames(harmonizedData$x)
     # Make sure weight and harmonizedData have the same yearly resolution. Sometimes x has more years than weigth,
     # thus the intersect operation. Then if weight has more years than x, only years that exist in x are used.
     # (this applies specifically to the noCovid and ISIMIP scenarios)
-    harmonizedData$x <- harmonizedData$x[, intersect(getYears(harmonizedData$x), getYears(weight)), ]
-    weight <- weight[, getYears(harmonizedData$x), ]
+    harmonizedData$x <- harmonizedData$x[, intersect(getYears(harmonizedData$x), getYears(weight$x)), ]
+    weight$x <- weight$x[, getYears(harmonizedData$x), ]
+    description <- glue("{description} Associated {weight$description}")
+  }
+
+  if (extension2150 == "bezier") {
+    description <- glue("{description} Extended from 2100 to 2150 using bezier curves, resulting in a smooth \\
+                        flattening of the scenario (the slope in 2150 is equal to half of that in 2100).")
+  } else if (extension2150 == "constant") {
+    description <- glue("{description} Extended from 2100 to 2150 using the constant 2100 value.")
   }
 
   list(x = harmonizedData$x,
-       weight = weight,
+       weight = weight$x,
        unit = harmonizedData$unit,
-       description = glue("{driver} {scenario} data: {harmonizedData$description}"))
+       description = glue("{driver} {scenario} scenarios: {description}"))
 }
 
 
diff --git a/R/calcGDP.R b/R/calcGDP.R
index d743e96..a523f54 100644
--- a/R/calcGDP.R
+++ b/R/calcGDP.R
@@ -44,6 +44,7 @@
 #' @examples \dontrun{
 #' # Return default scenarios
 #' calcOutput("GDP")
+#' calcOutput("GDPpc")
 #'
 #' # Return only the SSP2EU GDP scenario
 #' calcOutput("GDP", scenario = "SSP2EU")
@@ -79,7 +80,7 @@ calcGDP <- function(scenario = c("SSPs", "SDPs", "SSP2EU"),
     average2020 <- FALSE
   }
   if (average2020) {
-    # For REMIND, the concensus is to avergae the 2020 value so as to dampen the effect of the COVID shock. (The
+    # For REMIND, the consensus is to average the 2020 value so as to dampen the effect of the COVID shock. (The
     # reasoning being that REMIND uses 5-year time steps, and that the year-in-itself should represent the 2,5 years
     # before and after.)
     xNew2020 <- (gdp$x[, 2018, ] + gdp$x[, 2019, ] + gdp$x[, 2020, ] + gdp$x[, 2021, ] + gdp$x[, 2022, ]) / 5
@@ -89,7 +90,7 @@ calcGDP <- function(scenario = c("SSPs", "SDPs", "SSP2EU"),
     gdp$description <- paste(gdp$description, "|| 2020 value averaged over 2018-2022 time period.")
     # Return only 5 year time steps, since the yearly data around 2020 is not connected to the 2020 value anymore.
     years5ts <- getYears(gdp$x, as.integer = TRUE)[getYears(gdp$x, as.integer = TRUE) %% 5 == 0 &
-                                                   getYears(gdp$x, as.integer = TRUE) != 1960]
+                                                     getYears(gdp$x, as.integer = TRUE) != 1960]
     gdp$x <- gdp$x[, years5ts, ]
     gdp$weight <- gdp$weight[, years5ts, ]
     gdp$description <- paste(gdp$description, "5 year time steps only.")
diff --git a/R/calcGDPFuture.R b/R/calcGDPFuture.R
index 515ff39..b7feab8 100644
--- a/R/calcGDPFuture.R
+++ b/R/calcGDPFuture.R
@@ -26,16 +26,15 @@ calcGDPFuture <- function(GDPFuture = "SSPs-MI", unit = "constant 2005 Int$PPP")
 calcInternalGDPFuture <- function(GDPFuture, unit) { # nolint
   data <- switch(
     GDPFuture,
-    "SSPs"   = calcOutput("InternalGDPFutureSSPs", unit = unit, aggregate = FALSE),
-    "SSP2EU" = calcOutput("InternalGDPFutureSSP2EU", unit = unit, aggregate = FALSE),
-    "SDPs"   = calcOutput("InternalGDPFutureSDPs", unit = unit, aggregate = FALSE),
-    "MI"     = calcOutput("InternalGDPMI", unit = unit, aggregate = FALSE),
+    "SSPs"   = calcOutput("InternalGDPFutureSSPs", unit = unit, aggregate = FALSE, supplementary = TRUE),
+    "SSP2EU" = calcOutput("InternalGDPFutureSSP2EU", unit = unit, aggregate = FALSE, supplementary = TRUE),
+    "SDPs"   = calcOutput("InternalGDPFutureSDPs", unit = unit, aggregate = FALSE, supplementary = TRUE),
+    "MI"     = calcOutput("InternalGDPMI", unit = unit, aggregate = FALSE, supplementary = TRUE),
     stop("Bad input for calcGDPFuture. Invalid 'GDPFuture' argument.")
   )
 
-  data <- toolFinishingTouches(data)
-
-  list(x = data, weight = NULL, unit = glue("mil. {unit}"), description = glue("GDP data from {GDPFuture}"))
+  data$x <- toolFinishingTouches(data$x)
+  data
 }
 
 
@@ -77,7 +76,7 @@ calcInternalGDPFutureSSPs <- function(unit) {
       GDPuc::convertGDP("constant 2005 Int$PPP", unit, replace_NAs = c("linear", "no_conversion"))
 
     y2 <- getYears(data2005PPP)[getYears(data2005PPP, as.integer = TRUE) > c15 &
-                                 getYears(data2005PPP, as.integer = TRUE) < 2100]
+                                  getYears(data2005PPP, as.integer = TRUE) < 2100]
     dataFarFut <- data2005PPP[, y2, ] * NA
 
     # Convert to 2017 Int$PPP using the 2017 value of base 2005 GDP deflator
@@ -87,14 +86,14 @@ calcInternalGDPFutureSSPs <- function(unit) {
 
     data2017PPP <- mbind(dataNearFut, dataFarFut, data2100)
 
-    q <- data2005PPP / data2017PPP
+    ratio <- data2005PPP / data2017PPP
     # For interpolation to work, the last and first values have to be non-NA/non-NaN
-    q[, 2100, ][is.na(q[, 2100, ])] <- 0
+    ratio[, 2100, ][is.na(ratio[, 2100, ])] <- 0
     # The first 2 years of the SSP data set are incomplete. For countries that only lack data in these first 2 years,
     # set NaN to 0.
-    q[, 2000, ][is.nan(q[, 2000, ]) & !is.nan(q[, 2010, ])] <- 0
+    ratio[, 2000, ][is.nan(ratio[, 2000, ]) & !is.nan(ratio[, 2010, ])] <- 0
 
-    q <- as.data.frame(q, rev = 2) %>%
+    ratio <- as.data.frame(ratio, rev = 2) %>%
       dplyr::rename("value" = ".value") %>%
       dplyr::arrange(.data$year) %>%
       dplyr::group_by(.data$iso3c, .data$variable) %>%
@@ -102,7 +101,7 @@ calcInternalGDPFutureSSPs <- function(unit) {
       dplyr::ungroup() %>%
       as.magpie(tidy = TRUE)
 
-    data2017PPP <- data2005PPP / q
+    data2017PPP <- data2005PPP / ratio
     data2017PPP[is.na(data2017PPP)] <- data2005PPP[is.na(data2017PPP)]
     # Above should probably be "<- 0"
     ##################
@@ -111,10 +110,10 @@ calcInternalGDPFutureSSPs <- function(unit) {
 
   # If unit was in $MER
   if (constructUnit != unit) {
-     data <- GDPuc::convertGDP(data, constructUnit, unit, replace_NAs = c("linear", "no_conversion"))
+    data <- GDPuc::convertGDP(data, constructUnit, unit, replace_NAs = c("linear", "no_conversion"))
   }
 
-  list(x = data, weight = NULL, unit = unit, description = "GDP data from SSPs")
+  list(x = data, weight = NULL, unit = glue("mil. {unit}"), description = "SSP projections")
 }
 
 calcInternalGDPFutureSDPs <- function(unit) {
@@ -124,12 +123,12 @@ calcInternalGDPFutureSDPs <- function(unit) {
                      ~ setNames(dataSSP1, gsub("SSP1", .x, getNames(dataSSP1)))) %>%
     mbind()
 
-  list(x = data, weight = NULL, unit = unit, description = "GDP data from SDPs")
+  list(x = data, weight = NULL, unit = glue("mil. {unit}"), description = "SDP projections")
 }
 
 calcInternalGDPFutureSSP2EU <- function(unit) {
   dataSSP2EU <- readSource("ARIADNE", "gdp_corona") %>%
-      GDPuc::convertGDP("constant 2005 Int$PPP", unit, replace_NAs = c("linear", "no_conversion"))
+    GDPuc::convertGDP("constant 2005 Int$PPP", unit, replace_NAs = c("linear", "no_conversion"))
   dataSSP <- calcOutput("InternalGDPFutureSSPs", unit = unit, aggregate = FALSE)
 
   # Get EU-27 countries
@@ -143,11 +142,11 @@ calcInternalGDPFutureSSP2EU <- function(unit) {
   data <- dataSSP[, , "gdp_SSP2"] %>% setNames("gdp_SSP2EU")
   data[euCountries, , ] <- 0
   data[euCountries, cy, ] <- dataSSP2EU[euCountries, cy, ]
-  list(x = data, weight = NULL, unit = unit, description = "GDP data from ARIADNE")
+  list(x = data, weight = NULL, unit = glue("mil. {unit}"), description = "ARIADNE projections")
 }
 
 calcInternalGDPMI <- function(unit) {
   data <- readSource("MissingIslands", "gdp") %>%
     GDPuc::convertGDP("constant 2005 Int$PPP", unit, replace_NAs = c("linear", "no_conversion"))
-  list(x = data, weight = NULL, unit = unit, description = "GDP data from MI")
+  list(x = data, weight = NULL, unit = glue("mil. {unit}"), description = "MI projections")
 }
diff --git a/R/calcGDPHarmonized.R b/R/calcGDPHarmonized.R
index dfb58ea..ab256d5 100644
--- a/R/calcGDPHarmonized.R
+++ b/R/calcGDPHarmonized.R
@@ -46,28 +46,27 @@ toolGDPHarmonizeSSP2EU <- function(past, future, unit) {
                      supplementary = TRUE)
   ssp2Data <- ssps$x[, , "gdp_SSP2"]
 
-  # For SSP2EU: simply glue past with future
-  # Get EUR countries.
+  # For SSP2EU: EU countries are equal to past
   euCountries <- toolGetEUcountries(onlyWithARIADNEgdpData = TRUE)
-  ssp2EUData <- ssp2Data[, getYears(ssp2Data)[getYears(ssp2Data, as.integer = TRUE) <= 2100], ]
-  ssp2EUData[euCountries, , ] <- 0
-  ssp2EUData[euCountries, getYears(past$x), ] <- past$x[euCountries, , ]
+  ssp2EUData <- ssp2Data[euCountries, getYears(ssp2Data)[getYears(ssp2Data, as.integer = TRUE) <= 2100], ]
+  ssp2EUData[, , ] <- 0
+  ssp2EUData[, getYears(past$x), ] <- past$x[euCountries, , ]
 
   # Use GDP growth rates of eurostat for the years 2022, 2023, 2024
   gr <- readSource("EurostatPopGDP", "GDPgr_projections")
-  ssp2EUData[euCountries, 2022, ] <- ssp2EUData[euCountries, 2021, ] * (1 + gr[euCountries, 2022, ] / 100)
-  ssp2EUData[euCountries, 2023, ] <- ssp2EUData[euCountries, 2022, ] * (1 + gr[euCountries, 2023, ] / 100)
-  ssp2EUData[euCountries, 2024, ] <- ssp2EUData[euCountries, 2023, ] * (1 + gr[euCountries, 2024, ] / 100)
+  ssp2EUData[, 2022, ] <- ssp2EUData[, 2021, ] * (1 + gr[euCountries, 2022, ] / 100)
+  ssp2EUData[, 2023, ] <- ssp2EUData[, 2022, ] * (1 + gr[euCountries, 2023, ] / 100)
+  ssp2EUData[, 2024, ] <- ssp2EUData[, 2023, ] * (1 + gr[euCountries, 2024, ] / 100)
 
   # After 2024 use growth rates from future object
   pastYears <- getYears(ssp2EUData)[getYears(ssp2EUData, as.integer = TRUE) <= 2024]
   cy <- union(pastYears, getYears(future$x))
-  ssp2EUData[euCountries, cy, ] <- toolHarmonizePastGrFuture(ssp2EUData[euCountries, pastYears, ],
-                                                             future$x[euCountries, , ])
+  ssp2EUData[, cy, ] <- toolHarmonizePastGrFuture(ssp2EUData[, pastYears, ],
+                                                  future$x[euCountries, , ])
 
   # After 2070, transition to SSP2 values by 2150
   pastYears <- getYears(ssp2EUData)[getYears(ssp2EUData, as.integer = TRUE) <= 2070]
-  combinedSSP2EU <- toolHarmonizePastTransition(ssp2EUData[euCountries, pastYears, ],
+  combinedSSP2EU <- toolHarmonizePastTransition(ssp2EUData[, pastYears, ],
                                                 ssp2Data[euCountries, , ],
                                                 2150)
 
@@ -79,7 +78,9 @@ toolGDPHarmonizeSSP2EU <- function(past, future, unit) {
   combined <- combined[, getYears(combined)[getYears(combined, as.integer = TRUE) <= 2100], ]
 
   list(x = combined,
-      description = glue("Equal to SSP2 in all countries except for EUR countries. For EUR countries glue past data \\
-                          ({past$description}) with future data ({future$description}) and after 2070 converge to \\
-                          2150 SSP2 values."))
+       description = glue("equal to SSP2 in all countries except for EU countries. \\
+                          For EU countries use {past$description} until 2021, \\
+                          growth rates from Eurostat until 2024, \\
+                          growth rates from {future$description} until 2070, \\
+                          and converge to 2150 (bezier-extended) SSP2 values thereafter."))
 }
diff --git a/R/calcGDPPast.R b/R/calcGDPPast.R
index 492acd4..3431d50 100644
--- a/R/calcGDPPast.R
+++ b/R/calcGDPPast.R
@@ -57,7 +57,7 @@ calcInternalGDPPast <- function(GDPPast, unit) { # nolint
 
   data <- toolFinishingTouches(data)
 
-  list(x = data, weight = NULL, unit = glue("mil. {unit}"), description = glue("GDP data from {GDPPast}."))
+  list(x = data, weight = NULL, unit = glue("mil. {unit}"), description = glue("{GDPPast} data"))
 }
 
 
@@ -69,12 +69,12 @@ calcInternalGDPPast <- function(GDPPast, unit) { # nolint
 calcInternalGDPPastWDI <- function(unit) {
   # "NY.GDP.MKTP.PP.KD" = GDP in constant 2017 Int$PPP (as of time of writing this function)
   data <- readSource("WDI", "NY.GDP.MKTP.PP.KD") %>%
-     GDPuc::convertGDP("constant 2017 Int$PPP", unit, replace_NAs = c("linear", "no_conversion"))
+    GDPuc::convertGDP("constant 2017 Int$PPP", unit, replace_NAs = c("linear", "no_conversion"))
 
   data <- fillWithWBFromJames2019(data, unit)
 
   getNames(data) <- glue("gdp in {unit}")
-  list(x = data, weight = NULL, unit = unit, description = "GDP from WDI")
+  list(x = data, weight = NULL, unit = unit, description = "WDI data")
 }
 
 calcInternalGDPPastEurostat <- function(unit) {
@@ -89,7 +89,7 @@ calcInternalGDPPastEurostat <- function(unit) {
   data <- data %>% toolCountryFill(fill = 0) %>% suppressMessages()
 
   getNames(data) <- glue("gdp in {unit}")
-  list(x = data, weight = NULL, unit = unit, description = "GDP from Eurostat")
+  list(x = data, weight = NULL, unit = unit, description = "Eurostat data")
 }
 
 calcInternalGDPPastJames <- function(subtype) {
diff --git a/R/calcGDPpcFuture.R b/R/calcGDPpcFuture.R
index 63e4c70..c19a71a 100644
--- a/R/calcGDPpcFuture.R
+++ b/R/calcGDPpcFuture.R
@@ -28,7 +28,7 @@ calcGDPpcFuture <- function(GDPpcFuture = "SSPsOld-MI", # nolint
 
   data <- data[, getYears(weight), ]
 
-  list(x = data, weight = weight, unit = unit, description = glue("GDPpc data from {GDPpcFuture}"))
+  list(x = data, weight = weight, unit = unit, description = glue("{GDPpcFuture} projections"))
 }
 
 toolGDPpcFutureSSPsOld <- function(unit, mi = FALSE) {
diff --git a/R/calcGDPpcHarmonized.R b/R/calcGDPpcHarmonized.R
index 37ef42b..b45fef7 100644
--- a/R/calcGDPpcHarmonized.R
+++ b/R/calcGDPpcHarmonized.R
@@ -5,7 +5,8 @@
 #' @inherit madrat::calcOutput return
 #' @keywords internal
 calcGDPpcHarmonized <- function(args) {
-  harmonizedData <- switch(args$harmonization,
+  harmonizedData <- switch(
+    args$harmonization,
     "calibSSPs"       = toolGDPpcHarmonizeSSP(args$past, args$future, args$unit, yEnd = 2100),
     "calibSDPs"       = toolGDPpcHarmonizeSDP(args$unit),
     "GDPoverPop"      = toolDivideGDPbyPop(args$scenario, args$unit),
@@ -83,19 +84,22 @@ toolGDPpcHarmonizeSSP <- function(past, future, unit, yEnd, noCovid = FALSE) {
     as.magpie() %>%
     toolCountryFill(fill = 0)
 
+  lastPastYear <- max(getYears(past$x, as.integer = TRUE))
+
   list(x = combinedGDPpc,
-       description = glue("use past data ({past$description}) until {yStart}, short term growth rates from IMF \\
-                           and afterwards transition to future data ({future$description}) by {yEnd}."))
+       description = glue("use {past$description} until {lastPastYear}, \\
+                          short term growth rates from IMF until {yStart}, \\
+                          and transition to {future$description} by {yEnd}."))
 }
 
 toolGDPpcHarmonizeSDP <- function(unit) {
 
   gdppcapSSP1 <- calcOutput("GDPpc",
-                             scenario = "SSPs",
-                             unit = unit,
-                             extension2150 = "none",
-                             average2020 = FALSE,
-                             aggregate = FALSE)[, , "gdppc_SSP1"]
+                            scenario = "SSPs",
+                            unit = unit,
+                            extension2150 = "none",
+                            average2020 = FALSE,
+                            aggregate = FALSE)[, , "gdppc_SSP1"]
 
   # standard SDP inherits SSP1 GDP
   gdppcapSDP <- gdppcapSSP1
@@ -144,7 +148,7 @@ toolGDPpcHarmonizeShortCovid <- function(unit) {
                           average2020 = FALSE,
                           aggregate = FALSE)
 
-  gdppcNoCovid <- calcOutput("GDPpc",
+  gdppcNoCovid <- calcOutput("GDPpc", # nolint
                              scenario = "noCovid",
                              unit = unit,
                              extension2150 = "none",
@@ -170,7 +174,7 @@ toolGDPpcHarmonizeLongCovid <- function(unit) {
                           average2020 = FALSE,
                           aggregate = FALSE)
 
-  gdppcNoCovid <- calcOutput("GDPpc",
+  gdppcNoCovid <- calcOutput("GDPpc", # nolint
                              scenario = "noCovid",
                              unit = unit,
                              extension2150 = "none",
@@ -225,29 +229,29 @@ convergeSpecial <- function(x, yearStart, yearEnd) {
       ),
       # Fast convergence
       value = dplyr::if_else(
-         .data$year > yearStart &
-            ((.data$SSP %in% c("SSP1", "SSP5") & .data$d >= 0) |
+        .data$year > yearStart &
+          ((.data$SSP %in% c("SSP1", "SSP5") & .data$d >= 0) |
              (.data$SSP %in% c("SSP3", "SSP4") & .data$d < 0)),
-         dplyr::if_else(
-             .data$year <= yearEnd,
-             .data$iiasa_gdppc - .data$d * (yearEnd - .data$year) / (yearEnd - yearStart),
-             .data$iiasa_gdppc
-         ),
-         .data$value
+        dplyr::if_else(
+          .data$year <= yearEnd,
+          .data$iiasa_gdppc - .data$d * (yearEnd - .data$year) / (yearEnd - yearStart),
+          .data$iiasa_gdppc
+        ),
+        .data$value
       ),
       # Slow convergence
       value = dplyr::if_else(
         .data$year > yearStart &
-            ((.data$SSP %in% c("SSP3", "SSP4") & .data$d >= 0) |
+          ((.data$SSP %in% c("SSP3", "SSP4") & .data$d >= 0) |
              (.data$SSP %in% c("SSP1", "SSP5") & .data$d < 0)),
         dplyr::if_else(
-            .data$year <= y2,
-            .data$iiasa_gdppc - .data$d,
-            dplyr::if_else(
-              .data$year <= yearEnd,
-              .data$iiasa_gdppc - .data$d * (yearEnd - .data$year) / (yearEnd - y2),
-              .data$iiasa_gdppc
-            )
+          .data$year <= y2,
+          .data$iiasa_gdppc - .data$d,
+          dplyr::if_else(
+            .data$year <= yearEnd,
+            .data$iiasa_gdppc - .data$d * (yearEnd - .data$year) / (yearEnd - y2),
+            .data$iiasa_gdppc
+          )
         ),
         .data$value
       ),
@@ -310,9 +314,9 @@ computeSHAPEgrowth <- function(shapeGDPScenario, gdppcapSSP1, startFromYear) {
         modificationFactor <- logisticTransition(
           gdppcap[, yr, ], l0 = 1.15, l = 1, k = 20, x0 = 15e3, useLog10 = TRUE
         )
-      }
-      # service-driven (SDP_MC): growth rate reduced based on relative distance to technology frontier (given by the US)
-      else if (shapeGDPScenario == "gdppc_SDP_MC") {
+        # service-driven (SDP_MC): growth rate reduced based on relative distance to technology frontier
+        # (given by the US)
+      } else if (shapeGDPScenario == "gdppc_SDP_MC") {
         # define US as technology frontier
         frontier <- gdppcap["USA", yr, ]
         getItems(frontier, 1) <- "GLO"
@@ -321,9 +325,8 @@ computeSHAPEgrowth <- function(shapeGDPScenario, gdppcapSSP1, startFromYear) {
         modificationFactor <- logisticTransition(
           reldiff2frontier[, yr, ], l0 = 1, l = 0.5, k = -30, x0 = 0.2, useLog10 = FALSE
         )
-      }
-      # society-driven (SDP_RC): gradual transition to zero growth for high-income countries
-      else if (shapeGDPScenario == "gdppc_SDP_RC") {
+        # society-driven (SDP_RC): gradual transition to zero growth for high-income countries
+      } else if (shapeGDPScenario == "gdppc_SDP_RC") {
         modificationFactor <- logisticTransition(gdppcap[, yr, ], l0 = 1, l = 0, k = 10, x0 = 30e3, useLog10 = TRUE)
       } else {
         stop("cannot create SHAPE GDP scenarios: unknown scenario")
diff --git a/R/calcGDPpcPast.R b/R/calcGDPpcPast.R
index 00e253f..c4121fc 100644
--- a/R/calcGDPpcPast.R
+++ b/R/calcGDPpcPast.R
@@ -18,7 +18,7 @@ calcGDPpcPast <- function(GDPpcPast = "WDI-MI", unit = "constant 2005 Int$PPP")
                        PopulationPast = GDPpcPast,
                        aggregate = FALSE)
 
-  list(x = data, weight = weight, unit = unit, description = glue("GDPpc data from {GDPpcPast}."))
+  list(x = data, weight = weight, unit = unit, description = glue("{GDPpcPast} data"))
 }
 
 cGDPpcFromGDPAndPop <- function(GDPpcPast, unit) { # nolint
diff --git a/R/calcPopulationFuture.R b/R/calcPopulationFuture.R
index ed138f5..7394879 100644
--- a/R/calcPopulationFuture.R
+++ b/R/calcPopulationFuture.R
@@ -24,18 +24,19 @@ calcPopulationFuture <- function(PopulationFuture = "SSPs-UN_PopDiv-MI") { # nol
 # Internal Function
 ######################################################################################
 calcInternalPopulationFuture <- function(PopulationFuture) { # nolint
-  data <- switch(PopulationFuture,
-                 "SSPs"      = calcOutput("InternalPopulationFutureSSPs", aggregate = FALSE),
-                 "SSP2EU"    = calcOutput("InternalPopulationFutureSSP2EU", aggregate = FALSE),
-                 "SDPs"      = calcOutput("InternalPopulationFutureSDPs", aggregate = FALSE),
-                 "UN_PopDiv" = calcOutput("InternalPopulationFutureUN_PopDiv", aggregate = FALSE),
-                 "MI"        = readSource("MissingIslands", "pop"),
-                 "SSPsOld"   = calcOutput("InternalPopulationFutureSSPsOld", aggregate = FALSE),
-                 stop("Bad input for PopulationFuture. Invalid 'PopulationFuture' argument."))
-
-  data <- toolFinishingTouches(data)
-
-  list(x = data, weight = NULL, unit = "million", description = glue("Population {PopulationFuture} data"))
+  data <- switch(
+    PopulationFuture,
+    "SSPs"      = calcOutput("InternalPopulationFutureSSPs", aggregate = FALSE, supplementary = TRUE),
+    "SSP2EU"    = calcOutput("InternalPopulationFutureSSP2EU", aggregate = FALSE, supplementary = TRUE),
+    "SDPs"      = calcOutput("InternalPopulationFutureSDPs", aggregate = FALSE, supplementary = TRUE),
+    "UN_PopDiv" = calcOutput("InternalPopulationFutureUN_PopDiv", aggregate = FALSE, supplementary = TRUE),
+    "MI"        = calcOutput("InternalPopMI", aggregate = FALSE, supplementary = TRUE),
+    "SSPsOld"   = calcOutput("InternalPopulationFutureSSPsOld", aggregate = FALSE, supplementary = TRUE),
+    stop("Bad input for PopulationFuture. Invalid 'PopulationFuture' argument.")
+  )
+
+  data$x <- toolFinishingTouches(data$x)
+  data
 }
 
 
@@ -46,7 +47,7 @@ calcInternalPopulationFuture <- function(PopulationFuture) { # nolint
 calcInternalPopulationFutureSSPs <- function() { # nolint
   data <- readSource("SSP", "pop2018Update") * 1e-3
   getNames(data) <- paste0("pop_", getNames(data))
-  list(x = data, weight = NULL, unit = "million", description = "Population data from SSP")
+  list(x = data, weight = NULL, unit = "million", description = "SSP projections")
 }
 
 calcInternalPopulationFutureSDPs <- function() { # nolint
@@ -55,7 +56,7 @@ calcInternalPopulationFutureSDPs <- function() { # nolint
   data <- purrr::map(c("SDP", "SDP_EI", "SDP_RC", "SDP_MC"),
                      ~ setNames(data_SSP1, gsub("SSP1", .x, getNames(data_SSP1)))) %>%
     mbind()
-  list(x = data, weight = NULL, unit = "million", description = "Population data from SDP")
+  list(x = data, weight = NULL, unit = "million", description = "SSP1 projections")
 }
 
 calcInternalPopulationFutureSSP2EU <- function() { # nolint
@@ -74,7 +75,10 @@ calcInternalPopulationFutureSSP2EU <- function() { # nolint
     setNames("pop_SSP2EU")
   data[euCountries, , ] <- 0
   data[euCountries, cy, ] <- dataEurostat[euCountries, cy, ]
-  list(x = data, weight = NULL, unit = "million", description = "Population data from SSP2EU")
+  list(x = data,
+       weight = NULL,
+       unit = "million",
+       description = "SSP2 projections for non-EU countries, and EUROSTAT projections for EU countries")
 }
 
 calcInternalPopulationFutureSSPsOld <- function() { # nolint
@@ -85,11 +89,16 @@ calcInternalPopulationFutureSSPsOld <- function() { # nolint
   getNames(data) <- paste0("pop_", gsub("_v[[:alnum:],[:punct:]]*", "", getNames(data)))
   getNames(data) <- sub("SSP4d", "SSP4", getNames(data))
 
-  list(x = data, weight = NULL, unit = "million", description = "Population data from SSPsOld")
+  list(x = data, weight = NULL, unit = "million", description = "old SSP projections")
 }
 
 calcInternalPopulationFutureUN_PopDiv <- function() { # nolint
   data <- readSource("UN_PopDiv", "medium") * 1e-3
   getNames(data) <- "pop_medium_variant"
-  list(x = data, weight = NULL, unit = "million", description = "Population data from UN_PopDiv")
+  list(x = data, weight = NULL, unit = "million", description = "UN_PopDiv projections")
+}
+
+calcInternalPopMI <- function() {
+  data <- readSource("MissingIslands", "pop")
+  list(x = data, weight = NULL, unit = "million", description = "MI projections")
 }
diff --git a/R/calcPopulationHarmonized.R b/R/calcPopulationHarmonized.R
index 664a2e2..59fa692 100644
--- a/R/calcPopulationHarmonized.R
+++ b/R/calcPopulationHarmonized.R
@@ -27,9 +27,12 @@ toolHarmonizeWithPEAPandFuture <- function(past, future) {
 
   # Use PEAP growth rates until last year of IMF WEO data, and future growth rates after that
   x <- past$x %>% toolHarmonizePastGrFuture(shortTerm) %>% toolHarmonizePastGrFuture(future$x)
+
+  lastPastYear <- max(getYears(past$x, as.integer = TRUE))
   list(x = x,
-       description = glue("use past data from '{past$description}', then the growth rates from the Wolrld Bank's \\
-                           PEAP until {lastYearIMF}, and then the growth rates from '{future$description}'."))
+       description = glue("use {past$description} until {lastPastYear}, \\
+                          growth rates from the Wolrld Bank's PEAP until {lastYearIMF}, \\
+                          and growth rates from {future$description} thereafter."))
 }
 
 toolHarmonizeSSP2EU <- function(past, future) {
@@ -42,7 +45,9 @@ toolHarmonizeSSP2EU <- function(past, future) {
   harmonizedData$x[euCountries, , ] <- x[euCountries, getYears(harmonizedData$x), ]
 
   list(x = harmonizedData$x,
-       description = glue("{harmonizedData$description} For European countries, just use future growth rates."))
+       description = glue("equal to SSP2 in all countries except for EU countries. \\
+                          For EU countries use {past$description} until 2021, \\
+                          and growth rates from projections ({future$description}) thereafter."))
 }
 
 toolHarmonizeISIMIP <- function(past, future, yEnd) {
@@ -55,7 +60,6 @@ toolHarmonizeISIMIP <- function(past, future, yEnd) {
   x <- toolHarmonizePastTransition(past$x, future$x, yEnd)
 
   list(x = x,
-       description = glue("use past data from '{past$description}' - should be UN_PopDiv with data, currently, until \\
-                           2020. Add the 2021 projections from UN_PopDiv. Then converge towards \\
-                           '{future$description}' by {yEnd}."))
+       description = glue("use {past$description} until 2020, UN_PopDiv projections for 2021, \\
+                          and converge towards {future$description} by {yEnd}."))
 }
diff --git a/R/calcPopulationPast.R b/R/calcPopulationPast.R
index 4d840ae..b0d7aef 100644
--- a/R/calcPopulationPast.R
+++ b/R/calcPopulationPast.R
@@ -39,7 +39,7 @@ calcInternalPopulationPast <- function(PopulationPast) { # nolint
   getNames(data) <- "population"
   data <- toolFinishingTouches(data)
 
-  list(x = data, weight = NULL, unit = "million", description = glue("Population {PopulationPast} data"))
+  list(x = data, weight = NULL, unit = "million", description = glue("{PopulationPast} data"))
 }
 
 
@@ -84,5 +84,5 @@ calcInternalPopulationPastEurostat <- function() { # nolint
   }
   # Extrapolate in years outside of WDI range
   dataEurostat <- toolInterpolateAndExtrapolate(dataEurostat)
-  list(x = dataEurostat, weight = NULL, unit = "million", description = "Population data from Eurostat")
+  list(x = dataEurostat, weight = NULL, unit = "million", description = "Eurostat data")
 }
diff --git a/R/readIMF.R b/R/readIMF.R
index 1697741..c8df095 100644
--- a/R/readIMF.R
+++ b/R/readIMF.R
@@ -14,7 +14,7 @@
 readIMF <- function(subtype = "current_account", subset = "WEOOct2022all.xls") {
   # Check function input
   if (!subtype %in% c("current_account", "GDPpc")) {
-    stop("Bad input for readiMD. Invalid 'subtype' argument.")
+    stop("Bad input for readIMF. Invalid 'subtype' argument.")
   }
 
   # Define source file
@@ -34,9 +34,9 @@ readIMF <- function(subtype = "current_account", subset = "WEOOct2022all.xls") {
     dplyr::mutate(value = gsub(",", "", .data$value),
                   dplyr::across(.cols = c("year", "value"),
                                 ~ suppressWarnings(as.double(.x))),
-                                # The warnings that are being suppressed above, come from
-                                # character strings that can't be converted to numeric, and
-                                # are thus returned as NA.
+                  # The warnings that are being suppressed above, come from
+                  # character strings that can't be converted to numeric, and
+                  # are thus returned as NA.
                   value = tidyr::replace_na(.data$value, 0)) %>%
     tidyr::pivot_wider(names_from = "Subject Descriptor")
 
diff --git a/R/readUN_PopDiv.R b/R/readUN_PopDiv.R
index 54f31b2..8e5842f 100644
--- a/R/readUN_PopDiv.R
+++ b/R/readUN_PopDiv.R
@@ -13,7 +13,7 @@ readUN_PopDiv <- function(subtype = "estimates") { # nolint
   }
   file <- "WPP2022_POP_F01_1_POPULATION_SINGLE_AGE_BOTH_SEXES.xlsx"
   sheet <- if (subtype == "estimates") "Estimates" else "Medium variant"
-  readxl::read_xlsx(file, sheet = sheet, skip = 16, col_types = "text") %>%
+  readxl::read_xlsx(file, sheet = sheet, skip = 16, col_types = "text", progress = FALSE) %>%
     dplyr::select("ISO3 Alpha-code", "year" = "Year", dplyr::matches("^[0-9]*$")) %>%
     dplyr::filter(!is.na(.data$`ISO3 Alpha-code`)) %>%
     tidyr::pivot_longer(cols = dplyr::matches("^[0-9]*$"),
@@ -29,14 +29,14 @@ readUN_PopDiv <- function(subtype = "estimates") { # nolint
 #' @rdname readUN_PopDiv
 #' @order 3
 #' @param x MAgPIE object returned from readUN_PopDiv
-convertUN_PopDiv <- function(x) {
+convertUN_PopDiv <- function(x) { # nolint
   toolGeneralConvert(x, no_remove_warning = "XKX")
 }
 
 #' @rdname readUN_PopDiv
 #' @order 1
 downloadUN_PopDiv <- function() { # nolint
-  url <- "https://population.un.org/wpp/Download/Files/1_Indicators%20(Standard)/EXCEL_FILES/2_Population/WPP2022_POP_F01_1_POPULATION_SINGLE_AGE_BOTH_SEXES.xlsx"
+  url <- "https://population.un.org/wpp/Download/Files/1_Indicators%20(Standard)/EXCEL_FILES/2_Population/WPP2022_POP_F01_1_POPULATION_SINGLE_AGE_BOTH_SEXES.xlsx" # nolint
   utils::download.file(url, basename(url), quiet = TRUE)
 
   # Compose meta data
diff --git a/R/toolExtend2150.R b/R/toolExtend2150.R
index 26ab569..41478b2 100644
--- a/R/toolExtend2150.R
+++ b/R/toolExtend2150.R
@@ -1,8 +1,8 @@
 # Extend until 2150 in 5 year time steps. Either with bezierExtension or constant.
-# The bezier extension is only possible if there is data until 2100, and only concerns
-# the years between 2100 and 2150.
 toolExtend2150 <- function(data, extension2150) {
   if (extension2150 != "none") {
+    # The bezier extension is only possible if there is data until 2100, and only affects years between 2100 and 2150.
+    # It extends the time series in such a way as for the slope in 2105 to be half of that in 2100.
     if (extension2150 == "bezier" && "y2100" %in% getYears(data)) {
       data <- bezierExtension(data, seq(2105, 2150, 5))
     } else {
@@ -18,63 +18,54 @@ toolExtend2150 <- function(data, extension2150) {
   data
 }
 
+
 bezierExtension <- function(data, timeExtend) {
+  # Define bezier coordinates
+  bc <- new.magpie(getItems(data, 1), c(2100, 2110, 2140, 2150), getNames(data), fill = 0)
 
-  extension <- new.magpie(getItems(data, 1), timeExtend, getNames(data), fill = 0)
+  slopeStart <- (data[, 2100, ] - data[, 2095, ]) / 5
+  slopeEnd <- slopeStart / 2
 
-  slope <- (data[, 2100, ] - data[, 2095, ]) / 5
+  bc[, 2100, ] <- data[, 2100, ]
+  bc[, 2110, ] <- data[, 2100, ] + slopeStart * 10
+  bc[, 2140, ] <- data[, 2100, ] + slopeEnd * 40
+  bc[, 2150, ] <- data[, 2100, ] + slopeEnd * 50
 
   nr <- nregions(data)
+  nd <- ndata(data)
 
-  for (scen in getNames(data)) {
-    yStart <- data[, 2100, scen] %>% as.matrix()
-    xStart <- matrix(2100, nrow = 249, ncol = 1)
-    slopeStart <- slope[, , scen] %>% as.matrix()
-
-    slopeEnd <- slopeStart / 2
-    yEnd <- yStart + slopeEnd * 50
-    xEnd <- matrix(2150, nrow = nr, ncol = 1)
-
-    x1 <- matrix(2110, nrow = nr, ncol = 1)
-    y1 <- yStart + slopeStart * 10
-    x2 <- matrix(2140, nrow = nr, ncol = 1)
-    y2 <- yEnd - slopeEnd * 10
-
-    for (i in 1:nr) {
-      # If Bezier extension would lead to negative GDP, keep last value constant instead
-      if (yStart[i] == 0 || yEnd[i] < 0) {
-        extension[i, , scen] <- data[i, 2100, scen]
-        next
-      }
+  # If Bezier extension would lead to negative GDP, set bezier coordinates equal to start point
+  # (comes down to a constant extension instead)
+  for (i in 1:nr) for (j in 1:nd) if (bc[i, 2100, j] == 0 || bc[i, 2150, j] < 0) bc[i, , j] <- bc[i, 2100, j]
 
-      p <- matrix(c(xStart[i], yStart[i], x1[i], y1[i], x2[i], y2[i], xEnd[i], yEnd[i]),
-                  nrow = 4,
-                  ncol = 2,
-                  byrow = TRUE)
+  x <- rep(c(2100, 2110, 2140, 2150), nr * nd)
+  y <- purrr::reduce(purrr::map(1:nd, ~purrr::reduce(purrr::map(1:nr, function(y) bc[y, , .x]), c)), c)
+  z <- purrr::reduce(purrr::map(1:(nr * nd), ~rep(.x, 4)), c)
 
-      # Get Bezier curve (Use max_dist method because it's super fast)
-      # Keep this code for later: bp <- bezier::bezier(t = seq(0, 1, length = 10), p = p)
-      pob <- bezier::pointsOnBezier(p = p,
-                                    method = "max_dist",
-                                    max.dist = max(abs(yEnd[i] - yStart[i]), 50) / 100,
-                                    print.progress = FALSE) %>%
-        suppressWarnings()
-      # The above Warnings which are suppressed are
-      # : "essentially perfect fit: summary may be unreliable"
+  # grid::bezierGrob returns the point in a weird graphical unit, and "only" returns 48 points, but is super fast.
+  bezierPoints <- grid::bezierGrob(x, y, id = z) %>% grid::bezierPoints()
+  cfx <- x[1] / as.numeric(bezierPoints[[1]]$x[[1]])
+  cfy <- y[1] / as.numeric(bezierPoints[[1]]$y[[1]])
+  id <- paste(purrr::reduce(purrr::map(getNames(data), ~rep(.x, nr)), c),
+              rep(getItems(data, 1), nd),
+              sep = "-")
 
-      # Get y coordinates of points with x coordinates = timeExtend
-      myBezierOutflow <- pob$points %>%
-        tibble::as_tibble(.name_repair = ~ paste0("V", seq_along(.x))) %>%
-        # Complicatd / elegant use of function factories to get closest points to timeExtend coordinates
-        dplyr::mutate(dplyr::across("V1", purrr::map(timeExtend, ~ function(y) {
-          abs(y - .x)
-        }))) %>%
-        dplyr::filter(dplyr::if_any(tidyselect::contains("_"), ~ .x == min(.x))) %>%
-        dplyr::pull(.data$V2)
-
-      extension[i, , scen] <- myBezierOutflow
-    }
-  }
+  extension <- purrr::map2(bezierPoints, id,
+                           ~.x %>%
+                             tibble::as_tibble() %>%
+                             dplyr::mutate(x = as.numeric(x) * cfx,
+                                           y = as.numeric(y) * cfy,
+                                           id = .y)) %>%
+    purrr::list_rbind() %>%
+    # Complicated / elegant use of function factories to get closest points to timeExtend coordinates
+    # First create columns with distance to timeExtend points
+    dplyr::mutate(dplyr::across("x", purrr::map(timeExtend, ~ function(y) abs(y - .x)))) %>%
+    # Then keep only the rows with the min values of the newly created columns
+    dplyr::filter(dplyr::if_any(tidyselect::contains("_"), ~ .x == min(.x)), .by = "id") %>%
+    dplyr::mutate(year = round(x)) %>%
+    tidyr::separate_wider_delim("id", names = c("data", "iso3c"), delim = "-") %>%
+    dplyr::select("iso3c", "year", "data", "y") %>%
+    as.magpie()
 
   mbind(data, extension)
 }
diff --git a/R/toolGetEUcountries.R b/R/toolGetEUcountries.R
index d156671..27afc84 100644
--- a/R/toolGetEUcountries.R
+++ b/R/toolGetEUcountries.R
@@ -1,4 +1,4 @@
-# These countries are the ones (essentiall the EU-27) that receive special treatment
+# These countries are the ones (essential the EU-27) that receive special treatment
 # in the SSP2EU scenario.
 toolGetEUcountries <- function(onlyWithARIADNEgdpData = FALSE) {
   x <- toolGetMapping("regionmappingH12.csv", type = "regional", where = "mappingfolder") %>%
@@ -7,7 +7,7 @@ toolGetEUcountries <- function(onlyWithARIADNEgdpData = FALSE) {
     dplyr::pull(.data$CountryCode)
 
   if (onlyWithARIADNEgdpData) {
-     x <- x[x %in% where(readSource("ARIADNE", "gdp_corona") != 0)$true$regions]
+    x <- x[x %in% where(readSource("ARIADNE", "gdp_corona") != 0)$true$regions]
   }
   x
 }
diff --git a/R/toolReduce.R b/R/toolReduce.R
index f9d9314..0beb8e7 100644
--- a/R/toolReduce.R
+++ b/R/toolReduce.R
@@ -2,15 +2,24 @@ toolReduce <- function(x, mbindOrFillWith = "mbind") {
   # Combine list elements using mbind and glue.
   if (mbindOrFillWith == "mbind") {
     x %>%
-       purrr::reduce(~ list(x = mbind(.x$x, .y$x),
-                            weight = mbind(.x$weight, .y$weight),
-                            unit = glue("{.x$unit} || {.y$unit}"),
-                            description = glue("{.x$description} || {.y$description}")))
+      purrr::reduce(~ list(x = mbind(.x$x, .y$x),
+                           weight = mbind(.x$weight, .y$weight),
+                           unit = glue("{.x$unit} || {.y$unit}"),
+                           description = glue("{.x$description} || {.y$description}")))
   } else if (mbindOrFillWith == "fillWith") {
-    x %>%
-       purrr::reduce(~ list(x = toolFillWith(.x$x, .y$x),
-                            weight = toolFillWith(.x$weight, .y$weight),
-                            unit = glue("{.x$unit}"),
-                            description = glue("{.x$description} completed with {.y$description}")))
+
+    if (length(x) > 1) {
+      sep <- if (length(x) == 2) " (completed with" else c(" (completed with", c(rep(",", length(x) - 3), " and"))
+      closer <- c(rep("", length(sep) - 1), ")")
+      helper <- purrr::map2(sep, closer, c)
+    } else {
+      helper <- NULL
+    }
+
+    purrr::reduce2(x, helper,
+                   ~ list(x = toolFillWith(.x$x, .y$x),
+                          weight = toolFillWith(.x$weight, .y$weight),
+                          unit = glue("{.x$unit}"),
+                          description = glue("{.x$description}{..3[1]} {.y$description}{..3[2]}")))
   }
 }
diff --git a/README.md b/README.md
index 017eeaf..4831933 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # Create GDP and Population Scenarios
 
-R package **mrdrivers**, version **1.4.0**
+R package **mrdrivers**, version **1.5.0**
 
 [![CRAN status](https://www.r-pkg.org/badges/version/mrdrivers)](https://cran.r-project.org/package=mrdrivers)  [![R build status](https://pik-piam.github.io/mrdrivers/workflows/check/badge.svg)](https://pik-piam.github.io/mrdrivers/actions) [![codecov](https://codecov.io/gh/mrdrivers/branch/master/graph/badge.svg)](https://app.codecov.io/gh/mrdrivers) [![r-universe](https://pik-piam.r-universe.dev/badges/mrdrivers)](https://pik-piam.r-universe.dev/builds)
 
@@ -103,7 +103,7 @@ In case of questions / problems please contact Johannes Koch <jokoch@pik-potsdam
 
 To cite package **mrdrivers** in publications use:
 
-Koch J, Soergel B, Leip D, Benke F, Dietrich J (2023). _mrdrivers: Create GDP and Population Scenarios_. R package version 1.4.0, <https://pik-piam.github.io/mrdrivershttps://github.com/pik-piam/mrdrivers>.
+Koch J, Soergel B, Leip D, Benke F, Dietrich J (2023). _mrdrivers: Create GDP and Population Scenarios_. R package version 1.5.0, <https://pik-piam.github.io/mrdrivershttps://github.com/pik-piam/mrdrivers>.
 
 A BibTeX entry for LaTeX users is
 
@@ -112,7 +112,7 @@ A BibTeX entry for LaTeX users is
   title = {mrdrivers: Create GDP and Population Scenarios},
   author = {Johannes Koch and Bjoern Soergel and Deborra Leip and Falk Benke and Jan Philipp Dietrich},
   year = {2023},
-  note = {R package version 1.4.0},
+  note = {R package version 1.5.0},
   url = {https://pik-piam.github.io/mrdrivers},
   url = {https://github.com/pik-piam/mrdrivers},
 }
diff --git a/man/calcGDP.Rd b/man/calcGDP.Rd
index 7b20f42..ed1186f 100644
--- a/man/calcGDP.Rd
+++ b/man/calcGDP.Rd
@@ -63,6 +63,7 @@ See the vignette: \code{vignette("scenarios")} for scenario options, definitions
 \dontrun{
 # Return default scenarios
 calcOutput("GDP")
+calcOutput("GDPpc")
 
 # Return only the SSP2EU GDP scenario
 calcOutput("GDP", scenario = "SSP2EU")