Enable global mean surface temperature as Hector output (#701)

Save global mean surface temperature and return it as Hector output. Update the Hector version in preparation for a new release.

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: hector
 Title: The Hector Simple Climate Model
-Version: 3.0.1
+Version: 3.1.0
 Authors@R: c(person("Kalyn", "Dorheim", 
 	        email = "[email protected]",
 	        role = c("aut", "cre"),

NAMESPACE

@@ -107,6 +107,7 @@ export(F_LITTERD)
 export(F_NPPD)
 export(F_NPPV)
 export(GLOBAL_TAS)
+export(GMST)
 export(HALON1211_CONSTRAIN)
 export(HALON1301_CONSTRAIN)
 export(HALON2402_CONSTRAIN)

NEWS.md

@@ -1,3 +1,7 @@
+# hector 3.1.0 
+
+Return global mean surface temperature as a Hector output.
+
 # hector 3.0.1
 
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.7617326.svg)](https://doi.org/10.5281/zenodo.7617326)

R/RcppExports.R

@@ -1498,6 +1498,12 @@ GLOBAL_TAS <- function() {
     .Call('_hector_GLOBAL_TAS', PACKAGE = 'hector')
 }
 
+#' @describeIn temperature Global mean surface temperature anomaly
+#' @export
+GMST <- function() {
+    .Call('_hector_GMST', PACKAGE = 'hector')
+}
+
 #' @describeIn temperature Average sea surface temperature anomaly
 #' @export
 SST <- function() {

inst/include/component_data.hpp

@@ -394,6 +394,7 @@
 #define D_LAND_TAS "land_tas"
 #define D_SST "sst"
 #define D_OCEAN_TAS "ocean_tas"
+#define D_GMST "gmst"
 #define D_LO_WARMING_RATIO "lo_warming_ratio"
 #define D_DIFFUSIVITY "diff"
 #define D_AERO_SCALE "alpha"

inst/include/temperature_component.hpp

@@ -133,6 +133,7 @@ class TemperatureComponent : public IModelComponent {
 
   // Time series arrays that are updated with each DOECLIM time-step
   std::vector<double> temp;
+  std::vector<double> temp_surface;
   std::vector<double> temp_landair;
   std::vector<double> temp_sst;
   std::vector<double> heatflux_mixed;
@@ -148,7 +149,8 @@ class TemperatureComponent : public IModelComponent {
   unitval volscl; //!< volcanic forcing scaling factor, unitless
 
   // Model outputs
-  unitval tas;       //!< global average surface air temperature anomaly, deg C
+  unitval tas;       //!< global average air temperature anomaly, deg C
+  unitval gmst_val; //!< global average surface temperature anomaly, deg C
   unitval tas_land;  //!< average air temperature anomaly over land, deg C
   unitval tas_ocean; //!< global average ocean surface air temperature anomaly,
                      //!< deg C

src/RcppExports.cpp

@@ -2507,6 +2507,16 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// GMST
+String GMST();
+RcppExport SEXP _hector_GMST() {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    rcpp_result_gen = Rcpp::wrap(GMST());
+    return rcpp_result_gen;
+END_RCPP
+}
 // SST
 String SST();
 RcppExport SEXP _hector_SST() {
@@ -2993,6 +3003,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_hector_AERO_SCALE", (DL_FUNC) &_hector_AERO_SCALE, 0},
     {"_hector_VOLCANIC_SCALE", (DL_FUNC) &_hector_VOLCANIC_SCALE, 0},
     {"_hector_GLOBAL_TAS", (DL_FUNC) &_hector_GLOBAL_TAS, 0},
+    {"_hector_GMST", (DL_FUNC) &_hector_GMST, 0},
     {"_hector_SST", (DL_FUNC) &_hector_SST, 0},
     {"_hector_OCEAN_TAS", (DL_FUNC) &_hector_OCEAN_TAS, 0},
     {"_hector_LAND_TAS", (DL_FUNC) &_hector_LAND_TAS, 0},

src/csv_outputstream_visitor.cpp

@@ -231,6 +231,7 @@ void CSVOutputStreamVisitor::visit(TemperatureComponent *c) {
   if (!core->outputEnabled(c->getComponentName()))
     return;
   STREAM_MESSAGE(csvFile, c, D_GLOBAL_TAS);
+  STREAM_MESSAGE(csvFile, c, D_GMST);
   STREAM_MESSAGE(csvFile, c, D_FLUX_MIXED);
   STREAM_MESSAGE(csvFile, c, D_FLUX_INTERIOR)
   STREAM_MESSAGE(csvFile, c, D_HEAT_FLUX);

src/rcpp_constants.cpp

@@ -1444,6 +1444,11 @@ String VOLCANIC_SCALE() { return D_VOLCANIC_SCALE; }
 // [[Rcpp::export]]
 String GLOBAL_TAS() { return D_GLOBAL_TAS; }
 
+//' @describeIn temperature Global mean surface temperature anomaly
+//' @export
+// [[Rcpp::export]]
+String GMST() { return D_GMST; }
+
 //' @describeIn temperature Average sea surface temperature anomaly
 //' @export
 // [[Rcpp::export]]

src/temperature_component.cpp

@@ -114,6 +114,7 @@ void TemperatureComponent::init(Core *coreptr) {
 
   // Register the data we can provide
   core->registerCapability(D_GLOBAL_TAS, getComponentName());
+  core->registerCapability(D_GMST, getComponentName());
   core->registerCapability(D_LAND_TAS, getComponentName());
   core->registerCapability(D_OCEAN_TAS, getComponentName());
   core->registerCapability(D_SST, getComponentName());
@@ -241,6 +242,7 @@ void TemperatureComponent::prepareToRun() {
   Ker.resize(ns);
 
   temp.resize(ns);
+  temp_surface.resize(ns);
   temp_landair.resize(ns);
   temp_sst.resize(ns);
   heatflux_mixed.resize(ns);
@@ -488,6 +490,7 @@ void TemperatureComponent::run(const double runToDate) {
 
   // Reset the endogenous varibales for this time step
   temp[tstep] = 0.0;
+  temp_surface[tstep] = 0.0;
   temp_landair[tstep] = 0.0;
   temp_sst[tstep] = 0.0;
   heat_mixed[tstep] = 0.0;
@@ -540,6 +543,8 @@ void TemperatureComponent::run(const double runToDate) {
   }
   temp[tstep] =
       flnd * temp_landair[tstep] + (1.0 - flnd) * bsi * temp_sst[tstep];
+  temp_surface[tstep] =
+      flnd * temp_landair[tstep] + (1.0 - flnd) * temp_sst[tstep];
 
   // If the user has supplied temperature data, use that instead
   if (tas_constrain.size() && runToDate >= tas_constrain.firstdate() &&
@@ -554,6 +559,9 @@ void TemperatureComponent::run(const double runToDate) {
         (temp[tstep] - (1.0 - flnd) * bsi * temp_sst[tstep]) / flnd;
     temp_sst[tstep] =
         (temp[tstep] - flnd * temp_landair[tstep]) / ((1.0 - flnd) * bsi);
+    temp_surface[tstep] =
+        flnd * temp_landair[tstep] + (1.0 - flnd) * temp_sst[tstep];
+
   }
 
   // Calculate ocean heat uptake [W/m^2]
@@ -608,6 +616,12 @@ unitval TemperatureComponent::getData(const std::string &varName,
     } else {
       returnval = unitval(temp[tstep], U_DEGC);
     }
+  } else if (varName == D_GMST) {
+      if (date == Core::undefinedIndex()) {
+          returnval = tas;
+      } else {
+          returnval = unitval(temp_surface[tstep], U_DEGC);
+      }
   } else if (varName == D_LAND_TAS) {
     if (lo_warming_ratio != 0) {
       if (date == Core::undefinedIndex()) {
@@ -745,6 +759,7 @@ void TemperatureComponent::setoutputs(int tstep) {
   heatflux.set(heatflux_mixed[tstep] + fso * heatflux_interior[tstep], U_W_M2,
                0.0);
   tas.set(temp[tstep], U_DEGC, 0.0);
+  gmst_val.set(temp_surface[tstep], U_DEGC, 0.0);
   tas_land.set(temp_landair[tstep], U_DEGC, 0.0);
   sst.set(temp_sst[tstep], U_DEGC, 0.0);
   temp_oceanair = bsi * temp_sst[tstep];

tests/testthat/test_atmosphere.R

@@ -47,7 +47,7 @@ test_that("Checking RF values", {
 test_that("Check Temp", {
 
     # Define the comparison dates
-    t_dates <- 1850:2100
+    t_dates <- 2020:2100
 
     # Set and run Hector
     hc <- newcore(ssp245)
@@ -62,4 +62,9 @@ test_that("Check Temp", {
     weighted_sum <- flnd * land + ocean * (1 - flnd)
     global_temp_values <- fetchvars(hc, t_dates, vars = GLOBAL_TAS())[["value"]]
     expect_equal(global_temp_values, weighted_sum, tolerance = 1e-5)
+
+    #  future global air temperature should be larger than surface temperature
+    gmst <- fetchvars(hc, t_dates, vars = GMST())[["value"]]
+    expect_true(all(global_temp_values >  gmst))
+
 })

vignettes/articles/ex_GMAT-GMST.Rmd

@@ -22,14 +22,11 @@ Although changes in GMST are relatively small in terms of numbers, it takes a
 great deal of heat to warm the entire planet by even a fraction of a degree.
 Measurements of the temperature changes provide evidence for anthropogenic climate change.
 
-However, Hector does not produce GMST as an output at this time. Hector does
-calculate the global mean air temperature (GMAT), which we can use to back-calculate
-GMST. GMAT is the air temperature measured 2 meters above the ground and is what is
+Both GMST and global mean air temperature (GMAT) are reported as Hector outputs.
+GMAT is the air temperature measured 2 meters above the ground and is what is
 presented in a daily weather report.^[Taylor, Jessica. "Analyzing Surface Air
 Temperatures by Latitude." NASA, 2018. https://mynasadata.larc.nasa.gov/mini-lesson/analyzing-surface-air-temperatures-latitude-student-activity]
-
-This vignette demonstrates how to set up the calculation between GMAT and GMST and
-how to plot results to observe differences over time.
+This vignette demonstrates how to set up the calculate GMST and compares the two temperature results. 
 
 (Note that Hector assumes that land air temperature and land surface temperature
 are equivalent.^[Kriegler, Elmar. "Imprecise Probability Analysis