From 31a31b7feb0710443755862e992022edcbf63166 Mon Sep 17 00:00:00 2001 From: kdorheim Date: Tue, 3 Feb 2026 18:31:43 -0500 Subject: [PATCH] example implementation of using a consistent atmospheric moles --- inst/include/ch4_component.hpp | 1 - inst/include/component_data.hpp | 2 +- inst/include/core.hpp | 5 ++++ inst/input/hector_ssp245.ini | 2 +- src/ch4_component.cpp | 45 +++++++++++++++++++++------------ src/core.cpp | 8 +++++- 6 files changed, 43 insertions(+), 20 deletions(-) diff --git a/inst/include/ch4_component.hpp b/inst/include/ch4_component.hpp index 589b04b2c..eceff0e84 100644 --- a/inst/include/ch4_component.hpp +++ b/inst/include/ch4_component.hpp @@ -61,7 +61,6 @@ class CH4Component : public IModelComponent { tseries CH4_constrain; // CH4 concentration constraint, ppbv CH4 tseries 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 diff --git a/inst/include/component_data.hpp b/inst/include/component_data.hpp index de2a9efc4..f117d5c97 100644 --- a/inst/include/component_data.hpp +++ b/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" diff --git a/inst/include/core.hpp b/inst/include/core.hpp index e50f82cee..612af2d17 100644 --- a/inst/include/core.hpp +++ b/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. diff --git a/inst/input/hector_ssp245.ini b/inst/input/hector_ssp245.ini index 2bd7259af..96e6c84e8 100644 --- a/inst/input/hector_ssp245.ini +++ b/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 diff --git a/src/ch4_component.cpp b/src/ch4_component.cpp index 8d89bcfa9..eb756078a 100644 --- a/src/ch4_component.cpp +++ b/src/ch4_component.cpp @@ -61,7 +61,6 @@ 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()); @@ -69,6 +68,9 @@ void CH4Component::init(Core *coreptr) { 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) diff --git a/src/core.cpp b/src/core.cpp index 6a2f64a56..7d4543153 100644 --- a/src/core.cpp +++ b/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); }