Consistent Unit Conversion & Parameters

inst/include/ch4_component.hpp

@@ -61,7 +61,6 @@ class CH4Component : public IModelComponent {
   tseries<unitval> CH4_constrain; // CH4 concentration constraint, ppbv CH4
   tseries<unitval> CH4N;          // annual natural emissions, Tg CH4/yr
   unitval M0;                     // preindustrial methane, ppbv CH4
-  unitval UC_CH4; // conversion factor between emissions and concentration
   unitval Tsoil;  // lifetime of soil sink, yr
   unitval Tstrat; // lifetime of tropospheric sink, yr
 

inst/include/component_data.hpp

@@ -29,6 +29,7 @@
 #define D_MAX_SPINUP "max_spinup"
 #define D_ENABLED "enabled"
 #define D_OUTPUT_ENABLED "output"
+#define D_ATM_MOLS "atmosphere_mol"
 
 // bc component
 #define D_EMISSIONS_BC "BC_emissions"
@@ -262,7 +263,6 @@
 #define D_EMISSIONS_CH4 "CH4_emissions"
 #define D_CONSTRAINT_CH4 "CH4_constrain"
 #define D_NATURAL_CH4 "CH4N"
-#define D_CONVERSION_CH4 "UC_CH4"
 #define D_LIFETIME_SOIL "Tsoil"
 #define D_LIFETIME_STRAT "Tstrat"
 

inst/include/core.hpp

@@ -83,6 +83,7 @@ class Core : public IVisitable {
   double getStartDate() const { return startDate; };
   double getEndDate() const { return endDate; };
   double getTrackingDate() const { return trackingDate; };
+  double getAtmMols() const { return atm_mols; };
   std::string getTrackingData() const;
   double getCurrentDate() const { return lastDate; }
   std::string getRun_name() const { return run_name; };
@@ -163,6 +164,10 @@ class Core : public IVisitable {
   //! Maximum number of spinup steps allowed.
   int max_spinup;
 
+  //------------------------------------------------------------------------------
+  //! Total moles in the atmopshere based on the atmosperic fry air constant
+  double atm_mols;
+
   //------------------------------------------------------------------------------
   //! A comparison object to ensure modelComponents are ordered according to
   //! dependencies.

inst/input/hector_ssp245.ini

@@ -7,6 +7,7 @@ endDate=2300
 trackingDate=9999	; year to start tracking (only carbon currently)
 do_spinup=1			; if 1, spin up model before running (default=1)
 max_spinup=2000		; maximum steps allowed for spinup (default=2000)
+atmosphere_mol=1.727e20	; number of moles in the atmosphere, mols 
 
 ;------------------------------------------------------------------------
 [ocean]
@@ -120,7 +121,6 @@ SV=csv:tables/core_inputs.csv   		; volcanic radiative forcing time series
 M0= 731.41  		; 731.41 ppb preindustrial methane IPCC AR6 Table 7.SM.1, the CH4 forcing equations is calibrated to a M0 of 731.41 ppb
 Tsoil=120 			; lifetime of soil sink (years) Myhre et al., 2013, doi: 10.1017/cbo9781107415324.018
 Tstrat=150          ; lifetime of tropospheric sink (years) Myhre et al., 2013, doi: 10.1017/cbo9781107415324.018
-UC_CH4=2.78			; Tg(CH4)/ppb unit conversion between emissions and concentrations
 CH4N=csv:tables/core_inputs.csv			; Natural CH4 emissions (Tgrams) see PR #786
 CH4_emissions=csv:tables/ssp245_emiss-constraints_rf.csv     ; emissions time series
 ;CH4_constrain=csv:tables/ssp245_emiss-constraints_rf.csv    ; CH4 concentration constraint

src/ch4_component.cpp

@@ -61,14 +61,16 @@ void CH4Component::init(Core *coreptr) {
   core->registerCapability(D_PREINDUSTRIAL_CH4, getComponentName());
   core->registerCapability(D_LIFETIME_STRAT, getComponentName());
   core->registerCapability(D_LIFETIME_SOIL, getComponentName());
-  core->registerDependency(D_LIFETIME_OH, getComponentName());
   // ...and what input data that we can accept
   core->registerInput(D_EMISSIONS_CH4, getComponentName());
   core->registerInput(D_NATURAL_CH4, getComponentName());
   core->registerInput(D_CONSTRAINT_CH4, getComponentName());
   core->registerInput(D_PREINDUSTRIAL_CH4, getComponentName());
   core->registerInput(D_LIFETIME_STRAT, getComponentName());
   core->registerInput(D_LIFETIME_SOIL, getComponentName());
+  // ... and what is needed from else where
+  core->registerDependency(D_LIFETIME_OH, getComponentName());
+
 }
 
 //------------------------------------------------------------------------------
@@ -113,9 +115,6 @@ void CH4Component::setData(const string &varName, const message_data &data) {
     } else if (varName == D_LIFETIME_STRAT) {
       H_ASSERT(data.date == Core::undefinedIndex(), "date not allowed");
       Tstrat = data.getUnitval(U_YRS);
-    } else if (varName == D_CONVERSION_CH4) {
-      H_ASSERT(data.date == Core::undefinedIndex(), "date not allowed");
-      UC_CH4 = data.getUnitval(U_TG_PPBV);
     } else if (varName == D_NATURAL_CH4) {
       H_ASSERT(data.date != Core::undefinedIndex(), "date required");
       CH4N.set(data.date, data.getUnitval(U_TG_CH4));
@@ -157,24 +156,38 @@ void CH4Component::run(const double runToDate) {
     CH4.set(runToDate, CH4_constrain.get(runToDate));
   } else {
 
-    // modified from Wigley et al, 2002
-    // https://doi.org/10.1175/1520-0442(2002)015%3C2690:RFDTRG%3E2.0.CO;2
-    const double current_ch4em = CH4_emissions.get(runToDate).value(U_TG_CH4);
+   // When emission driven hector uses a modified implementation of
+   // Wigley et al, 2002 (https://doi.org/10.1175/1520-0442(2002)015%3C2690:RFDTRG%3E2.0.CO;2)
+
+    // Unit conversion realted terms
+    #define PG_C_TO_TG_CH4 (1000.0 * 16.04 / 12.01) // should this be #define or would it make sense to do something like constexpr?
+
+    // Determine the conversion factor used to convert from Tg CH4 --> ppbv
+    double atm_mols = core->getAtmMols();
+    const double Tg_to_mol = 1 * 1e12 * (1/16.04);         // 1 Tg CH4 --> moles of CH4
+    const double mol_ratio = Tg_to_mol / atm_mols;         // ratio of CH4 moles to total moles in atmosphere
+    double UC_CH4 = 1/(mol_ratio * 1e9);               // convert from decimal ratio to parts per billion
+    UC_CH4 = std::round(UC_CH4 * 1000.0)/1000.0;       // control the number of digits
+
+    // Get the OH methane sink
     const double current_toh =
-        core->sendMessage(M_GETDATA, D_LIFETIME_OH, runToDate).value(U_YRS);
-    H_LOG(logger, Logger::DEBUG)
-        << "Year " << runToDate << " current_toh = " << current_toh
-        << std::endl;
-
-// Permafrost thaw produces CH4 emissions
-#define PG_C_TO_TG_CH4 (1000.0 * 16.04 / 12.01)
+      core->sendMessage(M_GETDATA, D_LIFETIME_OH, runToDate).value(U_YRS);
+      H_LOG(logger, Logger::DEBUG)
+          << "Year " << runToDate << " current_toh = " << current_toh
+          << std::endl;
+
+    // Get the different sources of CH4
+    // Anthropogenic emissions
+    const double current_ch4em = CH4_emissions.get(runToDate).value(U_TG_CH4);
+    // Permafrost thaw produces CH4 emissions
     const double rh_ch4 =
         core->sendMessage(M_GETDATA, D_RH_CH4).value(U_PGC_YR) * PG_C_TO_TG_CH4;
-
     // Additional, background CH4 natural emissions
     const double ch4n = CH4N.get(runToDate).value(U_TG_CH4);
+
+
     const double emisTocon =
-        (current_ch4em + rh_ch4 + ch4n) / UC_CH4.value(U_TG_PPBV);
+        (current_ch4em + rh_ch4 + ch4n) / UC_CH4;
     const double previous_ch4 = CH4.get(oldDate);
 
     H_LOG(logger, Logger::DEBUG)

src/core.cpp

@@ -93,7 +93,8 @@ void Core::init() {
   // We handle this input; need to register it here so that R users can change
   // it
   registerInput(D_TRACKING_DATE, CORE_COMPONENT_NAME);
-
+  registerInput(D_ATM_MOLS, CORE_COMPONENT_NAME);
+
   IModelComponent *temp;
 
   temp = new SimpleNbox();
@@ -232,6 +233,9 @@ void Core::setData(const string &componentName, const string &varName,
       } else if (varName == D_TRACKING_DATE) {
         H_ASSERT(data.date == undefinedIndex(), "date not allowed");
         trackingDate = data.getUnitval(U_UNITLESS);
+      } else if (varName == D_ATM_MOLS) {
+        H_ASSERT(data.date == undefinedIndex(), "date not allowed");
+        atm_mols = data.getUnitval(U_MOL);
       } else if (varName == D_DO_SPINUP) {
         H_ASSERT(data.date == undefinedIndex(), "date not allowed");
         do_spinup = (data.getUnitval(U_UNDEFINED) > 0);
@@ -699,6 +703,8 @@ unitval Core::getData(const std::string &varName, const double date) {
     H_ASSERT(date == Core::undefinedIndex(),
              "Date not allowed for tracking date");
     returnval = unitval(trackingDate, U_UNITLESS);
+  } else if (varName == D_ATM_MOLS) {
+    returnval = unitval(atm_mols, U_MOL);
   } else {
     H_THROW("Caller is requesting unknown variable: " + varName);
   }