Add volcanic scaling factor

NAMESPACE

@@ -129,6 +129,7 @@ export(TID)
 export(TT)
 export(TU)
 export(TWI)
+export(VOLCANIC_SCALE)
 export(VOLCANIC_SO2)
 export(Y2000_SO2)
 export(enddate)

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() {

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/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_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)
+})