Proper handling of CO2

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: growR
 Type: Package
-Version: 1.3.0.9006
+Version: 1.3.0.9007
 Date: 2024-05-23
 Authors@R: person(
     given = "Kevin",

NEWS.md

@@ -8,6 +8,7 @@
 * New parameter for irrigation.
 * `growR_run_loop` allows to suppress returning the simulation objects. This 
   is useful to prevent overkill memory usage.
+* CO2 affects plant transpiration, according to `fCO2_transpiration_mod`.
 
 ## Changed
 
@@ -17,6 +18,10 @@
 * Performance improvements to ModvegeSite$un() in the order of 25%.
 * Small performance improvement of `WeatherData` initialization through 
   vectorization of snow model.
+* Set `CO2_growth_factor` from 0.5 to the much more realistic value of 0.1.
+* Atmospheric CO2 for years beyond 2020 is now estimated according to common 
+  representative concentration pathways (RCPs). See, e.g. `WeatherData$set_RCP()`
+  or `WeatherData$aCO2()`.
 
 ## Fixed
 

R/modvegesite.R

@@ -249,6 +249,9 @@ ModvegeSite = R6Class(
         self$year = year
         # Calculate everything that can be done before the loop
         P = self$parameters
+        private$CO2_growth_mod = fCO2_growth_mod(weather$aCO2,
+                                                 P$CO2_growth_factor)
+        private$CO2_transpiration_mod = fCO2_transpiration_mod(weather$aCO2)
         private$initialize_state_variables()
         private$calculate_temperature_sum()
         private$calculate_PETeff()
@@ -550,6 +553,8 @@ ModvegeSite = R6Class(
     minBMGV = NULL,
     minBMGR = NULL,
     WRmin = NULL,
+    CO2_growth_mod = NULL,
+    CO2_transpiration_mod = NULL,
 
     #-Private-methods-----------------------------------------------------------
 
@@ -736,7 +741,7 @@ ModvegeSite = R6Class(
       # Then impose water stress limitation, fW. In the case of soil
       # evaporation use an fW appropriate for high values of PET (PETmx).
       PETeff = private$PETeff[j]
-      PTr = PETeff * (1. - exp(-0.6 * LAI.ET))
+      PTr = PETeff * (1. - exp(-0.6 * LAI.ET)) * private$CO2_transpiration_mod
 
       ATr = PTr * fW(self$WRp / P$WHC, PETeff)
       PEv = PETeff -  PTr
@@ -762,8 +767,7 @@ ModvegeSite = R6Class(
         self$GRO[j] = 0
       } else {
         self$PGRO[j] = W$PAR[j] * P$RUEmax * 
-          (1. - exp(-0.6 * self$LAIGV[j])) * 10. * 
-          fCO2_growth_mod(W$aCO2, P$CO2_growth_factor)
+          (1. - exp(-0.6 * self$LAIGV[j])) * 10. * private$CO2_growth_mod
         self$GRO[j]  = P$NI * self$PGRO[j] * self$ENV[j] * 
           SEA(self$ST[j], P$minSEA, P$maxSEA, P$ST1, P$ST2)
       }

R/parameters.R

@@ -40,7 +40,7 @@ parameter_defaults = list(
   KlDR = 0.0005,
   OMDDV = 0.45,
   OMDDR = 0.4,
-  CO2_growth_factor = 0.5,
+  CO2_growth_factor = 0.1,
   crop_coefficient = 1.15,
   stubble_height = 0.02,
   SGS_method = "MTD",

R/weather.R

@@ -103,10 +103,11 @@ WeatherData = R6Class(
       self$years = years
       # Load weather data
       if (file.exists(weather_file)) {
-        logger(sprintf("Loading weather data from %s.", weather_file), 
+        logger(sprintf("[WeatherData]Loading weather data from %s.", 
+                       weather_file), 
                level = DEBUG)
       } else {
-        stop(sprintf("Weather file `%s` not found.", weather_file))
+        stop(sprintf("[WeatherData]Weather file `%s` not found.", weather_file))
       }
       weather = read.table(weather_file, header = TRUE)
       weather = self$ensure_file_integrity(weather)
@@ -255,7 +256,7 @@ WeatherData = R6Class(
       }
 
       W = list()
-      W[["aCO2"]]  = atmospheric_CO2(year)
+      W[["aCO2"]]  = self$aCO2(year, self$RCP)
       W[["year"]]  = self$year_vec[iW]
       W[["DOY"]]   = self$DOY_vec[iW]
       W[["Ta"]]    = self$Ta_vec[iW]
@@ -270,6 +271,72 @@ WeatherData = R6Class(
 
       self[["W"]] = W
       return(W)
+    },
+
+    #' @description Return the value of `private$RCP`. This is used for the 
+    #' calculation of atmospheric CO2.
+    #'
+    #' @return String representing the currently selected representative 
+    #'   concentration pathway (RCP).
+    #'
+    get_RCP = function() {
+      return(private$RCP)
+    },
+
+    #' @description Choose the RCP to be used for the calculation of 
+    #' atmospheric CO2 concentration.
+    #'
+    #' @param RCP float or string. One of {2.6, 4.5, 8.5}.
+    #'
+    set_RCP = function(RCP) {
+      RCP = as.character(RCP)
+      allowed_RCPs = private$cCO2_table$RCP
+      if (!RCP %in% allowed_RCPs) {
+        message = "[WeatherData]RCP `%s` not recognized. Allowed values:"
+        message = paste(message, paste(allowed_RCPs, collapse = ", "))
+        warning(sprintf(message, RCP))
+      } else {
+        logger(sprintf("[WeatherData]Setting RCP to `%s`.", RCP), level = TRACE)
+        private$RCP = RCP
+      }
+    },
+
+    #' @description Return average atmospheric CO2 concentration for given year.
+    #'
+    #' @param year int The year to consider. For years between 1980 and 2020, 
+    #'   the empirical NOAA fit is taken (see [atmospheric_CO2]). From 2020 
+    #'   to 2100 a rough estimation is made based on the specified *RCP*.
+    #' @param RCP string Name of the RCP in the format "X.Y". Has to be valid 
+    #'   for [WeatherData]`$set_RCP()`.
+    #' @return float The atmospheric CO2 concentration in ppm.
+    #'
+    #' @seealso [atmospheric_CO2]
+    aCO2 = function(year, RCP = NULL) {
+      if (year < 2020) {
+        # Use the empirical function from NOAA.
+        return(atmospheric_CO2(year))
+      } else {
+        # Otherwise, use simple approximations, depending on RCP.
+        # Just interpolate linearly between 2020, 2050 and 2100.
+        if (is.null(RCP)) {
+          RCP = private$RCP
+        }
+        row = private$cCO2_table[private$cCO2_table$RCP == RCP, ]
+        c2020 = atmospheric_CO2(2020)
+        c2050 = row[["c2050"]]
+        c2100 = row[["c2100"]]
+        if (year < 2050) {
+          return((c2050 - c2020)/(2050 - 2020) * (year - 2020) + c2020)
+        } else {
+          return((c2100 - c2050)/(2100 - 2050) * (year - 2050) + c2050)
+        }
+      }
     }
+  ),
+  private = list(
+    RCP = "2.6",
+    cCO2_table = data.frame(RCP = c("2.6", "4.5", "8.5"),
+                            c2050 = c(450, 500, 550),
+                            c2100 = c(420, 600, 900))
   )
 )