diff --git a/DESCRIPTION b/DESCRIPTION index d9a2b6877..01f2823ac 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: hector Title: The Hector Simple Climate Model -Version: 2.2.1 +Version: 2.2.2 Authors@R: c(person("Robert", "Link", email = "robert.link@pnnl.gov", role = c("aut", "cre")), person("Corinne", "Hartin", email = "corinne.hartin@pnnl.gov", role = "aut"), person("Ben", "Bond-Lamberty", email = "BondLamberty@pnnl.gov", role = "aut"), @@ -14,7 +14,7 @@ License: GPL-3 Encoding: UTF-8 LazyData: true LinkingTo: Rcpp, BH -Imports: Rcpp (>= 0.12), BH (>= 1.65) +Imports: Rcpp (>= 0.12), BH (>= 1.69) Suggests: testthat, nleqslv, diff --git a/R/RcppExports.R b/R/RcppExports.R index 7670ab519..a6cbfe3f7 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -2,7 +2,7 @@ # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 #' @include aadoc.R -#' @describeIn msgtype Message type for retrieving data from a component +#' @describeIn msgtype Message type for retrieving data from a component #' @export #' @keywords internal GETDATA <- function() { @@ -40,7 +40,7 @@ LL_SEVERE <- function() { .Call('_hector_LL_SEVERE', PACKAGE = 'hector') } -#' @describeIn emissions Black carbon emissions (\code{"Tg"}) +#' @describeIn emissions Black carbon emissions (\code{"Tg"}) #' @export EMISSIONS_BC <- function() { .Call('_hector_EMISSIONS_BC', PACKAGE = 'hector') @@ -118,6 +118,14 @@ RF_VOL <- function() { .Call('_hector_RF_VOL', PACKAGE = 'hector') } +RFADJ_PREFIX <- function() { + .Call('_hector_RFADJ_PREFIX', PACKAGE = 'hector') +} + +RF_PREFIX <- function() { + .Call('_hector_RF_PREFIX', PACKAGE = 'hector') +} + #' @describeIn haloforcings Radiative forcing due to CF4 #' @export RF_CF4 <- function() { diff --git a/R/aadoc.R b/R/aadoc.R index 666fcb4a4..85a8a2a29 100644 --- a/R/aadoc.R +++ b/R/aadoc.R @@ -76,7 +76,9 @@ NULL #' These identifiers specify forcing values that can be provided by hector via #' one of the myriad halocarbon components. All of the values corresponding to #' these identifiers are read-only (\emph{i.e.}, they can only appear in -#' \code{\link{GETDATA}} messages.) +#' \code{\link{GETDATA}} messages.) The forcings returned are the +#' \emph{relative} forcings, with the base year (typically 1750) values +#' subtracted off. #' #' @inheritSection msgtype Note #' diff --git a/R/messages.R b/R/messages.R index d0fe9b788..4a884a55e 100644 --- a/R/messages.R +++ b/R/messages.R @@ -53,6 +53,9 @@ fetchvars <- function(core, dates, vars=NULL, scenario=NULL) lapply(vars, function(v) { sendmessage(core, GETDATA(), v, dates, NA, '') })) + ## Fix the variable name for the adjusted halocarbon forcings so that they are + ## consistent with other forcings. + rslt$variable <- sub(paste0('^',RFADJ_PREFIX()), RF_PREFIX(), rslt$variable) cols <- names(rslt) rslt$scenario <- scenario ## reorder the columns to put the scenario name first diff --git a/changelog.txt b/changelog.txt index cdf11e8cd..caf5f4622 100644 --- a/changelog.txt +++ b/changelog.txt @@ -1,3 +1,19 @@ +Hector 2.2.2 +============ +* Fix bug that was causing requests for H2O forcing in the R interface + to return N2O forcing instead (the model internals were + unaffected). +* Fix bug that was causing API requests for halocarbon forcing to + return absolute forcing values, rather than values relative to the + base year (which is what is done for all other forcings). +* Add missing RF_VOL (volcanic forcing) dependency to forcing component. In practice the + missing dependency had no effect because the forcing component + already had a dependency on the SO2 component, which also provides + the volcanic forcing, but in the event that we ever split them up, + this will ensure that the dependency graph is still valid. + +(All 2.2.2 changes merged in PR #303) + Hector 2.2.1 ============ * Report global averages for land temperature, ocean air diff --git a/inst/include/component_data.hpp b/inst/include/component_data.hpp index 7f3fc08ce..3de3444a3 100644 --- a/inst/include/component_data.hpp +++ b/inst/include/component_data.hpp @@ -82,6 +82,41 @@ #define D_RF_halon2402 D_RF_PREFIX halon2402_COMPONENT_BASE #define D_RF_CH3Cl D_RF_PREFIX CH3Cl_COMPONENT_BASE #define D_RF_CH3Br D_RF_PREFIX CH3Br_COMPONENT_BASE + +// Adjusted (i.e., relative) halocarbon forcings +// Forcings are tracked relative to the base year forcings, but the +// halocarbon components don't know that. These capabilities allow +// a caller to fetch the adjusted values from the forcing component. +#define D_RFADJ_PREFIX "Fadj" +#define D_RFADJ_CF4 D_RFADJ_PREFIX CF4_COMPONENT_BASE +#define D_RFADJ_C2F6 D_RFADJ_PREFIX C2F6_COMPONENT_BASE +#define D_RFADJ_HFC23 D_RFADJ_PREFIX HFC23_COMPONENT_BASE +#define D_RFADJ_HFC32 D_RFADJ_PREFIX HFC32_COMPONENT_BASE +#define D_RFADJ_HFC4310 D_RFADJ_PREFIX HFC4310_COMPONENT_BASE +#define D_RFADJ_HFC125 D_RFADJ_PREFIX HFC125_COMPONENT_BASE +#define D_RFADJ_HFC134a D_RFADJ_PREFIX HFC134a_COMPONENT_BASE +#define D_RFADJ_HFC143a D_RFADJ_PREFIX HFC143a_COMPONENT_BASE +#define D_RFADJ_HFC227ea D_RFADJ_PREFIX HFC227ea_COMPONENT_BASE +#define D_RFADJ_HFC245fa D_RFADJ_PREFIX HFC245fa_COMPONENT_BASE +#define D_RFADJ_SF6 D_RFADJ_PREFIX SF6_COMPONENT_BASE +#define D_RFADJ_CFC11 D_RFADJ_PREFIX CFC11_COMPONENT_BASE +#define D_RFADJ_CFC12 D_RFADJ_PREFIX CFC12_COMPONENT_BASE +#define D_RFADJ_CFC113 D_RFADJ_PREFIX CFC113_COMPONENT_BASE +#define D_RFADJ_CFC114 D_RFADJ_PREFIX CFC114_COMPONENT_BASE +#define D_RFADJ_CFC115 D_RFADJ_PREFIX CFC115_COMPONENT_BASE +#define D_RFADJ_CCl4 D_RFADJ_PREFIX CCl4_COMPONENT_BASE +#define D_RFADJ_CH3CCl3 D_RFADJ_PREFIX CH3CCl3_COMPONENT_BASE +#define D_RFADJ_HCF22 D_RFADJ_PREFIX HCF22_COMPONENT_BASE +#define D_RFADJ_HCF141b D_RFADJ_PREFIX HCF141b_COMPONENT_BASE +#define D_RFADJ_HCF142b D_RFADJ_PREFIX HCF142b_COMPONENT_BASE +#define D_RFADJ_halon1211 D_RFADJ_PREFIX halon1211_COMPONENT_BASE +#define D_RFADJ_halon1301 D_RFADJ_PREFIX halon1301_COMPONENT_BASE +#define D_RFADJ_halon2402 D_RFADJ_PREFIX halon2402_COMPONENT_BASE +#define D_RFADJ_CH3Cl D_RFADJ_PREFIX CH3Cl_COMPONENT_BASE +#define D_RFADJ_CH3Br D_RFADJ_PREFIX CH3Br_COMPONENT_BASE +#define N_HALO_FORCINGS 26 + + // halocarbon emissions #define D_EMISSIONS_CF4 CF4_COMPONENT_BASE EMISSIONS_EXTENSION #define D_EMISSIONS_C2F6 C2F6_COMPONENT_BASE EMISSIONS_EXTENSION diff --git a/inst/include/forcing_component.hpp b/inst/include/forcing_component.hpp index 4fe489415..4e9fb922f 100644 --- a/inst/include/forcing_component.hpp +++ b/inst/include/forcing_component.hpp @@ -89,6 +89,10 @@ class ForcingComponent : public IModelComponent { Core* core; //! Core Logger logger; //! Logger + + static const char *adjusted_halo_forcings[]; //! Capability strings for halocarbon forcings + static const char *halo_forcing_names[]; //! Internal names of halocarbon forcings + static std::map forcing_name_map; }; } diff --git a/inst/include/h_util.hpp b/inst/include/h_util.hpp index 7152a53cc..0f938d36c 100644 --- a/inst/include/h_util.hpp +++ b/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.2.1" +#define MODEL_VERSION "2.2.2" #define OUTPUT_DIRECTORY "output/" diff --git a/man/haloforcings.Rd b/man/haloforcings.Rd index b5f1b2259..61f4302b1 100644 --- a/man/haloforcings.Rd +++ b/man/haloforcings.Rd @@ -86,7 +86,9 @@ RF_CH3BR() These identifiers specify forcing values that can be provided by hector via one of the myriad halocarbon components. All of the values corresponding to these identifiers are read-only (\emph{i.e.}, they can only appear in -\code{\link{GETDATA}} messages.) +\code{\link{GETDATA}} messages.) The forcings returned are the +\emph{adjusted} forcings, with the base year (typically 1750) values +subtracted off. } \section{Functions}{ \itemize{ diff --git a/src/Makevars b/src/Makevars index db5aa825b..9e1740bad 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1 +1,2 @@ -PKG_CPPFLAGS = -I../inst/include -DUSE_RCPP +CXX_STD = CXX11 +PKG_CPPFLAGS = -I../inst/include -DUSE_RCPP diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 4105f31f7..0ae948efa 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -195,6 +195,26 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// RFADJ_PREFIX +String RFADJ_PREFIX(); +RcppExport SEXP _hector_RFADJ_PREFIX() { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = Rcpp::wrap(RFADJ_PREFIX()); + return rcpp_result_gen; +END_RCPP +} +// RF_PREFIX +String RF_PREFIX(); +RcppExport SEXP _hector_RF_PREFIX() { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = Rcpp::wrap(RF_PREFIX()); + return rcpp_result_gen; +END_RCPP +} // RF_CF4 String RF_CF4(); RcppExport SEXP _hector_RF_CF4() { @@ -1493,6 +1513,8 @@ static const R_CallMethodDef CallEntries[] = { {"_hector_RF_SO2I", (DL_FUNC) &_hector_RF_SO2I, 0}, {"_hector_RF_SO2", (DL_FUNC) &_hector_RF_SO2, 0}, {"_hector_RF_VOL", (DL_FUNC) &_hector_RF_VOL, 0}, + {"_hector_RFADJ_PREFIX", (DL_FUNC) &_hector_RFADJ_PREFIX, 0}, + {"_hector_RF_PREFIX", (DL_FUNC) &_hector_RF_PREFIX, 0}, {"_hector_RF_CF4", (DL_FUNC) &_hector_RF_CF4, 0}, {"_hector_RF_C2F6", (DL_FUNC) &_hector_RF_C2F6, 0}, {"_hector_RF_HFC23", (DL_FUNC) &_hector_RF_HFC23, 0}, diff --git a/src/core.cpp b/src/core.cpp index ed7f2e477..87c8923d6 100644 --- a/src/core.cpp +++ b/src/core.cpp @@ -47,6 +47,7 @@ using namespace std; * \sa init() */ Core::Core(Logger::LogLevel loglvl, bool echotoscreen, bool echotofile) : + setup_complete(false), run_name( "" ), startDate( -1.0 ), endDate( -1.0 ), @@ -54,7 +55,6 @@ Core::Core(Logger::LogLevel loglvl, bool echotoscreen, bool echotofile) : isInited( false ), do_spinup( true ), max_spinup( 2000 ), - setup_complete(false), in_spinup( false ) { glog.open(string(MODEL_NAME), echotoscreen, echotofile, loglvl); @@ -565,6 +565,8 @@ void Core::registerCapability(const string& capabilityName, const string& compon // Only add the capability if it doesn't already exist. // Adding a duplicate capability has no useful effect anyhow. componentCapabilities.insert( pair( capabilityName, componentName ) ); + H_LOG(glog, Logger::DEBUG) << capabilityName << " registered to component " << componentName << "\n"; + } } diff --git a/src/forcing_component.cpp b/src/forcing_component.cpp index ab3d2351e..16524cd22 100644 --- a/src/forcing_component.cpp +++ b/src/forcing_component.cpp @@ -23,7 +23,86 @@ #include "avisitor.hpp" namespace Hector { - + +/* These next two arrays and the map that connects them are a + * workaround for the problems created by storing the halocarbon + * forcings in the halocarbon components. Because the halocarbon + * components don't know about the base year adjustments, they can't + * provide the forcings relative to the base year, which is what + * outside callers will generally want. Internally, however, we still + * need to be able to get the raw forcings from the halocarbon + * components, so we can't just change everything to point at the + * forcing component (which would return the base year adjusted value). + * + * The solution we adopted was to create a second set of capabilities + * to refer to the adjusted values, and we let the forcing component + * intercept those. However, the forcing values themselves are stored + * under the names used for the unadjusted values, so we have to have + * a name translation table so that we can find the data that the + * message is asking for. In the end, the whole process winds up + * being a little ugly, but it gets the job done. + */ + +const char *ForcingComponent::adjusted_halo_forcings[N_HALO_FORCINGS] = { + D_RFADJ_CF4, + D_RFADJ_C2F6, + D_RFADJ_HFC23, + D_RFADJ_HFC32, + D_RFADJ_HFC4310, + D_RFADJ_HFC125, + D_RFADJ_HFC134a, + D_RFADJ_HFC143a, + D_RFADJ_HFC227ea, + D_RFADJ_HFC245fa, + D_RFADJ_SF6, + D_RFADJ_CFC11, + D_RFADJ_CFC12, + D_RFADJ_CFC113, + D_RFADJ_CFC114, + D_RFADJ_CFC115, + D_RFADJ_CCl4, + D_RFADJ_CH3CCl3, + D_RFADJ_HCF22, + D_RFADJ_HCF141b, + D_RFADJ_HCF142b, + D_RFADJ_halon1211, + D_RFADJ_halon1301, + D_RFADJ_halon2402, + D_RFADJ_CH3Cl, + D_RFADJ_CH3Br +}; + +const char *ForcingComponent::halo_forcing_names[N_HALO_FORCINGS] = { + D_RF_CF4, + D_RF_C2F6, + D_RF_HFC23, + D_RF_HFC32, + D_RF_HFC4310, + D_RF_HFC125, + D_RF_HFC134a, + D_RF_HFC143a, + D_RF_HFC227ea, + D_RF_HFC245fa, + D_RF_SF6, + D_RF_CFC11, + D_RF_CFC12, + D_RF_CFC113, + D_RF_CFC114, + D_RF_CFC115, + D_RF_CCl4, + D_RF_CH3CCl3, + D_RF_HCF22, + D_RF_HCF141b, + D_RF_HCF142b, + D_RF_halon1211, + D_RF_halon1301, + D_RF_halon2402, + D_RF_CH3Cl, + D_RF_CH3Br +}; + +std::map ForcingComponent::forcing_name_map; + using namespace std; //------------------------------------------------------------------------------ @@ -74,6 +153,10 @@ void ForcingComponent::init( Core* coreptr ) { core->registerCapability( D_RF_SO2i, getComponentName()); core->registerCapability( D_RF_SO2, getComponentName()); core->registerCapability( D_RF_VOL, getComponentName()); + for(int i=0; iregisterCapability(adjusted_halo_forcings[i], getComponentName()); + forcing_name_map[adjusted_halo_forcings[i]] = halo_forcing_names[i]; + } // Register our dependencies @@ -407,7 +490,15 @@ unitval ForcingComponent::getData( const std::string& varName, } } } else { - std::map::const_iterator forcing = forcings.find( varName ); + std::string forcing_name; + auto forcit = forcing_name_map.find(varName); + if(forcit != forcing_name_map.end()) { + forcing_name = forcing_name_map[varName]; + } + else { + forcing_name = varName; + } + std::map::const_iterator forcing = forcings.find(forcing_name); if ( forcing != forcings.end() ) { // from the forcing map returnval = forcing->second; diff --git a/src/makefile.standalone b/src/makefile.standalone index a23bc28f1..6278959c8 100644 --- a/src/makefile.standalone +++ b/src/makefile.standalone @@ -1,10 +1,10 @@ ifeq ($(strip $(CXX)),) CXX = g++ endif -CXXFLAGS = -g $(INCLUDES) $(OPTFLAGS) $(CXXEXTRA) $(CXXPROF) $(WFLAGS) -MMD +CXXFLAGS = -g $(INCLUDES) $(OPTFLAGS) $(CXXEXTRA) $(CXXPROF) $(WFLAGS) -MMD -std=c++11 CFLAGS = -g $(INCLUDES) $(OPTFLAGS) -MMD INCLUDES = -I$(BOOSTROOT) -I$(HDRDIR) -WFLAGS = -Wall -Wno-unused-local-typedefs -Wno-c++11-extensions # Turn on warnings, turn off one particularly annoying one that infests Boost libs, also turn off warning about c++-11 extensions. +WFLAGS = -Wall -Wno-unused-local-typedefs # Turn on warnings, turn off one particularly annoying one that infests Boost libs OPTFLAGS = -O3 LDFLAGS = $(CXXPROF) -L$(BOOSTLIB) -L. -Wl,-rpath $(BOOSTLIB) diff --git a/src/rcpp_constants.cpp b/src/rcpp_constants.cpp index d4442f26c..f6abdb036 100644 --- a/src/rcpp_constants.cpp +++ b/src/rcpp_constants.cpp @@ -8,7 +8,7 @@ using namespace Rcpp; * Functions to expose constants defined in C++ headers to R * * The only way I've been able to find to get these constants into R is to provide - * Rcpp exported functions that return them. + * Rcpp exported functions that return them. *********************************************************************************/ // [[Rcpp::plugins("cpp11")]] @@ -17,7 +17,7 @@ using namespace Rcpp; * Message types */ //' @include aadoc.R -//' @describeIn msgtype Message type for retrieving data from a component +//' @describeIn msgtype Message type for retrieving data from a component //' @export //' @keywords internal // [[Rcpp::export]] @@ -69,7 +69,7 @@ return (int) Hector::Logger::SEVERE; *****************************************************************/ /* BC component */ -//' @describeIn emissions Black carbon emissions (\code{"Tg"}) +//' @describeIn emissions Black carbon emissions (\code{"Tg"}) //' @export // [[Rcpp::export]] String EMISSIONS_BC() { @@ -109,7 +109,7 @@ return D_RF_N2O; //' @export // [[Rcpp::export]] String RF_H2O() { -return D_RF_N2O; +return D_RF_H2O; } //' @describeIn forcings Radiative forcing due to ozone @@ -170,186 +170,201 @@ return D_RF_VOL; // return D_RF_halocarbons; // } +// This is the prefix for adjusted radiative forcing. We use it +// internally to correct the reported names of the halocarbon +// forcings. +// [[Rcpp::export]] +String RFADJ_PREFIX() { +return D_RFADJ_PREFIX; +} + +// Prefix for unadjusted forcing. This is what we substitute for +// the adjusted forcing prefix +// [[Rcpp::export]] +String RF_PREFIX() { +return D_RF_PREFIX; +} + //' @describeIn haloforcings Radiative forcing due to CF4 //' @export // [[Rcpp::export]] String RF_CF4() { -return D_RF_CF4; +return D_RFADJ_CF4; } //' @describeIn haloforcings Radiative forcing due to C2F6 //' @export // [[Rcpp::export]] String RF_C2F6() { -return D_RF_C2F6; +return D_RFADJ_C2F6; } //' @describeIn haloforcings Radiative forcing due to HFC-23 //' @export // [[Rcpp::export]] String RF_HFC23() { -return D_RF_HFC23; +return D_RFADJ_HFC23; } //' @describeIn haloforcings Radiative forcing due to HFC-32 //' @export // [[Rcpp::export]] String RF_HFC32() { -return D_RF_HFC32; +return D_RFADJ_HFC32; } //' @describeIn haloforcings Radiative forcing due to HFC-4310 //' @export // [[Rcpp::export]] String RF_HFC4310() { -return D_RF_HFC4310; +return D_RFADJ_HFC4310; } //' @describeIn haloforcings Radiative forcing due to HFC-125 //' @export // [[Rcpp::export]] String RF_HFC125() { -return D_RF_HFC125; +return D_RFADJ_HFC125; } //' @describeIn haloforcings Radiative forcing due to HFC-134a //' @export // [[Rcpp::export]] String RF_HFC134A() { -return D_RF_HFC134a; +return D_RFADJ_HFC134a; } //' @describeIn haloforcings Radiative forcing due to HFC-143a //' @export // [[Rcpp::export]] String RF_HFC143A() { -return D_RF_HFC143a; +return D_RFADJ_HFC143a; } //' @describeIn haloforcings Radiative forcing due to HFC227ea //' @export // [[Rcpp::export]] String RF_HFC227EA() { -return D_RF_HFC227ea; +return D_RFADJ_HFC227ea; } //' @describeIn haloforcings Radiative forcing due to HFC-245fa //' @export // [[Rcpp::export]] String RF_HFC245FA() { -return D_RF_HFC245fa; +return D_RFADJ_HFC245fa; } //' @describeIn haloforcings Radiative forcing due to sulfur hexafluoride //' @export // [[Rcpp::export]] String RF_SF6() { -return D_RF_SF6; +return D_RFADJ_SF6; } //' @describeIn haloforcings Radiative forcing due to CFC-11 //' @export // [[Rcpp::export]] String RF_CFC11() { -return D_RF_CFC11; +return D_RFADJ_CFC11; } //' @describeIn haloforcings Radiative forcing due to CFC-12 //' @export // [[Rcpp::export]] String RF_CFC12() { -return D_RF_CFC12; +return D_RFADJ_CFC12; } //' @describeIn haloforcings Radiative forcing due to CFC-113 //' @export // [[Rcpp::export]] String RF_CFC113() { -return D_RF_CFC113; +return D_RFADJ_CFC113; } //' @describeIn haloforcings Radiative forcing due to CFC-114 //' @export // [[Rcpp::export]] String RF_CFC114() { -return D_RF_CFC114; +return D_RFADJ_CFC114; } //' @describeIn haloforcings Radiative forcing due to CFC-115 //' @export // [[Rcpp::export]] String RF_CFC115() { -return D_RF_CFC115; +return D_RFADJ_CFC115; } //' @describeIn haloforcings Radiative forcing due to carbon tetrachloride //' @export // [[Rcpp::export]] String RF_CCL4() { -return D_RF_CCl4; +return D_RFADJ_CCl4; } //' @describeIn haloforcings Radiative forcing due to trichloroethane //' @export // [[Rcpp::export]] String RF_CH3CCL3() { -return D_RF_CH3CCl3; +return D_RFADJ_CH3CCl3; } //' @describeIn haloforcings Radiative forcing due to HCFC-22 //' @export // [[Rcpp::export]] String RF_HCF22() { -return D_RF_HCF22; +return D_RFADJ_HCF22; } //' @describeIn haloforcings Radiative forcing due to HCFC-141b //' @export // [[Rcpp::export]] String RF_HCF141B() { -return D_RF_HCF141b; +return D_RFADJ_HCF141b; } //' @describeIn haloforcings Radiative forcing due to HCFC-142b //' @export // [[Rcpp::export]] String RF_HCF142B() { -return D_RF_HCF142b; +return D_RFADJ_HCF142b; } //' @describeIn haloforcings Radiative forcing due to halon-1211 //' @export // [[Rcpp::export]] String RF_HALON1211() { -return D_RF_halon1211; +return D_RFADJ_halon1211; } //' @describeIn haloforcings Radiative forcing due to halon-1301 //' @export // [[Rcpp::export]] String RF_HALON1301() { -return D_RF_halon1301; +return D_RFADJ_halon1301; } //' @describeIn haloforcings Radiative forcing due to halon-2402 //' @export // [[Rcpp::export]] String RF_HALON2402() { -return D_RF_halon2402; +return D_RFADJ_halon2402; } //' @describeIn haloforcings Radiative forcing due to chloromethane //' @export // [[Rcpp::export]] String RF_CH3CL() { -return D_RF_CH3Cl; +return D_RFADJ_CH3Cl; } //' @describeIn haloforcings Radiative forcing due to bromomethane //' @export // [[Rcpp::export]] String RF_CH3BR() { -return D_RF_CH3Br; +return D_RFADJ_CH3Br; } /* halocarbon emissions */ @@ -533,7 +548,7 @@ return D_EMISSIONS_CH3Cl; // [[Rcpp::export]] String EMISSIONS_CH3BR() { return D_EMISSIONS_CH3Br; -} +} /* Methane component */ diff --git a/src/rcpp_hector.cpp b/src/rcpp_hector.cpp index 7a133d778..0c6556810 100644 --- a/src/rcpp_hector.cpp +++ b/src/rcpp_hector.cpp @@ -296,7 +296,9 @@ DataFrame sendmessage(Environment core, String msgtype, String capability, Numer // Assemble a data frame with the results: date, var, value, units DataFrame result = DataFrame::create(Named("year")=date, Named("variable")=capability, - Named("value")=valueout, Named("units")=unitsout); + Named("value")=valueout, + Named("units")=unitsout, + Named("stringsAsFactors")=false); return result; } diff --git a/src/temperature_component.cpp b/src/temperature_component.cpp index 53e9a2fad..b7757f888 100644 --- a/src/temperature_component.cpp +++ b/src/temperature_component.cpp @@ -24,7 +24,7 @@ #include "avisitor.hpp" namespace Hector { - + using namespace std; @@ -39,12 +39,12 @@ TemperatureComponent::TemperatureComponent() { */ TemperatureComponent::~TemperatureComponent(){ } - + //------------------------------------------------------------------------------ // documentation is inherited string TemperatureComponent::getComponentName() const { const string name = TEMPERATURE_COMPONENT_NAME; - + return name; } @@ -58,7 +58,7 @@ string TemperatureComponent::getComponentName() const { */ void TemperatureComponent::invert_1d_2x2_matrix(double * x, double * y) { double temp_d = (x[0]*x[3] - x[1]*x[2]); - + if(temp_d == 0) { H_THROW("Temperature: Matrix inversion divide by zero."); } @@ -67,33 +67,33 @@ void TemperatureComponent::invert_1d_2x2_matrix(double * x, double * y) { y[1] = temp * -1 * x[1]; y[2] = temp * -1 * x[2]; y[3] = temp * x[0]; - + return; } - + //------------------------------------------------------------------------------ // documentation is inherited void TemperatureComponent::init( Core* coreptr ) { logger.open( getComponentName(), false, coreptr->getGlobalLogger().getEchoToFile(), coreptr->getGlobalLogger().getMinLogLevel() ); H_LOG( logger, Logger::DEBUG ) << "hello " << getComponentName() << std::endl; - + tgaveq.set( 0.0, U_DEGC, 0.0 ); tgav.set( 0.0, U_DEGC, 0.0 ); flux_mixed.set( 0.0, U_W_M2, 0.0 ); flux_interior.set( 0.0, U_W_M2, 0.0 ); heatflux.set( 0.0, U_W_M2, 0.0 ); - + core = coreptr; - + tgav_constrain.allowInterp( true ); tgav_constrain.name = D_TGAV_CONSTRAIN; - + // Define the doeclim parameters 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() ); core->registerCapability( D_LAND_TEMP, getComponentName() ); @@ -110,6 +110,7 @@ void TemperatureComponent::init( Core* coreptr ) { core->registerDependency( D_RF_OC, getComponentName() ); core->registerDependency( D_RF_SO2d, getComponentName() ); core->registerDependency( D_RF_SO2i, getComponentName() ); + core->registerDependency( D_RF_VOL, getComponentName() ); // Register the inputs we can receive from outside core->registerInput(D_ECS, getComponentName()); @@ -125,15 +126,15 @@ unitval TemperatureComponent::sendMessage( const std::string& message, const message_data info ) throw ( h_exception ) { unitval returnval; - + if( message==M_GETDATA ) { //! Caller is requesting data - return getData( datum, info.date ); + return getData( datum, info.date ); } else if( message==M_SETDATA ) { //! Caller is requesting to set data setData(datum, info); } else { //! We don't handle any other messages H_THROW( "Caller sent unknown message: "+message ); } - + return returnval; } @@ -143,9 +144,9 @@ void TemperatureComponent::setData( const string& varName, const message_data& data ) throw ( h_exception ) { using namespace boost; - + H_LOG( logger, Logger::DEBUG ) << "Setting " << varName << "[" << data.date << "]=" << data.value_str << std::endl; - + try { if( varName == D_ECS ) { H_ASSERT( data.date == Core::undefinedIndex() , "date not allowed" ); @@ -178,17 +179,17 @@ void TemperatureComponent::setData( const string& varName, // documentation is inherited // TO DO: should we put these in the ini file instead? void TemperatureComponent::prepareToRun() throw ( h_exception ) { - + H_LOG( logger, Logger::DEBUG ) << "prepareToRun " << std::endl; - + if( tgav_constrain.size() ) { Logger& glog = core->getGlobalLogger(); H_LOG( glog, Logger::WARNING ) << "Temperature will be overwritten by user-supplied values!" << std::endl; } - + // Initializing all model components that depend on the number of timesteps (ns) ns = core->getEndDate() - core->getStartDate() + 1; - + KT0 = std::vector(ns, 0.0); KTA1 = std::vector(ns, 0.0); KTB1 = std::vector(ns, 0.0); @@ -196,9 +197,9 @@ void TemperatureComponent::prepareToRun() throw ( h_exception ) { KTB2 = std::vector(ns, 0.0); KTA3 = std::vector(ns, 0.0); KTB3 = std::vector(ns, 0.0); - + Ker.resize(ns); - + temp.resize(ns); temp_landair.resize(ns); temp_sst.resize(ns); @@ -207,12 +208,12 @@ void TemperatureComponent::prepareToRun() throw ( h_exception ) { heat_mixed.resize(ns); heat_interior.resize(ns); forcing.resize(ns); - + for(int i=0; i<3; i++) { B[i] = 0.0; C[i] = 0.0; } - + // DOECLIM parameters calculated from constants set in header ocean_area = (1.0 - flnd) * earth_area; // m2 cnum = rlam * flnd + bsi * (1.0 - flnd); // factor from sea-surface climate sensitivity to global mean @@ -228,48 +229,48 @@ void TemperatureComponent::prepareToRun() throw ( h_exception ) { taudif = pow(cas,2) / pow(csw,2) * M_PI / keff; // interior ocean heat uptake time scale, yr tauksl = (1.0 - flnd) * cas / kls; // sea-land heat exchange time scale, yr taukls = flnd * cal / kls; // land-sea heat exchange time scale, yr - + // Components of the analytical solution to the integral found in the temperature difference equation // Third order bottom correction terms will be "more than sufficient" for simulations out to 2500 // (Equation A.25, EK05, or 2.3.23, TK07) - + // First order KT0[ns-1] = 4.0 - 2.0 * pow(2.0, 0.5); KTA1[ns-1] = -8.0 * exp(-taubot / double(dt)) + 4.0 * pow(2.0, 0.5) * exp(-0.5 * taubot / double(dt)); KTB1[ns-1] = 4.0 * pow((M_PI * taubot / double(dt)), 0.5) * (1.0 + erf(pow(0.5 * taubot / double(dt), 0.5)) - 2.0 * erf(pow(taubot / double(dt), 0.5))); - + // Second order KTA2[ns-1] = 8.0 * exp(-4.0 * taubot / double(dt)) - 4.0 * pow(2.0, 0.5) * exp(-2.0 * taubot / double(dt)); KTB2[ns-1] = -8.0 * pow((M_PI * taubot / double(dt)), 0.5) * (1.0 + erf(pow((2.0 * taubot / double(dt)), 0.5)) - 2.0 * erf(2.0 * pow((taubot / double(dt)), 0.5)) ); - + // Third order KTA3[ns-1] = -8.0 * exp(-9.0 * taubot / double(dt)) + 4.0 * pow(2.0, 0.5) * exp(-4.5 * taubot / double(dt)); KTB3[ns-1] = 12.0 * pow((M_PI * taubot / double(dt)), 0.5) * (1.0 + erf(pow((4.5 * taubot / double(dt)), 0.5)) - 2.0 * erf(3.0 * pow((taubot / double(dt)), 0.5)) ); - + // Calculate the kernel component vectors for(int i=0; i<(ns-1); i++) { - + // First order KT0[i] = 4.0 * pow((double(ns-i)), 0.5) - 2.0 * pow((double(ns+1-i)), 0.5) - 2.0 * pow(double(ns-1-i), 0.5); KTA1[i] = -8.0 * pow(double(ns-i), 0.5) * exp(-taubot / double(dt) / double(ns-i)) + 4.0 * pow(double(ns+1-i), 0.5) * exp(-taubot / double(dt) / double(ns+1-i)) + 4.0 * pow(double(ns-1-i), 0.5) * exp(-taubot/double(dt) / double(ns-1-i)); KTB1[i] = 4.0 * pow((M_PI * taubot / double(dt)), 0.5) * ( erf(pow((taubot / double(dt) / double(ns-1-i)), 0.5)) + erf(pow((taubot / double(dt) / double(ns+1-i)), 0.5)) - 2.0 * erf(pow((taubot / double(dt) / double(ns-i)), 0.5)) ); - + // Second order KTA2[i] = 8.0 * pow(double(ns-i), 0.5) * exp(-4.0 * taubot / double(dt) / double(ns-i)) - 4.0 * pow(double(ns+1-i), 0.5) * exp(-4.0 * taubot / double(dt) / double(ns+1-i)) - 4.0 * pow(double(ns-1-i), 0.5) * exp(-4.0 * taubot / double(dt) / double(ns-1-i)); KTB2[i] = -8.0 * pow((M_PI * taubot / double(dt)), 0.5) * ( erf(2.0 * pow((taubot / double(dt) / double(ns-1-i)), 0.5)) + erf(2.0 * pow((taubot / double(dt) / double(ns+1-i)), 0.5)) - 2.0 * erf(2.0 * pow((taubot / double(dt) / double(ns-i)), 0.5)) ); - + // Third order KTA3[i] = -8.0 * pow(double(ns-i), 0.5) * exp(-9.0 * taubot / double(dt) / double(ns-i)) + 4.0 * pow(double(ns+1-i), 0.5) * exp(-9.0 * taubot / double(dt) / double(ns+1-i)) + 4.0 * pow(double(ns-1-i), 0.5) * exp(-9.0 * taubot / double(dt) / double(ns-1-i)); KTB3[i] = 12.0 * pow((M_PI * taubot / double(dt)), 0.5) * ( erf(3.0 * pow((taubot / double(dt) / double(ns-1-i)), 0.5)) + erf(3.0 * pow((taubot / double(dt) / double(ns+1-i)), 0.5)) - 2.0 * erf(3.0 * pow((taubot / double(dt) / double(ns-i)), 0.5)) ); } - + // Sum up the kernel components for(int i=0; igetStartDate(); 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_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 - (1.0 - volscl) * volcanic_forcing; - + // Initialize variables for time-stepping through the model double DQ1 = 0.0; double DQ2 = 0.0; @@ -349,7 +350,7 @@ void TemperatureComponent::run( const double runToDate ) throw ( h_exception ) { double DPAST2 = 0.0; double DTEAUX1 = 0.0; double DTEAUX2 = 0.0; - + // Reset the endogenous varibales for this time step temp[tstep] = 0.0; temp_landair[tstep] = 0.0; @@ -358,22 +359,22 @@ void TemperatureComponent::run( const double runToDate ) throw ( h_exception ) { heat_interior[tstep] = 0.0; heatflux_mixed[tstep] = 0.0; heatflux_interior[tstep] = 0.0; - + // Assume land and ocean forcings are equal to global forcing std::vector QL = forcing; std::vector QO = forcing; - + if (tstep > 0) { - + DelQL = QL[tstep] - QL[tstep - 1]; DelQO = QO[tstep] - QO[tstep - 1]; - + // Assume linear forcing change between tstep and tstep+1 QC1 = (DelQL/cal*(1.0/taucfl+1.0/taukls)-bsi*DelQO/cas/taukls); QC2 = (DelQO/cas*(1.0/taucfs+bsi/tauksl)-DelQL/cal/tauksl); QC1 = QC1 * pow(double(dt), 2.0)/12.0; QC2 = QC2 * pow(double(dt), 2.0)/12.0; - + // ----------------- Initial Conditions -------------------- // Initialization of temperature and forcing vector: // Factor 1/2 in front of Q in Equation A.27, EK05, and Equation 2.3.27, TK07 is a typo! @@ -382,17 +383,17 @@ void TemperatureComponent::run( const double runToDate ) throw ( h_exception ) { DQ2 = 0.5*double(dt)/cas*(QO[tstep]+QO[tstep-1]); DQ1 = DQ1 + QC1; DQ2 = DQ2 + QC2; - + // ---------- SOLVE MODEL ------------------ // Calculate temperatures for(int i = 0; i <= tstep; i++) { DPAST2 = DPAST2 + temp_sst[i] * Ker[ns-tstep+i-1]; } DPAST2 = DPAST2 * fso * pow((double(dt)/taudif), 0.5); - + DTEAUX1 = A[0] * temp_landair[tstep-1] + A[1] * temp_sst[tstep-1]; DTEAUX2 = A[2] * temp_landair[tstep-1] + A[3] * temp_sst[tstep-1]; - + temp_landair[tstep] = IB[0] * (DQ1 + DPAST1 + DTEAUX1) + IB[1] * (DQ2 + DPAST2 + DTEAUX2); temp_sst[tstep] = IB[2] * (DQ1 + DPAST1 + DTEAUX1) + IB[3] * (DQ2 + DPAST2 + DTEAUX2); } @@ -401,7 +402,7 @@ void TemperatureComponent::run( const double runToDate ) throw ( h_exception ) { temp_sst[0] = 0.0; } temp[tstep] = flnd * temp_landair[tstep] + (1.0 - flnd) * bsi * temp_sst[tstep]; - + // Calculate ocean heat uptake [W/m^2] // heatflux[tstep] captures in the heat flux in the period between tstep-1 and tstep. // Numerical implementation of Equation 2.7, EK05, or Equation 2.3.13, TK07) @@ -415,7 +416,7 @@ void TemperatureComponent::run( const double runToDate ) throw ( h_exception ) { heat_mixed[tstep] = heat_mixed[tstep-1] + heatflux_mixed[tstep] * (powtoheat*dt); heat_interior[tstep] = heat_interior[tstep-1] + heatflux_interior[tstep] * (fso*powtoheat*dt); } - + else { // Handle the initial conditions heatflux_mixed[0] = 0.0; heatflux_interior[0] = 0.0; @@ -431,11 +432,11 @@ void TemperatureComponent::run( const double runToDate ) throw ( h_exception ) { // documentation is inherited unitval TemperatureComponent::getData( const std::string& varName, const double date ) throw ( h_exception ) { - + unitval returnval; - + if(date == Core::undefinedIndex()) { - // If no date is supplied, return the current value + // If no date is supplied, return the current value if( varName == D_GLOBAL_TEMP ) { returnval = tgav; } else if( varName == D_GLOBAL_TEMPEQ ) { @@ -461,7 +462,7 @@ unitval TemperatureComponent::getData( const std::string& varName, } else if(varName == D_VOLCANIC_SCALE) { returnval = volscl; } else { - H_THROW( "Caller is requesting unknown variable: " + varName ); + H_THROW( "Caller is requesting unknown variable: " + varName ); } } else { @@ -491,7 +492,7 @@ unitval TemperatureComponent::getData( const std::string& varName, returnval = unitval(value, U_W_M2); } } - + return returnval; } @@ -504,7 +505,7 @@ void TemperatureComponent::reset(double time) throw(h_exception) // vectors. if(time < core->getStartDate()) // in this case, reset to the starting value. time = core->getStartDate(); - + int tstep = time - core->getStartDate(); setoutputs(tstep); H_LOG(logger, Logger::NOTICE) @@ -539,7 +540,7 @@ void TemperatureComponent::setoutputs(int tstep) tgav_land.set(temp_landair[tstep], U_DEGC, 0.0); tgav_sst.set(temp_sst[tstep], U_DEGC, 0.0); temp_oceanair = bsi * temp_sst[tstep]; - tgav_oceanair.set(temp_oceanair, U_DEGC, 0.0); + tgav_oceanair.set(temp_oceanair, U_DEGC, 0.0); } } diff --git a/tests/testthat/test_hector.R b/tests/testthat/test_hector.R index ff42afa9c..0c92de20f 100644 --- a/tests/testthat/test_hector.R +++ b/tests/testthat/test_hector.R @@ -185,3 +185,35 @@ test_that("Setting past or parameter values does trigger a reset.", { shutdown(hc) }) + + +test_that('Test RF output.', { + + hc <- newcore(file.path(inputdir, 'hector_rcp45.ini'), suppresslogging = TRUE) + run(hc, 2100) + + # Extract the total cliamte RF. + total_rf <- fetchvars(hc, dates = 1750:2100, RF_TOTAL()) + + # The vector of all the individual componets that contribute to RF. + rf_componets <- c(RF_BC(), RF_C2F6(), RF_CCL4(), RF_CF4(), RF_CFC11(), RF_CFC113(), RF_CFC114(), RF_CFC115(), RF_CH3BR(), RF_CH3CCL3(), + RF_CH3CL(), RF_CH4(), RF_CO2(), RF_H2O(), RF_HALON1211(), RF_CFC12(), RF_HALON1301(), RF_HALON2402(), RF_HCF141B(), + RF_HCF142B(), RF_HCF22(), RF_HFC125(), RF_HFC134A(), RF_HFC143A(), RF_HFC227EA(), RF_HFC23(), RF_T_ALBEDO(), + RF_HFC245FA(), RF_HFC32(), RF_HFC4310(), RF_N2O(), RF_O3(), RF_OC(), RF_SF6(), RF_SO2D(), RF_SO2I(), RF_VOL()) + individual_rf <- fetchvars(hc, 1750:2100, rf_componets) + + # The RF value should be equal to 0 in the base year (1750) for all of the RF agents. + values_1750 <- dplyr::filter(individual_rf, year == 1750) + expect_equal(values_1750$value, expected = rep(0, nrow(values_1750)), info = 'RF values in the base year must be 0') + + # That the sum of the RF agents should equal the total climate RF. + individual_rf %>% + dplyr::group_by(year) %>% + dplyr::summarise(value = sum(value)) %>% + dplyr::ungroup() -> + rf_aggregate + + expect_equal(rf_aggregate$value, total_rf$value) + + shutdown(hc) +}) diff --git a/tests/testthat/test_scenarios.R b/tests/testthat/test_scenarios.R index 6ef409b34..e3e38f546 100644 --- a/tests/testthat/test_scenarios.R +++ b/tests/testthat/test_scenarios.R @@ -20,7 +20,7 @@ test_that('RCP scenarios are correct', { ## Get the comparison data sampleoutfile <- sprintf('sample_outputstream_rcp%s.csv', rcp) - sampledata <- read.csv(file.path(sampledir, sampleoutfile), comment.char = '#') + sampledata <- read.csv(file.path(sampledir, sampleoutfile), comment.char = '#', stringsAsFactors = FALSE) sampledata$scenario <- as.character(sampledata$run_name) samplekeep <- sampledata$variable %in% testvars & sampledata$year %in% dates sampledata <- sampledata[samplekeep, c('scenario', 'year','variable','value', 'units')]