Merge pull request #303 from JGCRI/krd_RF_H2O

Fix bugs in reporting for halocarbon forcings.

Version 2.2.2

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,

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

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
 #'

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

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

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

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<std::string, std::string> forcing_name_map;
 };
 
 }

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/"
 

src/Makevars

@@ -1 +1,2 @@
-PKG_CPPFLAGS = -I../inst/include -DUSE_RCPP 
+CXX_STD = CXX11
+PKG_CPPFLAGS = -I../inst/include -DUSE_RCPP

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},

src/core.cpp

@@ -47,14 +47,14 @@ 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 ),
     lastDate( -1.0), 
     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<string, string>( capabilityName, componentName ) );
+        H_LOG(glog, Logger::DEBUG) << capabilityName << " registered to component " << componentName << "\n";
+
     }
 }
 

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<std::string, std::string> 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; i<N_HALO_FORCINGS; ++i) {
+        core->registerCapability(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<std::string, unitval>::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<std::string, unitval>::const_iterator forcing = forcings.find(forcing_name);
         if ( forcing != forcings.end() ) {
             // from the forcing map
             returnval = forcing->second;

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)