SIP328 fix wood storage bookkeeping so handling of plantWoodCStorageDelta is consistent

src/sipnet/events.c

@@ -479,15 +479,18 @@ void processEvents(void) {
         const double fracRB = harvParams->fractionRemovedBelow;
         const double fracTB = harvParams->fractionTransferredBelow;
 
-        const double woodC = envi.plantWoodC + envi.plantWoodCStorageDelta;
+        const double aboveHarvestFrac = fracRA + fracTA;
+        const double storageWoodC = fmax(0.0, envi.plantWoodCStorageDelta);
+        const double woodC = getStorageBackedWoodCarbon();
         // Litter increase
         double litterAdd = fracTA * (envi.plantLeafC + woodC);
         double soilAdd = fracTB * (envi.fineRootC + envi.coarseRootC);
 
         // Pool reductions, counting both mass moved to litter and removed by
         // the harvest itself. Above-ground changes:
-        const double leafDelta = -envi.plantLeafC * (fracRA + fracTA);
-        const double woodDelta = -woodC * (fracRA + fracTA);
+        const double leafDelta = -envi.plantLeafC * aboveHarvestFrac;
+        const double woodDelta = -envi.plantWoodC * aboveHarvestFrac;
+        const double storageWoodDelta = -storageWoodC * aboveHarvestFrac;
         // Below-ground changes:
         const double fineDelta = -envi.fineRootC * (fracRB + fracTB);
         const double coarseDelta = -envi.coarseRootC * (fracRB + fracTB);
@@ -502,6 +505,7 @@ void processEvents(void) {
         fluxes.eventSoilC += soilAdd / climLen;
         fluxes.eventLeafC += leafDelta / climLen;
         fluxes.eventWoodC += woodDelta / climLen;
+        fluxes.eventWoodStorageDelta += storageWoodDelta / climLen;
         fluxes.eventFineRootC += fineDelta / climLen;
         fluxes.eventCoarseRootC += coarseDelta / climLen;
 
@@ -601,6 +605,7 @@ void updatePoolsForEvents(void) {
   // CARBON
   // Harvest and planting events
   envi.plantWoodC += fluxes.eventWoodC * climate->length;
+  envi.plantWoodCStorageDelta += fluxes.eventWoodStorageDelta * climate->length;
   envi.plantLeafC += fluxes.eventLeafC * climate->length;
 
   // Harvest and fertilization events

src/sipnet/sipnet.c

@@ -444,8 +444,7 @@ void outputHeader(FILE *out) {
 void outputState(FILE *out, int year, int day, double time) {
 
   fprintf(out, "%4d %3d %5.2f %10.2f %10.2f %12.2f ", year, day, time,
-          (envi.plantWoodC + envi.plantWoodCStorageDelta), envi.plantLeafC,
-          trackers.woodCreation);
+          getPlantWoodCTotal(), envi.plantLeafC, trackers.woodCreation);
   fprintf(out, "%8.2f ", envi.soilC);
   fprintf(out, "%11.2f %9.2f ", envi.coarseRootC, envi.fineRootC);
   fprintf(out, "%8.2f %10.3f %15.3f %8.2f ", envi.litterC, envi.soilWater,
@@ -1041,8 +1040,7 @@ void vegResp(double *folResp, double *woodResp, double baseFolResp) {
   // end snowpack addition
 
   // :: from [1], eq (A19)
-  *woodResp = params.baseVegResp *
-              (envi.plantWoodC + envi.plantWoodCStorageDelta) *
+  *woodResp = params.baseVegResp * getStorageBackedWoodCarbon() *
               pow(params.vegRespQ10, climate->tair / 10.0);
 }
 
@@ -1069,8 +1067,7 @@ void vegResp2(double *folResp, double *woodResp, double *growthResp,
     *folResp *= params.frozenSoilFolREff;  // allows foliar resp. to be shutdown
                                            // by a given fraction in winter
   }
-  *woodResp = params.baseVegResp *
-              (envi.plantWoodC + envi.plantWoodCStorageDelta) *
+  *woodResp = params.baseVegResp * getStorageBackedWoodCarbon() *
               pow(params.vegRespQ10, climate->tair / 10.0);
 
   // Rg is a fraction of the recent mean NPP
@@ -1329,8 +1326,7 @@ void calcRootAndWoodFluxes(void) {
 
   // Wood litter, in g C * m^-2 ground area * day^-1
   // turnover rate is fraction lost per day
-  fluxes.woodLitter =
-      (envi.plantWoodC + envi.plantWoodCStorageDelta) * params.woodTurnoverRate;
+  fluxes.woodLitter = getStorageBackedWoodCarbon() * params.woodTurnoverRate;
 
   // :: from [3], root model description
   calcRootResp(&fluxes.rCoarseRoot, params.coarseRootQ10,

src/sipnet/state.c

@@ -1,5 +1,7 @@
 // Definition of global state variables
 
+#include <math.h>
+
 #include "state.h"
 
 // linked list of climate variables
@@ -12,3 +14,11 @@ Envi envi;  // state variables
 Trackers trackers;
 PhenologyTrackers phenologyTrackers;
 Fluxes fluxes;
+
+double getPlantWoodCTotal(void) {
+  return envi.plantWoodC + envi.plantWoodCStorageDelta;
+}
+
+double getStorageBackedWoodCarbon(void) {
+  return fmax(envi.plantWoodC, getPlantWoodCTotal());
+}

src/sipnet/state.h

@@ -442,6 +442,9 @@ typedef struct Environment {
 // Global var
 extern Envi envi;  // state variables
 
+double getPlantWoodCTotal(void);
+double getStorageBackedWoodCarbon(void);
+
 // fluxes as per-day rates
 typedef struct FluxVars {
   // ****************************************
@@ -563,6 +566,8 @@ typedef struct FluxVars {
   double eventLeafC;
   // plantWoodC addition
   double eventWoodC;
+  // plantWoodCStorageDelta addition
+  double eventWoodStorageDelta;
   // plantFineRootC addition
   double eventFineRootC;
   // plantCoarseRootC addition

tests/sipnet/test_modeling/Makefile

@@ -8,7 +8,7 @@ LDFLAGS=-L$(ROOT_DIR)/libs
 LDLIBS=-lsipnet -lsipnet_common -lm
 
 # List test files in this directory here
-TEST_CFILES=testNitrogenCycle.c testDependencyFunctions.c testBalance.c testMethane.c testSoilMoisture.c
+TEST_CFILES=testNitrogenCycle.c testDependencyFunctions.c testBalance.c testMethane.c testSoilMoisture.c testWoodStorage.c
 
 # The rest is boilerplate, likely copyable as is to a new test directory
 TEST_OBJ_FILES=$(TEST_CFILES:%.c=%.o)
@@ -40,6 +40,6 @@ $(RUN_EXECUTABLES):
 	./$(basename $@)
 
 clean:
-	rm -f $(TEST_OBJ_FILES) $(TEST_EXECUTABLES) events.in events.out balance.out
+	rm -f $(TEST_OBJ_FILES) $(TEST_EXECUTABLES) events.in events.out balance.out wood_storage.out storage_harvest.in
 
 .PHONY: all tests clean run $(RUN_EXECUTABLES)

tests/sipnet/test_modeling/testWoodStorage.c

@@ -0,0 +1,260 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "common/logging.h"
+#include "utils/tUtils.h"
+
+#include "sipnet/events.c"
+#include "sipnet/sipnet.c"
+
+#define TEST_EPS 1e-8
+#define OUTPUT_FILE "wood_storage.out"
+#define HARVEST_FILE "storage_harvest.in"
+
+static ClimateNode testClimate;
+
+static int checkClose(const char *label, double actual, double expected) {
+  if (fabs(actual - expected) <= TEST_EPS) {
+    return 0;
+  }
+
+  logTest("%s is %f, expected %f\n", label, actual, expected);
+  return 1;
+}
+
+static void initStorageTestState(void) {
+  initContext();
+  updateIntContext("events", 1, CTX_TEST);
+
+  testClimate = (ClimateNode){0};
+  testClimate.year = 2026;
+  testClimate.day = 120;
+  testClimate.time = 12.0;
+  testClimate.length = 1.0;
+  testClimate.tair = 0.0;
+  testClimate.tsoil = 5.0;
+  climate = &testClimate;
+  firstClimate = &testClimate;
+
+  if (meanNPP == NULL) {
+    meanNPP = newMeanTracker(0.0, MEAN_NPP_DAYS, MEAN_NPP_MAX_ENTRIES);
+  } else {
+    resetMeanTracker(meanNPP, 0.0);
+  }
+
+  params.baseVegResp = 1.0;
+  params.vegRespQ10 = 1.0;
+  params.woodTurnoverRate = 0.5;
+  params.coarseRootTurnoverRate = 0.0;
+  params.fineRootTurnoverRate = 0.0;
+  params.coarseRootAllocation = 0.0;
+  params.fineRootAllocation = 0.0;
+  params.woodAllocation = 0.0;
+  params.leafCN = 20.0;
+  params.woodCN = 10.0;
+  params.fineRootCN = 30.0;
+
+  envi = (Envi){0};
+  envi.soilC = 10.0;
+  envi.soilWater = 5.0;
+
+  initTrackers();
+  resetFluxes();
+  initBalanceTracker();
+}
+
+static int getOutputValue(const char *name, double *value) {
+  FILE *out = fopen(OUTPUT_FILE, "w");
+  if (out == NULL) {
+    logTest("Unable to open %s for writing\n", OUTPUT_FILE);
+    return 1;
+  }
+
+  outputHeader(out);
+  outputState(out, testClimate.year, testClimate.day, testClimate.time);
+  fclose(out);
+
+  out = fopen(OUTPUT_FILE, "r");
+  if (out == NULL) {
+    logTest("Unable to open %s for reading\n", OUTPUT_FILE);
+    return 1;
+  }
+
+  char header[4096];
+  char values[4096];
+  if (fgets(header, sizeof(header), out) == NULL ||
+      fgets(values, sizeof(values), out) == NULL) {
+    fclose(out);
+    logTest("Unable to read output lines from %s\n", OUTPUT_FILE);
+    return 1;
+  }
+  fclose(out);
+
+  int targetIndex = -1;
+  int index = 0;
+  char *headerToken = strtok(header, " \t\n");
+  while (headerToken != NULL) {
+    if (strcmp(headerToken, name) == 0) {
+      targetIndex = index;
+      break;
+    }
+    index++;
+    headerToken = strtok(NULL, " \t\n");
+  }
+  if (targetIndex < 0) {
+    logTest("Unable to find output column %s\n", name);
+    return 1;
+  }
+
+  index = 0;
+  char *valueToken = strtok(values, " \t\n");
+  while (valueToken != NULL) {
+    if (index == targetIndex) {
+      *value = atof(valueToken);
+      return 0;
+    }
+    index++;
+    valueToken = strtok(NULL, " \t\n");
+  }
+
+  logTest("Unable to find value for output column %s\n", name);
+  return 1;
+}
+
+static int testOutputSemantics(void) {
+  int status = 0;
+  double reportedWood = 0.0;
+  double storage = 0.0;
+
+  logTest("Starting testOutputSemantics()\n");
+  initStorageTestState();
+
+  envi.plantWoodC = 100.0;
+  envi.plantWoodCStorageDelta = 25.0;
+  envi.plantLeafC = 5.0;
+
+  status |= getOutputValue("plantWoodC", &reportedWood);
+  status |= getOutputValue("nppStorage", &storage);
+  status |= checkClose("reported plantWoodC", reportedWood, 125.0);
+  status |= checkClose("nppStorage", storage, 25.0);
+
+  return status;
+}
+
+static int testNegativeStorageDoesNotReduceWoodFluxes(void) {
+  int status = 0;
+  double folResp = 0.0;
+  double woodResp = 0.0;
+
+  logTest("Starting testNegativeStorageDoesNotReduceWoodFluxes()\n");
+  initStorageTestState();
+
+  envi.plantWoodC = 3.0;
+  envi.plantWoodCStorageDelta = -5.0;
+
+  vegResp(&folResp, &woodResp, 0.0);
+  calcRootAndWoodFluxes();
+
+  status |= checkClose("woodResp", woodResp, 3.0);
+  status |= checkClose("woodLitter", fluxes.woodLitter, 1.5);
+  status |= checkClose("storage-backed wood carbon",
+                       getStorageBackedWoodCarbon(), 3.0);
+
+  return status;
+}
+
+static int writeHarvestEventFile(void) {
+  FILE *events = fopen(HARVEST_FILE, "w");
+  if (events == NULL) {
+    logTest("Unable to open %s for writing\n", HARVEST_FILE);
+    return 1;
+  }
+
+  fprintf(events, "%d %d harv 1 0 0 0\n", testClimate.year, testClimate.day);
+  fclose(events);
+  return 0;
+}
+
+static void runHarvestEvent(void) {
+  initEvents(HARVEST_FILE, "events.out", 0);
+  setupEvents();
+  processEvents();
+  updateBalanceTrackerPreUpdate();
+  updatePoolsForEvents();
+  updateBalanceTrackerPostUpdate();
+  ensureNonNegativeStocks();
+  updateBalanceTrackerPostClamp();
+  checkBalance();
+  closeEventOutFile();
+}
+
+static int testHarvestRemovesPositiveStorage(void) {
+  int status = 0;
+
+  logTest("Starting testHarvestRemovesPositiveStorage()\n");
+  initStorageTestState();
+  status |= writeHarvestEventFile();
+
+  envi.plantWoodC = 10.0;
+  envi.plantWoodCStorageDelta = 5.0;
+  envi.plantLeafC = 2.0;
+
+  runHarvestEvent();
+
+  status |= checkClose("plantWoodC", envi.plantWoodC, 0.0);
+  status |=
+      checkClose("plantWoodCStorageDelta", envi.plantWoodCStorageDelta, 0.0);
+  status |= checkClose("eventWoodC", fluxes.eventWoodC, -10.0);
+  status |=
+      checkClose("eventWoodStorageDelta", fluxes.eventWoodStorageDelta, -5.0);
+  status |= checkClose("eventOutputC", fluxes.eventOutputC, 17.0);
+  status |= checkClose("deltaC", balanceTracker.deltaC, 0.0);
+
+  return status;
+}
+
+static int testHarvestWithNegativeStorageLeavesDebtVisible(void) {
+  int status = 0;
+
+  logTest("Starting testHarvestWithNegativeStorageLeavesDebtVisible()\n");
+  initStorageTestState();
+  status |= writeHarvestEventFile();
+
+  envi.plantWoodC = 10.0;
+  envi.plantWoodCStorageDelta = -15.0;
+  envi.plantLeafC = 2.0;
+
+  runHarvestEvent();
+
+  status |= checkClose("plantWoodC", envi.plantWoodC, 0.0);
+  status |=
+      checkClose("plantWoodCStorageDelta", envi.plantWoodCStorageDelta, -15.0);
+  status |= checkClose("total plantWoodC", getPlantWoodCTotal(), -15.0);
+  status |= checkClose("eventWoodC", fluxes.eventWoodC, -10.0);
+  status |=
+      checkClose("eventWoodStorageDelta", fluxes.eventWoodStorageDelta, 0.0);
+  status |= checkClose("eventOutputC", fluxes.eventOutputC, 12.0);
+  status |= checkClose("deltaC", balanceTracker.deltaC, 0.0);
+
+  return status;
+}
+
+int main(void) {
+  int status = 0;
+
+  logTest("Starting testWoodStorage\n");
+
+  status |= testOutputSemantics();
+  status |= testNegativeStorageDoesNotReduceWoodFluxes();
+  status |= testHarvestRemovesPositiveStorage();
+  status |= testHarvestWithNegativeStorageLeavesDebtVisible();
+
+  if (status) {
+    logTest("FAILED testWoodStorage with status %d\n", status);
+    exit(status);
+  }
+
+  logTest("PASSED testWoodStorage\n");
+}