Skip to content

Commit

Permalink
Set current LUC, FFI, DACCS once at beginning of timestep
Browse files Browse the repository at this point in the history 
Remove old `ffi()`, etc. functions; see #643
  • Loading branch information
@bpbond
bpbond committed Oct 31, 2022
1 parent 797d272 commit b0c7663
Show file tree
Hide file tree
Showing 3 changed files with 35 additions and 77 deletions.
23 changes: 12 additions & 11 deletions inst/include/simpleNbox.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,9 +172,9 @@ class SimpleNbox : public CarbonCycleModel {
*****************************************************************/

// Carbon fluxes
tseries<unitval>
tseries<fluxpool>
ffiEmissions; //!< fossil fuels and industry emissions, Pg C/yr
tseries<unitval>
tseries<fluxpool>
daccsUptake; //!< direct air carbon capture and storage, Pg C/yr
tseries<fluxpool> lucEmissions; //!< land use change emissions, Pg C/yr
tseries<fluxpool> lucUptake; //!< land use change uptake, Pg C/yr
Expand Down Expand Up @@ -205,14 +205,19 @@ class SimpleNbox : public CarbonCycleModel {
// Initial fluxes
fluxpool_stringmap npp_flux0; //!< preindustrial NPP

// Atmospheric CO2, temperature, and their effects
// Atmospheric CO2, temperature, and their effects
fluxpool C0; //!< preindustrial [CO2], ppmv

double_stringmap beta; //!< shape of CO2 response
double_stringmap
warmingfactor; //!< regional warming relative to global (1.0=same)
// Slowly-changing variables
// These get computed only once per year, in slowparameval()
fluxpool current_luc_e, //!< Current year LUC emissions (/yr)
current_luc_u, //!< Current year LUC uptake (/yr)
current_ffi_e, //!< Current year FFI emissions (/yr)
current_daccs_u; //!< Current year DACCS uptake (/yr)
double_stringmap beta; //!< shape of CO2 response
double_stringmap
q10_rh; //!< Q10 for heterotrophic respiration (1.0=no response, unitless)
warmingfactor; //!< regional warming relative to global (1.0=same)
double_stringmap q10_rh; //!< Q10 for heterotrophic respiration (1.0=no response, unitless)

/*****************************************************************
* Functions computing sub-elements of the carbon cycle
Expand All @@ -234,10 +239,6 @@ class SimpleNbox : public CarbonCycleModel {
double time = Core::undefinedIndex()) const; //!< calculates RH for a biome
fluxpool sum_rh(double time = Core::undefinedIndex())
const; //!< calculates RH, global total
fluxpool ffi(double t, bool in_spinup) const;
fluxpool ccs(double t, bool in_spinup) const;
fluxpool luc_emission(double t, bool in_spinip) const;
fluxpool luc_uptake(double t, bool in_spinip) const;

/*****************************************************************
* Private helper functions
Expand Down
83 changes: 19 additions & 64 deletions src/simpleNbox-runtime.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,8 +223,8 @@ void SimpleNbox::stashCValues(double t, const double c[]) {
// This changes only when we start to apportion NPP and RH to biomes

// get the UNTRACKED earth emissions (ffi) and uptake (ccs)
fluxpool ffi_untracked = ffi(t, in_spinup);
fluxpool ccs_untracked = ccs(t, in_spinup);
fluxpool ffi_untracked = current_ffi_e;
fluxpool ccs_untracked = current_daccs_u;
// now construct the TRACKED versions
// because earth_c is tracked, ffi_flux and ccs_flux automatically
// become tracked as well
Expand All @@ -239,9 +239,8 @@ void SimpleNbox::stashCValues(double t, const double c[]) {
fluxpool ao_flux = omodel->get_aoflux();

// Land-use change emissions and uptake
fluxpool luc_e_untracked = luc_emission(t, in_spinup);
fluxpool luc_u_untracked = luc_uptake(t, in_spinup);

fluxpool luc_e_untracked = current_luc_e;
fluxpool luc_u_untracked = current_luc_u;
fluxpool luc_fav_flux = atmos_c.flux_from_fluxpool(luc_u_untracked * f_lucv);
fluxpool luc_fad_flux = atmos_c.flux_from_fluxpool(luc_u_untracked * f_lucd);
fluxpool luc_fas_flux =
Expand Down Expand Up @@ -554,51 +553,7 @@ fluxpool SimpleNbox::sum_rh(double time) const {
return total;
}

//------------------------------------------------------------------------------
/*! \brief Compute fossil fuel industrial (FFI) emissions (when input
* emissions > 0) \returns FFI flux from earth to atmosphere
*/
fluxpool SimpleNbox::ffi(double t, bool in_spinup) const {
if (!in_spinup) { // no perturbation allowed if in spinup
return fluxpool(ffiEmissions.get(t).value(U_PGC_YR), U_PGC_YR);
}
return fluxpool(0.0, U_PGC_YR);
}

//------------------------------------------------------------------------------
/*! \brief Compute carbon capture storage (CCS) flux (when input emissions
* < 0) \returns CCS flux from atmosphere to earth
*/
fluxpool SimpleNbox::ccs(double t, bool in_spinup) const {
if (!in_spinup) { // no perturbation allowed if in spinup
return fluxpool(daccsUptake.get(t).value(U_PGC_YR), U_PGC_YR);
}
return fluxpool(0.0, U_PGC_YR);
}

//------------------------------------------------------------------------------
/*! \brief Compute land use change emissions (when input emissions > 0)
* \returns luc flux from land to atmosphere
*/
fluxpool SimpleNbox::luc_emission(double t, bool in_spinup) const {
if (!in_spinup) { // no perturbation allowed if in spinup
return lucEmissions.get(t);
}
return fluxpool(0.0, U_PGC_YR);
}

//------------------------------------------------------------------------------
/*! \brief Compute land use change uptake (when input emissions < 0)
* \returns luc flux from atmosphere to land
*/
fluxpool SimpleNbox::luc_uptake(double t, bool in_spinup) const {
if (!in_spinup) { // no perturbation allowed if in spinup
return lucUptake.get(t);
}
return fluxpool(0.0, U_PGC_YR);
}

//------------------------------------------------------------------------------
///------------------------------------------------------------------------------
/*! \brief Compute model fluxes for a time step
* \param[in] t time
* \param[in] c carbon pools (no units)
Expand Down Expand Up @@ -659,15 +614,6 @@ int SimpleNbox::calcderivs(double t, const double c[], double dcdt[]) const {
fluxpool(detritus_c.at(biome).value(U_PGC) * 0.6, U_PGC_YR);
}

// Annual fossil fuels and industry emissions and atmosphere CO2 capture (CCS
// or whatever)
fluxpool ffi_flux_current = ffi(t, in_spinup);
fluxpool ccs_flux_current = ccs(t, in_spinup);

// Annual land use change emissions
fluxpool luc_emission_current = luc_emission(t, in_spinup);
fluxpool luc_uptake_current = luc_uptake(t, in_spinup);

// Land-use change emissions come from veg, detritus, and soil
fluxpool luc_fva = luc_emission_current * f_lucv;
fluxpool luc_fda = luc_emission_current * f_lucd;
Expand All @@ -686,7 +632,7 @@ int SimpleNbox::calcderivs(double t, const double c[], double dcdt[]) const {
if (!in_spinup && NBP_constrain.size() && NBP_constrain.exists(rounded_t)) {
// Compute how different we are from the user-specified constraint
const double nbp = npp_current.value(U_PGC_YR) - rh_current.value(U_PGC_YR) -
luc_emission_current.value(U_PGC_YR) + luc_uptake_current.value(U_PGC_YR);
current_luc_e.value(U_PGC_YR) + current_luc_u.value(U_PGC_YR);
const unitval diff = NBP_constrain.get(rounded_t) - unitval(nbp, U_PGC_YR);

// Adjust total NPP and total RH equally (but not LUC, which is an input)
Expand All @@ -709,9 +655,9 @@ int SimpleNbox::calcderivs(double t, const double c[], double dcdt[]) const {

// Compute fluxes
dcdt[SNBOX_ATMOS] = // change in atmosphere pool
ffi_flux_current.value(U_PGC_YR) - ccs_flux_current.value(U_PGC_YR) +
luc_emission_current.value(U_PGC_YR) -
luc_uptake_current.value(U_PGC_YR) + ch4ox_current.value(U_PGC_YR) -
current_ffi_e.value(U_PGC_YR) - current_daccs_u.value(U_PGC_YR) +
current_luc_e.value(U_PGC_YR) - current_luc_u.value(U_PGC_YR) +
ch4ox_current.value(U_PGC_YR) -
ocean_uptake.value(U_PGC_YR) + ocean_release.value(U_PGC_YR) -
npp_current.value(U_PGC_YR) + rh_current.value(U_PGC_YR);
dcdt[SNBOX_VEG] = // change in vegetation pool
Expand All @@ -728,7 +674,7 @@ int SimpleNbox::calcderivs(double t, const double c[], double dcdt[]) const {
dcdt[SNBOX_OCEAN] = // change in ocean pool
ocean_uptake.value(U_PGC_YR) - ocean_release.value(U_PGC_YR);
dcdt[SNBOX_EARTH] = // change in earth pool
-ffi_flux_current.value(U_PGC_YR) + ccs_flux_current.value(U_PGC_YR);
-current_ffi_e.value(U_PGC_YR) + current_daccs_u.value(U_PGC_YR);

return omodel_err;
}
Expand All @@ -745,6 +691,15 @@ int SimpleNbox::calcderivs(double t, const double c[], double dcdt[]) const {
void SimpleNbox::slowparameval(double t, const double c[]) {
omodel->slowparameval(t, c); // pass msg on to ocean model

// Set this year's LUC and FFI/DACCS emissions and uptake
// We do this here, and not allow interpolation of their time series,
// so that pulse tests work correctly (see #643)
fluxpool zero_flux(0.0, U_PGC_YR);
current_luc_e = in_spinup ? zero_flux : lucEmissions.get(t);
current_luc_u = in_spinup ? zero_flux : lucUptake.get(t);
current_ffi_e = in_spinup ? zero_flux : ffiEmissions.get(t);
current_daccs_u = in_spinup ? zero_flux : daccsUptake.get(t);

// Compute CO2 fertilization factor globally (and for each biome specified)
for (auto biome : biome_list) {
if (in_spinup) {
Expand Down
6 changes: 4 additions & 2 deletions src/simpleNbox.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -287,12 +287,14 @@ void SimpleNbox::setData(const std::string &varName, const message_data &data) {
H_ASSERT(data.date != Core::undefinedIndex(), "date required");
H_ASSERT(biome == SNBOX_DEFAULT_BIOME,
"fossil fuels and industry emissions must be global");
ffiEmissions.set(data.date, data.getUnitval(U_PGC_YR));
unitval ffi = data.getUnitval(U_PGC_YR);
ffiEmissions.set(data.date, fluxpool(ffi.value(U_PGC_YR), U_PGC_YR));
} else if (varNameParsed == D_DACCS_UPTAKE) {
H_ASSERT(data.date != Core::undefinedIndex(), "date required");
H_ASSERT(biome == SNBOX_DEFAULT_BIOME,
"direct air carbon capture and storage must be global");
daccsUptake.set(data.date, data.getUnitval(U_PGC_YR));
unitval daccs = data.getUnitval(U_PGC_YR);
daccsUptake.set(data.date, fluxpool(daccs.value(U_PGC_YR), U_PGC_YR));
} else if (varNameParsed == D_LUC_EMISSIONS) {
H_ASSERT(data.date != Core::undefinedIndex(), "date required");
unitval luc = data.getUnitval(U_PGC_YR);
Expand Down

0 comments on commit b0c7663

Please sign in to comment.