Merge pull request #307 from JGCRI/spinup-conc

Constrain CO2 conc to be equal to preindustrial during spinup.
JGCRI · May 20, 2019 · 16c480a · 16c480a
2 parents fddfed5 + 88fe909
commit 16c480a
Show file tree

Hide file tree

Showing 11 changed files with 108,809 additions and 108,773 deletions.
diff --git a/.Rbuildignore b/.Rbuildignore
@@ -27,7 +27,8 @@
 ^src/testing$
 ^misc$
 ^data-raw$
+^analysis$
 
 # pkgdown config and output
 ^inst/_pkgdown.yml
-^docs/
+^docs/
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: hector
 Title: The Hector Simple Climate Model
-Version: 2.2.2
+Version: 2.3.0
 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"),

diff --git a/changelog.txt b/changelog.txt
@@ -1,3 +1,11 @@
+Hector 2.3.0
+============
+* Constrain atmospheric CO2 concentration to be equal to the
+  preindustrial concentration during the model spinup.  The result of
+  this is that the concentration at the beginning of the simulation
+  will be equal to the value specified in the `PREINDUSTRIAL_CO2`
+  parameter, which was not the case previously.
+
 Hector 2.2.2
 ============
 * Fix bug that was causing requests for H2O forcing in the R interface

diff --git 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.2"
+#define MODEL_VERSION "2.3.0"
 
 #define OUTPUT_DIRECTORY "output/"
 

diff --git a/inst/output/sample_outputstream_rcp26.csv b/inst/output/sample_outputstream_rcp26.csv
diff --git a/inst/output/sample_outputstream_rcp45.csv b/inst/output/sample_outputstream_rcp45.csv
diff --git a/inst/output/sample_outputstream_rcp60.csv b/inst/output/sample_outputstream_rcp60.csv
diff --git a/inst/output/sample_outputstream_rcp85.csv b/inst/output/sample_outputstream_rcp85.csv
diff --git a/man/haloforcings.Rd b/man/haloforcings.Rd
diff --git a/src/simpleNbox.cpp b/src/simpleNbox.cpp
@@ -655,13 +655,24 @@ void SimpleNbox::stashCValues( double t, const double c[] )
     masstot = sum;
 
     // If user has supplied Ca values, adjust atmospheric C to match
-    if( !core->inSpinup() && Ca_constrain.size() && t <= Ca_constrain.lastdate() ) {
-
-        H_LOG( logger, Logger::WARNING ) << "** Constraining atmospheric CO2 to user-supplied value" << std::endl;
+    if(core->inSpinup() ||
+       ( Ca_constrain.size() && t <= Ca_constrain.lastdate() )) {
+
+        unitval atmos_cpool_to_match;
+        unitval atmppmv;
+        if(core->inSpinup()) {
+            atmos_cpool_to_match.set(C0.value(U_PPMV_CO2) / PGC_TO_PPMVCO2, U_PGC);
+            atmppmv.set(C0.value(U_PPMV_CO2), U_PPMV_CO2);
+        }
+        else {
+            H_LOG( logger, Logger::WARNING ) << "** Constraining atmospheric CO2 to user-supplied value" << std::endl; 
+            atmos_cpool_to_match.set(Ca_constrain.get(t).value( U_PPMV_CO2 ) /
+                                     PGC_TO_PPMVCO2, U_PGC);
+            atmppmv.set(Ca_constrain.get(t).value(U_PPMV_CO2), U_PPMV_CO2);
+        }
 
-        unitval atmos_cpool_to_match( Ca_constrain.get( t ).value( U_PPMV_CO2 ) / PGC_TO_PPMVCO2, U_PGC );
         residual = atmos_c - atmos_cpool_to_match;
-        H_LOG( logger,Logger::DEBUG ) << t << "- have " << Ca << " want " << Ca_constrain.get( t ).value( U_PPMV_CO2 ) << std::endl;
+        H_LOG( logger,Logger::DEBUG ) << t << "- have " << Ca << " want " <<  atmppmv.value( U_PPMV_CO2 ) << std::endl;
         H_LOG( logger,Logger::DEBUG ) << t << "- have " << atmos_c << " want " << atmos_cpool_to_match << "; residual = " << residual << std::endl;
 
         // Transfer C from atmosphere to deep ocean and update our C and Ca variables

diff --git a/tests/testthat/test_parameters.R b/tests/testthat/test_parameters.R
@@ -22,6 +22,22 @@ test_that('Rerunning spinup produces minimal change', {
     shutdown(hc)
 })
 
+test_that('Initial CO2 concentration equals preindustrial', {
+    hc <- newcore(rcp45, suppresslogging=TRUE)
+    run(hc, 1800)
+    initcval <- fetchvars(hc, 1745, ATMOSPHERIC_CO2())
+    preind <- fetchvars(hc, NA, PREINDUSTRIAL_CO2())
+    expect_equal(initcval$value, preind$value)
+
+    setvar(hc, NA, PREINDUSTRIAL_CO2(), 285, "ppmv CO2")
+    reset(hc)
+    run(hc, 1800)
+    initcval <- fetchvars(hc, 1745, ATMOSPHERIC_CO2())
+    preind <- fetchvars(hc, NA, PREINDUSTRIAL_CO2())
+    expect_equal(initcval$value, preind$value)
+})
+
+
 test_that('Lowering initial CO2 lowers output CO2', {
     hc <- newcore(rcp45, suppresslogging = TRUE)
     run(hc, 2100)
@@ -142,7 +158,7 @@ test_that('Increasing volcanic forcing scaling factor increases the effect of vo
 })
 
 test_that('Decreasing vegetation NPP fraction increases CO2 concentration', {
-    # More NPP to vegetation means less C to soil, where it decomposes. 
+    # More NPP to vegetation means less C to soil, where it decomposes.
     hc <- newcore(rcp45, suppresslogging = TRUE)
     run(hc, 2100)
     qdates <- 2000:2100                 # limit to years where temperature