Merge pull request #298 from JGCRI/land-ocean-branch

Report land and ocean temperatures

Version 2.2.1

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: hector
 Title: The Hector Simple Climate Model
-Version: 2.2.0
+Version: 2.2.1
 Authors@R: c(person("Robert", "Link", email = "[email protected]", role = c("aut", "cre")),
 	     person("Corinne", "Hartin", email = "[email protected]", role = "aut"),
 	     person("Ben", "Bond-Lamberty", email = "[email protected]", role = "aut"),
@@ -21,6 +21,8 @@ Suggests:
     knitr,
     rmarkdown,
     ggplot2,
+    tidyr (>= 0.7),
+    dplyr (>= 0.7),
     hectortools
 Remotes:
     jgcri/hectortools

NAMESPACE

@@ -64,6 +64,7 @@ export(GLOBAL_TEMP)
 export(GLOBAL_TEMPEQ)
 export(HEAT_FLUX)
 export(LAND_CFLUX)
+export(LAND_TEMP)
 export(LIFETIME_SOIL)
 export(LIFETIME_STRAT)
 export(LL_DEBUG)
@@ -73,12 +74,14 @@ export(LL_WARNING)
 export(LUC_EMISSIONS)
 export(NATURAL_CH4)
 export(NATURAL_SO2)
+export(OCEAN_AIR_TEMP)
 export(OCEAN_C)
 export(OCEAN_CFLUX)
 export(OCEAN_C_DO)
 export(OCEAN_C_HL)
 export(OCEAN_C_IO)
 export(OCEAN_C_LL)
+export(OCEAN_SURFACE_TEMP)
 export(PCO2_HL)
 export(PCO2_LL)
 export(PH_HL)

R/RcppExports.R

@@ -790,6 +790,24 @@ GLOBAL_TEMPEQ <- function() {
     .Call('_hector_GLOBAL_TEMPEQ', PACKAGE = 'hector')
 }
 
+#' @describeIn temperature Average ocean surface temperature anomaly
+#' @export
+OCEAN_SURFACE_TEMP <- function() {
+    .Call('_hector_OCEAN_SURFACE_TEMP', PACKAGE = 'hector')
+}
+
+#' @describeIn temperature Average ocean air temperature anomaly
+#' @export
+OCEAN_AIR_TEMP <- function() {
+    .Call('_hector_OCEAN_AIR_TEMP', PACKAGE = 'hector')
+}
+
+#' @describeIn temperature Average land temperature anomaly
+#' @export
+LAND_TEMP <- function() {
+    .Call('_hector_LAND_TEMP', PACKAGE = 'hector')
+}
+
 #' @describeIn parameters Ocean heat diffusivity (\code{"cm2/s"})
 #' @export
 DIFFUSIVITY <- function() {

changelog.txt

@@ -1,3 +1,9 @@
+Hector 2.2.1
+============
+* Report global averages for land temperature, ocean air
+  temperature, and ocean surface temperature.
+* Fix bug that prevented CH4 emissions from being read.
+
 Hector 2.2.0
 ============
 * Add a new parameter: VOLCANIC_SCALE.  This parameter adjusts the

inst/include/component_data.hpp

@@ -242,6 +242,9 @@
 // temperature component
 #define D_ECS                   "S"
 #define D_GLOBAL_TEMP           "Tgav"
+#define D_LAND_TEMP             "Tgav_land"
+#define D_OCEAN_SURFACE_TEMP    "Tgav_ocean_ST"
+#define D_OCEAN_AIR_TEMP        "Tgav_ocean_air"
 #define D_GLOBAL_TEMPEQ         "Tgaveq"
 #define D_TGAV_CONSTRAIN        "tgav_constrain"
 #define D_SO2D_B                "so2d_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.0"
+#define MODEL_VERSION "2.2.1"
 
 #define OUTPUT_DIRECTORY "output/"
 

inst/include/temperature_component.hpp

@@ -137,7 +137,10 @@ class TemperatureComponent : public IModelComponent {
     unitval volscl;        //!< volcanic forcing scaling factor, unitless
 
     // Model outputs
-    unitval tgav;          //!< global temperature delta, deg C
+    unitval tgav;          //!< global average surface air temperature anomaly, deg C
+    unitval tgav_land;     //!< global average land surface air temperature anomaly, deg C
+    unitval tgav_oceanair; //!< global average ocean surface air temperature anomaly, deg C
+    unitval tgav_sst;      //!< global average ocean surface (water) temperature anomaly, deg C
     unitval tgaveq;        //!< equilibrium temp without ocean heat flux, currently set = tgav
     unitval flux_mixed;    //!< heat flux into mixed layer of ocean, W/m2
     unitval flux_interior; //!< heat flux into interior layer of ocean, W/m2

src/RcppExports.cpp

@@ -1315,6 +1315,36 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// OCEAN_SURFACE_TEMP
+String OCEAN_SURFACE_TEMP();
+RcppExport SEXP _hector_OCEAN_SURFACE_TEMP() {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    rcpp_result_gen = Rcpp::wrap(OCEAN_SURFACE_TEMP());
+    return rcpp_result_gen;
+END_RCPP
+}
+// OCEAN_AIR_TEMP
+String OCEAN_AIR_TEMP();
+RcppExport SEXP _hector_OCEAN_AIR_TEMP() {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    rcpp_result_gen = Rcpp::wrap(OCEAN_AIR_TEMP());
+    return rcpp_result_gen;
+END_RCPP
+}
+// LAND_TEMP
+String LAND_TEMP();
+RcppExport SEXP _hector_LAND_TEMP() {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    rcpp_result_gen = Rcpp::wrap(LAND_TEMP());
+    return rcpp_result_gen;
+END_RCPP
+}
 // DIFFUSIVITY
 String DIFFUSIVITY();
 RcppExport SEXP _hector_DIFFUSIVITY() {
@@ -1575,6 +1605,9 @@ static const R_CallMethodDef CallEntries[] = {
     {"_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_OCEAN_SURFACE_TEMP", (DL_FUNC) &_hector_OCEAN_SURFACE_TEMP, 0},
+    {"_hector_OCEAN_AIR_TEMP", (DL_FUNC) &_hector_OCEAN_AIR_TEMP, 0},
+    {"_hector_LAND_TEMP", (DL_FUNC) &_hector_LAND_TEMP, 0},
     {"_hector_DIFFUSIVITY", (DL_FUNC) &_hector_DIFFUSIVITY, 0},
     {"_hector_FLUX_MIXED", (DL_FUNC) &_hector_FLUX_MIXED, 0},
     {"_hector_FLUX_INTERIOR", (DL_FUNC) &_hector_FLUX_INTERIOR, 0},

src/ch4_component.cpp

@@ -172,6 +172,9 @@ unitval CH4Component::getData( const std::string& varName,
     } else if( varName == D_PREINDUSTRIAL_CH4 ) {
         H_ASSERT( date == Core::undefinedIndex(), "Date not allowed for preindustrial CH4" );
         returnval = M0;
+    } else if( varName == D_EMISSIONS_CH4 ) {
+        H_ASSERT(  date != Core::undefinedIndex(), "Date not allowed for CH4 emissions" );
+        returnval = CH4_emissions.get( date );
     } else {
         H_THROW( "Caller is requesting unknown variable: " + varName );
     }

src/forcing_component.cpp

@@ -340,7 +340,10 @@ void ForcingComponent::run( const double runToDate ) throw ( h_exception ) {
         H_LOG( logger, Logger::DEBUG ) << "forcing total is " << forcings[ D_RF_TOTAL ] << std::endl;
 
         //---------- Change to relative forcing ----------
-        // At this point, we've computed all absolute forcings. If base year, save those values
+        // Note that the code below assumes model is always consistently run from base-year forward. 
+        // Results will not be consistent if parameters are changed but base-year is not re-run.
+
+       // At this point, we've computed all absolute forcings. If base year, save those values
         if( runToDate==baseyear ) {
             H_LOG( logger, Logger::DEBUG ) << "** At base year! Storing current forcing values" << std::endl;
             baseyear_forcings = forcings;

src/rcpp_constants.cpp

@@ -998,6 +998,27 @@ String GLOBAL_TEMPEQ() {
 return D_GLOBAL_TEMPEQ;
 }
 
+//' @describeIn temperature Average ocean surface temperature anomaly
+//' @export
+// [[Rcpp::export]]
+String OCEAN_SURFACE_TEMP() {
+return D_OCEAN_SURFACE_TEMP;
+}
+
+//' @describeIn temperature Average ocean air temperature anomaly
+//' @export
+// [[Rcpp::export]]
+String OCEAN_AIR_TEMP() {
+return D_OCEAN_AIR_TEMP;
+}
+
+//' @describeIn temperature Average land temperature anomaly
+//' @export
+// [[Rcpp::export]]
+String LAND_TEMP() {
+return D_LAND_TEMP;
+}
+
 
 //' @describeIn parameters Ocean heat diffusivity (\code{"cm2/s"})
 //' @export

src/temperature_component.cpp

@@ -96,6 +96,9 @@ void TemperatureComponent::init( Core* coreptr ) {
 
     // Register the data we can provide
     core->registerCapability( D_GLOBAL_TEMP, getComponentName() );
+    core->registerCapability( D_LAND_TEMP, getComponentName() );
+    core->registerCapability( D_OCEAN_AIR_TEMP, getComponentName() );
+    core->registerCapability( D_OCEAN_SURFACE_TEMP, getComponentName() );
     core->registerCapability( D_GLOBAL_TEMPEQ, getComponentName() );
     core->registerCapability( D_FLUX_MIXED, getComponentName() );
     core->registerCapability( D_FLUX_INTERIOR, getComponentName() );
@@ -437,6 +440,12 @@ unitval TemperatureComponent::getData( const std::string& varName,
             returnval = tgav;
         } else if( varName == D_GLOBAL_TEMPEQ ) {
             returnval = tgaveq;
+        } else if( varName == D_LAND_TEMP ) {
+            returnval = tgav_land;
+        } else if( varName == D_OCEAN_SURFACE_TEMP ) {
+            returnval = tgav_sst;
+        } else if( varName == D_OCEAN_AIR_TEMP ) {
+            returnval = tgav_oceanair;
         } else if( varName == D_DIFFUSIVITY ) {
             returnval = diff;
         } else if( varName == D_AERO_SCALE ) {
@@ -465,6 +474,12 @@ unitval TemperatureComponent::getData( const std::string& varName,
 
         if( varName == D_GLOBAL_TEMP ) {
             returnval = unitval(temp[tstep], U_DEGC);
+        } else if( varName == D_LAND_TEMP ) {
+            returnval = unitval(temp_landair[tstep], U_DEGC);
+        } else if( varName == D_OCEAN_SURFACE_TEMP ) {
+            returnval = unitval(temp_sst[tstep], U_DEGC);
+        } else if( varName == D_OCEAN_AIR_TEMP ) {
+            returnval = bsi * unitval(temp_sst[tstep], U_DEGC);
         } else if( varName == D_GLOBAL_TEMPEQ ) {
             returnval = unitval(temp[tstep], U_DEGC);
         } else if( varName == D_FLUX_MIXED ) {
@@ -514,11 +529,17 @@ void TemperatureComponent::accept( AVisitor* visitor ) {
 
 void TemperatureComponent::setoutputs(int tstep)
 {
+    double temp_oceanair;
+
     flux_mixed.set( heatflux_mixed[tstep], U_W_M2, 0.0 );
     flux_interior.set( heatflux_interior[tstep], U_W_M2, 0.0 );
     heatflux.set( heatflux_mixed[tstep] + fso * heatflux_interior[tstep], U_W_M2, 0.0 );
     tgav.set(temp[tstep], U_DEGC, 0.0);
     tgaveq.set(temp[tstep], U_DEGC, 0.0); // per comment line 140 of temperature_component.hpp
+    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);    
 }
 
 }

tests/testthat/test_scenarios.R

@@ -3,6 +3,8 @@ context('Verify correctness of hector scenarios')
 inputdir <- system.file('input', package='hector')
 sampledir <- system.file('output', package='hector')
 testvars <- c(ATMOSPHERIC_CO2(), RF_TOTAL(), GLOBAL_TEMP())
+tempvars <- c(GLOBAL_TEMP(), OCEAN_AIR_TEMP(), LAND_TEMP())
+fland <- 0.29        # Global land area fraction
 dates <- 2000:2300
 
 
@@ -14,6 +16,7 @@ test_that('RCP scenarios are correct', {
         expect_true(inherits(hc, 'hcore'))
         run(hc)
         outdata <- fetchvars(hc, dates, testvars, scen)
+        tempdata <- fetchvars(hc, dates, tempvars, scen)
 
         ## Get the comparison data
         sampleoutfile <- sprintf('sample_outputstream_rcp%s.csv', rcp)
@@ -29,6 +32,11 @@ test_that('RCP scenarios are correct', {
         outdata$value <- signif(outdata$value, 4)
 
         expect_equivalent(outdata, sampledata, info=sprintf("Output doesn't match for scenario rcp%s", rcp))
+
+        td <- dplyr::select(tempdata, year, variable, value) %>% tidyr::spread(variable, value)
+        tgcomp <- fland*td$Tgav_land + (1.0-fland)*td$Tgav_ocean_air
+        expect_equal(tgcomp, td$Tgav, info=sprintf("Global temperature doesn't add up for scenario rcp%s", rcp))
+
         hc <- shutdown(hc)
     }
 })