Add volcanic scaling factor

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: hector
 Title: Run the Hector Simple Climate Model
-Version: 0.3.0
+Version: 0.4.0
 Authors@R: person("Robert", "Link", email = "robert.link@pnnl.gov", role = c("aut", "cre"))
 Description: Provides R bindings for the Hector Simple Climate Model.
 Depends: R (>= 3.3)

NAMESPACE

@@ -129,11 +129,13 @@ export(TID)
 export(TT)
 export(TU)
 export(TWI)
+export(VOLCANIC_SCALE)
 export(VOLCANIC_SO2)
 export(Y2000_SO2)
 export(enddate)
 export(fetchvars)
 export(getdate)
+export(getname)
 export(getunits)
 export(isactive)
 export(newcore)

R/RcppExports.R

@@ -736,12 +736,18 @@ ECS <- function() {
     .Call('_hector_ECS', PACKAGE = 'hector')
 }
 
-#' @describeIn parameters Aerosol scaling factor (\code{"(unitless)"})
+#' @describeIn parameters Aerosol forcing scaling factor (\code{"(unitless)"})
 #' @export
 AERO_SCALE <- function() {
     .Call('_hector_AERO_SCALE', PACKAGE = 'hector')
 }
 
+#' @describeIn parameters Volcanic forcing scaling factor (\code{"(unitless)"})
+#' @export
+VOLCANIC_SCALE <- function() {
+    .Call('_hector_VOLCANIC_SCALE', PACKAGE = 'hector')
+}
+
 #' @describeIn temperature Global mean temperature
 #' @export
 GLOBAL_TEMP <- function() {

R/hector.R

@@ -133,6 +133,13 @@ enddate <- function(core)
     core$enddate
 }
 
+#' \strong{getname}: Report the name for a Hector instance
+#' @rdname hectorutil
+#' @export
+getname <- function(core)
+{
+    core$name
+}
 
 
 #### Methods for hcore objects

changelog.txt

@@ -1,3 +1,17 @@
+Hector 2.2.0
+============
+* Add a new parameter: VOLCANIC_SCALE.  This parameter adjusts the
+  strength of the response to volcanic forcing.  (PR #291)
+* Add getname function to return the name of a Hector core.
+
+
+Hector 2.1.0
+============
+* Add reset() function to the API.  Calling this function resets a
+  hector core to an earlier time, rerunning the spinup if
+  appropriate.  (PR #243)
+
+
 Hector 2.0.0
 ============
 * Incorporated 1-D diffusive ocean heat model as new temperature component (DOECLIM) (PR #206)

data-raw/units-data.R

@@ -1,14 +1,14 @@
 library(hector)
 
-hvars <- c(ECS(), PREINDUSTRIAL_CO2(), Q10_RH(), BETA(), AERO_SCALE(), DIFFUSIVITY(),
+hvars <- c(ECS(), PREINDUSTRIAL_CO2(), Q10_RH(), BETA(), AERO_SCALE(), VOLCANIC_SCALE(), DIFFUSIVITY(),
            FFI_EMISSIONS(), LUC_EMISSIONS(),
            EMISSIONS_BC(), EMISSIONS_N2O(), EMISSIONS_NOX(), EMISSIONS_CO(),
            EMISSIONS_NMVOC(), EMISSIONS_OC(),
            EMISSIONS_CH4(), PREINDUSTRIAL_CH4(), NATURAL_CH4(), LIFETIME_SOIL(), LIFETIME_STRAT(),
            EMISSIONS_SO2(), VOLCANIC_SO2())
 
 
-hunits <- c('degC', 'ppmv CO2', '(unitless)', '(unitless)', '(unitless)', 'cm2/s',
+hunits <- c('degC', 'ppmv CO2', '(unitless)', '(unitless)', '(unitless)', '(unitless)', 'cm2/s',
             'Pg C/yr', 'Pg C/yr',
             'Tg', 'Tg N', 'Tg N', 'Tg CO',
             'Tg NMVOC', 'Tg',

inst/include/component_data.hpp

@@ -252,6 +252,7 @@
 // temperature component
 #define D_DIFFUSIVITY           "diff"
 #define D_AERO_SCALE            "alpha"
+#define D_VOLCANIC_SCALE        "volscl"
 #define D_FLUX_MIXED            "flux_mixed"
 #define D_FLUX_INTERIOR         "flux_interior"
 #define D_HEAT_FLUX             "heatflux"

inst/include/h_util.hpp

@@ -26,7 +26,7 @@
  * \brief The model version number to be included in logs and outputs.
  * \note This must be updated manually when the model version changes.
  */
-#define MODEL_VERSION "2.1.0"
+#define MODEL_VERSION "2.2.0"
 
 #define OUTPUT_DIRECTORY "output/"
 

inst/include/temperature_component.hpp

@@ -134,6 +134,7 @@ class TemperatureComponent : public IModelComponent {
     unitval S;             //!< climate sensitivity for 2xCO2, deg C
     unitval diff;          //!< ocean heat diffusivity, cm2/s
     unitval alpha;	       //!< aerosol forcing factor, unitless
+    unitval volscl;        //!< volcanic forcing scaling factor, unitless
 
     // Model outputs
     unitval tgav;          //!< global temperature delta, deg C

inst/input/hector-gcam.ini

@@ -139,6 +139,7 @@ enabled=1
 S=3.0 				; equilibrium climate sensitivity for 2xCO2, degC
 diff=2.3			; ocean heat diffusivity, cm2/s
 alpha=1.0     ; scaling factor for aerosol forcing
+volscl=1.0    ; scaling factor for volcanic forcing
 ; Optional global temperature constraint
 ; If supplied, the model will use these data, ignoring what it calculates
 ; tgav_constrain=csv:constraints/tgav_historical.csv

inst/input/hector_rcp26.ini

@@ -158,6 +158,7 @@ enabled=1
 S=3.0 				; equilibrium climate sensitivity for 2xCO2, degC
 diff=2.3			; ocean heat diffusivity, cm2/s
 alpha=1.0     ; scaling factor for aerosol forcing
+volscl=1.0    ; scaling factor for volcanic forcing
 ; Optional global temperature constraint
 ; If supplied, the model will use these data, ignoring what it calculates
 ; tgav_constrain=csv:constraints/tgav_historical.csv

inst/input/hector_rcp26_constrained.ini

@@ -159,6 +159,7 @@ enabled=1
 S=3.0 				; equilibrium climate sensitivity for 2xCO2, degC
 diff=2.3			; ocean heat diffusivity, cm2/s
 alpha=1.0     ; scaling factor for aerosol forcing
+volscl=1.0    ; scaling factor for volcanic forcing
 ; Optional global temperature constraint
 ; If supplied, the model will use these data, ignoring what it calculates
 ; tgav_constrain=csv:constraints/tgav_historical.csv

inst/input/hector_rcp26_histconstrain.ini

@@ -160,6 +160,7 @@ enabled=1
 S=3.0 				; equilibrium climate sensitivity for 2xCO2, degC
 diff=2.3			; ocean heat diffusivity, cm2/s
 alpha=1.0     ; scaling factor for aerosol forcing
+volscl=1.0    ; scaling factor for volcanic forcing
 ; Optional global temperature constraint
 ; If supplied, the model will use these data, ignoring what it calculates
 ; tgav_constrain=csv:constraints/tgav_historical.csv

inst/input/hector_rcp45.ini

@@ -158,6 +158,7 @@ enabled=1
 S=3.0 				; equilibrium climate sensitivity for 2xCO2, degC
 diff=2.3			; ocean heat diffusivity, cm2/s
 alpha=1.0     ; scaling factor for aerosol forcing
+volscl=1.0    ; scaling factor for volcanic forcing
 ; Optional global temperature constraint
 ; If supplied, the model will use these data, ignoring what it calculates
 ; tgav_constrain=csv:constraints/tgav_historical.csv

inst/input/hector_rcp45_constrained.ini

@@ -159,6 +159,7 @@ enabled=1
 S=3.0 				; equilibrium climate sensitivity for 2xCO2, degC
 diff=2.3			; ocean heat diffusivity, cm2/s
 alpha=1.0     ; scaling factor for aerosol forcing
+volscl=1.0    ; scaling factor for volcanic forcing
 ; Optional global temperature constraint
 ; If supplied, the model will use these data, ignoring what it calculates
 ; tgav_constrain=csv:constraints/tgav_historical.csv

inst/input/hector_rcp60.ini

@@ -158,6 +158,7 @@ enabled=1
 S=3.0 				; equilibrium climate sensitivity for 2xCO2, degC
 diff=2.3			; ocean heat diffusivity, cm2/s
 alpha=1.0     ; scaling factor for aerosol forcing
+volscl=1.0    ; scaling factor for volcanic forcing
 ; Optional global temperature constraint
 ; If supplied, the model will use these data, ignoring what it calculates
 ; tgav_constrain=csv:constraints/tgav_historical.csv

inst/input/hector_rcp60_constrained.ini

@@ -159,6 +159,7 @@ enabled=1
 S=3.0 				; equilibrium climate sensitivity for 2xCO2, degC
 diff=2.3			; ocean heat diffusivity, cm2/s
 alpha=1.0     ; scaling factor for aerosol forcing
+volscl=1.0    ; scaling factor for volcanic forcing
 ; Optional global temperature constraint
 ; If supplied, the model will use these data, ignoring what it calculates
 ; tgav_constrain=csv:constraints/tgav_historical.csv

inst/input/hector_rcp85.ini

@@ -158,6 +158,7 @@ enabled=1
 S=3.0 				; equilibrium climate sensitivity for 2xCO2, degC
 diff=2.3			; ocean heat diffusivity, cm2/s
 alpha=1.0     ; scaling factor for aerosol forcing
+volscl=1.0    ; scaling factor for volcanic forcing
 ; Optional global temperature constraint
 ; If supplied, the model will use these data, ignoring what it calculates
 ; tgav_constrain=csv:constraints/tgav_historical.csv

inst/input/hector_rcp85_constrained.ini

@@ -159,6 +159,7 @@ enabled=1
 S=3.0 				; equilibrium climate sensitivity for 2xCO2, degC
 diff=2.3			; ocean heat diffusivity, cm2/s
 alpha=1.0     ; scaling factor for aerosol forcing
+volscl=1.0    ; scaling factor for volcanic forcing
 ; Optional global temperature constraint
 ; If supplied, the model will use these data, ignoring what it calculates
 ; tgav_constrain=csv:constraints/tgav_historical.csv

src/RcppExports.cpp

@@ -1235,6 +1235,16 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// VOLCANIC_SCALE
+String VOLCANIC_SCALE();
+RcppExport SEXP _hector_VOLCANIC_SCALE() {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    rcpp_result_gen = Rcpp::wrap(VOLCANIC_SCALE());
+    return rcpp_result_gen;
+END_RCPP
+}
 // GLOBAL_TEMP
 String GLOBAL_TEMP();
 RcppExport SEXP _hector_GLOBAL_TEMP() {
@@ -1507,6 +1517,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_hector_VOLCANIC_SO2", (DL_FUNC) &_hector_VOLCANIC_SO2, 0},
     {"_hector_ECS", (DL_FUNC) &_hector_ECS, 0},
     {"_hector_AERO_SCALE", (DL_FUNC) &_hector_AERO_SCALE, 0},
+    {"_hector_VOLCANIC_SCALE", (DL_FUNC) &_hector_VOLCANIC_SCALE, 0},
     {"_hector_GLOBAL_TEMP", (DL_FUNC) &_hector_GLOBAL_TEMP, 0},
     {"_hector_GLOBAL_TEMPEQ", (DL_FUNC) &_hector_GLOBAL_TEMPEQ, 0},
     {"_hector_DIFFUSIVITY", (DL_FUNC) &_hector_DIFFUSIVITY, 0},

src/rcpp_constants.cpp

@@ -936,13 +936,20 @@ String ECS() {
 return D_ECS;
 }
 
-//' @describeIn parameters Aerosol scaling factor (\code{"(unitless)"})
+//' @describeIn parameters Aerosol forcing scaling factor (\code{"(unitless)"})
 //' @export
 // [[Rcpp::export]]
 String AERO_SCALE() {
   return D_AERO_SCALE;
 }
 
+//' @describeIn parameters Volcanic forcing scaling factor (\code{"(unitless)"})
+//' @export
+// [[Rcpp::export]]
+String VOLCANIC_SCALE() {
+return D_VOLCANIC_SCALE;
+}
+
 //' @describeIn temperature Global mean temperature
 //' @export
 // [[Rcpp::export]]

src/temperature_component.cpp

@@ -92,6 +92,7 @@ void TemperatureComponent::init( Core* coreptr ) {
     diff.set( 0.55, U_CM2_S );    // default ocean heat diffusivity, cm2/s. value is CDICE default (varname is kappa there).
     S.set( 3.0, U_DEGC );         // default climate sensitivity, K (varname is t2co in CDICE).
     alpha.set( 1.0, U_UNITLESS);  // default aerosol scaling, unitless (similar to alpha in CDICE).
+    volscl.set(1.0, U_UNITLESS);  // Default volcanic scaling, unitless (works the same way as alpha)
 
     // Register the data we can provide
     core->registerCapability( D_GLOBAL_TEMP, getComponentName() );
@@ -111,6 +112,7 @@ void TemperatureComponent::init( Core* coreptr ) {
     core->registerInput(D_ECS, getComponentName());
     core->registerInput(D_DIFFUSIVITY, getComponentName());
     core->registerInput(D_AERO_SCALE, getComponentName());
+    core->registerInput(D_VOLCANIC_SCALE, getComponentName());
 }
 
 //------------------------------------------------------------------------------
@@ -151,6 +153,9 @@ void TemperatureComponent::setData( const string& varName,
 	} else if( varName == D_AERO_SCALE ) {
             H_ASSERT( data.date == Core::undefinedIndex(), "date not allowed" );
             alpha = data.getUnitval(U_UNITLESS);
+        } else if(varName == D_VOLCANIC_SCALE) {
+            H_ASSERT( data.date == Core::undefinedIndex(), "date not allowed" );
+            volscl = data.getUnitval(U_UNITLESS);
         } else if( varName == D_TGAV_CONSTRAIN ) {
             H_ASSERT( data.date != Core::undefinedIndex(), "date required" );
             tgav_constrain.set(data.date, data.getUnitval(U_DEGC));
@@ -324,8 +329,11 @@ void TemperatureComponent::run( const double runToDate ) throw ( h_exception ) {
     double aero_forcing =
         double(core->sendMessage( M_GETDATA, D_RF_BC ).value( U_W_M2 )) + double(core->sendMessage( M_GETDATA, D_RF_OC).value( U_W_M2 )) + 
         double(core->sendMessage( M_GETDATA, D_RF_SO2d ).value( U_W_M2 )) + double(core->sendMessage( M_GETDATA, D_RF_SO2i ).value( U_W_M2 ));
+    double volcanic_forcing = double(core->sendMessage(M_GETDATA, D_RF_VOL));
 
-    forcing[tstep] = double(core->sendMessage( M_GETDATA, D_RF_TOTAL ).value( U_W_M2 )) - ( 1.0 - alpha ) * aero_forcing;
+    forcing[tstep] = double(core->sendMessage(M_GETDATA, D_RF_TOTAL).value(U_W_M2))
+                      - (1.0 - alpha) * aero_forcing
+                      - (1.0 - volscl) * volcanic_forcing;
 
     // Initialize variables for time-stepping through the model
     double DQ1 = 0.0;
@@ -441,8 +449,8 @@ unitval TemperatureComponent::getData( const std::string& varName,
             returnval = heatflux;
         } else if( varName == D_ECS ) {
             returnval = S;
-        } else if(varName == D_AERO_SCALE) {
-            returnval = alpha;
+        } else if(varName == D_VOLCANIC_SCALE) {
+            returnval = volscl;
         } else {
             H_THROW( "Caller is requesting unknown variable: " + varName ); 
         }

tests/testthat/test_hector.R

@@ -6,13 +6,14 @@ testvars <- c(ATMOSPHERIC_CO2(), RF_TOTAL(), GLOBAL_TEMP())
 dates <- 2000:2300
 
 test_that('Basic hcore functionality works', {
-    hc <- newcore(file.path(inputdir, 'hector_rcp45.ini'), suppresslogging = TRUE)
+    hc <- newcore(file.path(inputdir, 'hector_rcp45.ini'), name='RCP45', suppresslogging = TRUE)
     run(hc, 2100)
     expect_true(inherits(hc, 'hcore'))
     expect_true(isactive(hc))
     expect_equal(startdate(hc), 1745)
     expect_equal(enddate(hc), 2300)
     expect_equal(getdate(hc), 2100)
+    expect_equal(getname(hc), 'RCP45')
 
     expect_error(run(hc, 2050), "is prior to the current date")
     expect_silent(run(hc, 2100))

tests/testthat/test_parameters.R

@@ -102,8 +102,7 @@ test_that('Lowering diffusivity increases temperature', {
     shutdown(hc)
 })
 
-test_that('Lowering aerosol scaling factor increases temperature',
-      {
+test_that('Lowering aerosol forcing scaling factor increases temperature', {
     hc <- newcore(file.path(inputdir, 'hector_rcp45.ini'), suppresslogging =
                     TRUE)
     run(hc, 2100)
@@ -124,3 +123,20 @@ test_that('Lowering aerosol scaling factor increases temperature',
     shutdown(hc)
 })
 
+test_that('Increasing volcanic forcing scaling factor increases the effect of volcanism', {
+    ## Use the difference between 1960 temperature and 1965 temperature as a
+    ## measure of the volcanic effect.
+    hc <- newcore(file.path(inputdir, 'hector_rcp45.ini'), suppresslogging=TRUE)
+    run(hc, 1971)
+    dates <- c(1960, 1965)
+    tbase <- fetchvars(hc, dates, GLOBAL_TEMP())
+    vsibase <- tbase$value[1] - tbase$value[2]
+
+    setvar(hc, NA, VOLCANIC_SCALE(), 1.5, getunits(VOLCANIC_SCALE()))
+    reset(hc)
+    run(hc, 1971)
+    tscl <- fetchvars(hc, dates, GLOBAL_TEMP())
+    vsiscl <- tscl$value[1] - tscl$value[2]
+
+    expect_gt(vsiscl, vsibase)
+})