Rename nppStorage to woodCStorageDelta for clarity

.gitignore

@@ -9,6 +9,8 @@ transpose
 .mkdocs.stamp
 site/*
 site_preview/*
+_codeql_build_dir/
+_codeql_detected_source_root
 
 # documentation
 docs/api

CMakeLists.txt

@@ -37,6 +37,7 @@ add_library(sipnetlib
         src/sipnet/runmean.c
         src/sipnet/sipnet.c
         src/sipnet/state.c
+        src/sipnet/balance.c
 )
 
 add_library(tests
@@ -52,4 +53,6 @@ add_library(tests
         tests/sipnet/test_events_types/testEventTillage.c
         tests/sipnet/test_bugfixes/testEventFileOrderChecks.c
         tests/sipnet/test_modeling/testNitrogenCycle.c
+        tests/sipnet/test_modeling/testDependencyFunctions.c
+        tests/sipnet/test_modeling/testBalance.c
 )

Makefile

@@ -11,7 +11,7 @@ COMMON_CFILES:=context.c logging.c modelParams.c util.c
 COMMON_CFILES:=$(addprefix src/common/, $(COMMON_CFILES))
 COMMON_OFILES=$(COMMON_CFILES:.c=.o)
 
-SIPNET_CFILES:=sipnet.c cli.c events.c frontend.c outputItems.c runmean.c state.c
+SIPNET_CFILES:=sipnet.c cli.c events.c frontend.c outputItems.c runmean.c state.c balance.c
 SIPNET_CFILES:=$(addprefix src/sipnet/, $(SIPNET_CFILES))
 SIPNET_OFILES=$(SIPNET_CFILES:.c=.o)
 SIPNET_LIBS=-lsipnet_common

src/sipnet/balance.c

@@ -0,0 +1,101 @@
+#include "balance.h"
+#include "common/context.h"
+#include "state.h"
+
+// Definition of global balance tracker struct
+BalanceTracker balanceTracker;
+
+void getMassTotals(double *carbon, double *nitrogen) {
+  *carbon = envi.plantWoodC + envi.plantLeafC + envi.fineRootC +
+            envi.coarseRootC + envi.soilC + envi.woodCStorageDelta;
+  if (ctx.litterPool) {
+    *carbon += envi.litterC;
+  }
+
+  if (ctx.nitrogenCycle) {
+    *nitrogen =
+        envi.plantWoodC / params.woodCN + envi.plantLeafC / params.leafCN +
+        envi.fineRootC / params.fineRootCN + envi.coarseRootC / params.woodCN +
+        envi.soilOrgN + envi.litterN + envi.minN;
+  } else {
+    *nitrogen = 0.0;
+  }
+}
+
+void updateBalanceTrackerPreUpdate(void) {
+  // Set the pre-update pool totals
+  getMassTotals(&balanceTracker.preTotalC, &balanceTracker.preTotalN);
+}
+
+void updateBalanceTrackerPostUpdate(void) {
+  // Set the post-update pool totals
+  getMassTotals(&balanceTracker.postTotalC, &balanceTracker.postTotalN);
+
+  // Calculate the system inputs and outputs
+  // CARBON
+  balanceTracker.inputsC = fluxes.photosynthesis +  // GPP
+                           fluxes.eventInputC;  // agro event additions
+  balanceTracker.outputsC = fluxes.rVeg + fluxes.rFineRoot +
+                            fluxes.rCoarseRoot +  // R_a
+                            fluxes.rSoil +  // R_h
+                            fluxes.eventOutputC;  // agro event removals
+  if (ctx.litterPool) {
+    balanceTracker.outputsC += fluxes.rLitter;
+  }
+  // Account for climate length
+  balanceTracker.inputsC *= climate->length;
+  balanceTracker.outputsC *= climate->length;
+
+  // NITROGEN
+  if (ctx.nitrogenCycle) {
+    balanceTracker.inputsN =
+        // TODO: fluxes.fixation +
+        fluxes.eventInputN;
+    balanceTracker.outputsN =
+        fluxes.nLeaching + fluxes.nVolatilization + fluxes.eventOutputN;
+
+    // Account for climate length
+    balanceTracker.inputsN *= climate->length;
+    balanceTracker.outputsN *= climate->length;
+  }
+}
+
+void initBalanceTracker(void) {
+  // Initialize all to zero
+  balanceTracker.preTotalC = 0.0;
+  balanceTracker.postTotalC = 0.0;
+  balanceTracker.inputsC = 0.0;
+  balanceTracker.outputsC = 0.0;
+  balanceTracker.preTotalN = 0.0;
+  balanceTracker.postTotalN = 0.0;
+  balanceTracker.inputsN = 0.0;
+  balanceTracker.outputsN = 0.0;
+  balanceTracker.deltaC = 0.0;
+  balanceTracker.deltaN = 0.0;
+}
+
+void checkBalance(void) {
+  // CARBON
+  // Pool delta
+  double poolCDelta = balanceTracker.postTotalC - balanceTracker.preTotalC;
+  // System delta
+  double systemCDelta = balanceTracker.inputsC - balanceTracker.outputsC;
+  balanceTracker.deltaC = poolCDelta - systemCDelta;
+
+  // NITROGEN
+  // Pool delta
+  double poolNDelta = balanceTracker.postTotalN - balanceTracker.preTotalN;
+  // System delta
+  double systemNDelta = balanceTracker.inputsN - balanceTracker.outputsN;
+  balanceTracker.deltaN = poolNDelta - systemNDelta;
+
+  // To avoid weird negative-zero issues...
+  if (balanceTracker.deltaC < 1e-8) {
+    balanceTracker.deltaC = 0.0;
+  }
+  if (balanceTracker.deltaN < 1e-8) {
+    balanceTracker.deltaN = 0.0;
+  }
+
+  // TBD: warn if balance off?
+}

src/sipnet/balance.h

@@ -0,0 +1,41 @@
+#ifndef BALANCE_H
+#define BALANCE_H
+
+typedef struct BalanceTrackerStruct {
+  // Mass balance checks:
+  //   X_t - X_(t-1) = inputs - outputs + tolerance
+  // or, we check that
+  //   (pool delta) + (system delta) < tolerance
+  // where
+  //   pool delta = X_t - X_(t-1)
+  //   system delta = outputs - inputs
+
+  // Carbon balance
+  double preTotalC;
+  double postTotalC;
+  double inputsC;
+  double outputsC;
+
+  // Nitrogen balance
+  double preTotalN;
+  double postTotalN;
+  double inputsN;
+  double outputsN;
+
+  // Checks
+  double deltaC;
+  double deltaN;
+} BalanceTracker;
+
+// Global var
+extern BalanceTracker balanceTracker;
+
+void updateBalanceTrackerPreUpdate(void);
+
+void updateBalanceTrackerPostUpdate(void);
+
+void initBalanceTracker(void);
+
+void checkBalance(void);
+
+#endif  // BALANCE_H

src/sipnet/cli.c

@@ -89,7 +89,7 @@ void usage(char *progName) {
   printf("\n");
   printf("Output flags: (prepend flag with 'no-' to force off, eg '--no-print-header')\n");
   printf("  --do-main-output     Print time series of all output variables to <file-name>.out (1)\n");
-  printf("  --do-single-outputs  Print outputs one variable per file (e.g. <file-name>.NEE)\n");
+  printf("  --do-single-outputs  Print selection* of outputs one variable per file (e.g. <file-name>.NEE)\n");
   printf("  --dump-config        Print final config to <file-name>.config (0)\n");
   printf("  --print-header       Whether to print header row in output files (1)\n");
   printf("  --quiet              Suppress info and warning message (0)\n");
@@ -98,6 +98,8 @@ void usage(char *progName) {
   printf("  -h, --help           Print this message and exit\n");
   printf("  -v, --version        Print version information and exit\n");
   printf("\n");
+  printf("*do-single-outputs option outputs are: NEE, NEE_cum, GPP, GPP_cum\n");
+  printf("\n");
   printf("Configuration options are read from <input_file>. Other options specified on the command\n");
   printf("line override settings from that file.\n");
   printf("\n");

src/sipnet/events.c

@@ -192,7 +192,7 @@ EventNode *readEventData(char *eventFile) {
   if (access(eventFile, F_OK) != 0) {
     // no file found, which is fine; we're done, a vector of NULL is what we
     // want for newEvents
-    logInfo("No event file found, assuming no newEvents\n");
+    logInfo("No event file found, assuming no input events\n");
     return newEvents;
   }
 
@@ -271,24 +271,38 @@ void openEventOutFile(int printHeader) {
   }
 }
 
-void writeEventOut(EventNode *oneEvent, int numParams, ...) {
-  va_list args;
+void doWriteEventOut(int year, int day, const char *type, int numParams,
+                     va_list args) {
   int ind = 0;
 
   // Spec:
   // year day event_type <param name=delta>[,<param name>=<delta>,...]
 
   // Standard prefix for all
-  fprintf(eventOutFile, "%4d  %3d  %-7s  ", oneEvent->year, oneEvent->day,
-          eventTypeToString(oneEvent->type));
+  fprintf(eventOutFile, "%4d  %3d  %-7s  ", year, day, type);
+
   // Variable output per oneEvent type
-  va_start(args, numParams);
   for (ind = 0; ind < numParams - 1; ind++) {
     fprintf(eventOutFile, "%s=%-.2f,", va_arg(args, char *),
             va_arg(args, double));
   }
   fprintf(eventOutFile, "%s=%-.2f\n", va_arg(args, char *),
           va_arg(args, double));
+}
+
+void writeEventOut(EventNode *oneEvent, int numParams, ...) {
+  va_list args;
+  va_start(args, numParams);
+  doWriteEventOut(oneEvent->year, oneEvent->day,
+                  eventTypeToString(oneEvent->type), numParams, args);
+  va_end(args);
+}
+
+void writeComputedEventOut(int year, int day, const char *type, int numParams,
+                           ...) {
+  va_list args;
+  va_start(args, numParams);
+  doWriteEventOut(year, day, type, numParams, args);
   va_end(args);
 }
 
@@ -317,13 +331,14 @@ void resetEventFluxes(void) {
   fluxes.eventLitterC = 0.0;
   fluxes.eventMinN = 0.0;
   fluxes.eventOrgN = 0.0;
+  // mass balance
+  fluxes.eventInputC = 0.0;
+  fluxes.eventOutputC = 0.0;
+  fluxes.eventInputN = 0.0;
+  fluxes.eventOutputN = 0.0;
 }
 
 void processEvents(void) {
-  // This should be in events.h/c, but with all the global state defined in
-  // this file, let's leave it here for now. Maybe someday we will factor that
-  // out.
-
   // Set all event fluxes to zero, as these have no memory from one step to
   // the next
   resetEventFluxes();
@@ -390,7 +405,18 @@ void processEvents(void) {
         fluxes.eventFineRootC += fineRootC / climLen;
         fluxes.eventCoarseRootC += coarseRootC / climLen;
 
-        // FUTURE: allocate to N pools
+        // No need to allocate to biomass N pools, we don't track that N
+        // explicitly
+
+        // MASS BALANCE: this is a system input
+        fluxes.eventInputC += leafC + woodC + fineRootC + coarseRootC;
+        fluxes.eventInputC /= climLen;
+        if (ctx.nitrogenCycle) {
+          fluxes.eventInputN += leafC / params.leafCN + woodC / params.woodCN +
+                                fineRootC / params.fineRootCN +
+                                coarseRootC / params.woodCN;
+          fluxes.eventInputN /= climLen;
+        }
 
         writeEventOut(gEvent, 4, "fluxes.eventLeafC", leafC / climLen,
                       "fluxes.eventWoodC", woodC / climLen,
@@ -423,7 +449,24 @@ void processEvents(void) {
         fluxes.eventFineRootC += fineDelta / climLen;
         fluxes.eventCoarseRootC += coarseDelta / climLen;
 
-        // FUTURE: move/remove biomass in N pools
+        // No need to allocate to biomass N pools, we don't track that N
+        // explicitly. We do need to handle litter N, though.
+        // FUTURE (ticket #244): move/remove biomass in litter N
+
+        // MASS BALANCE: removed fractions are system outputs
+        fluxes.eventOutputC += (envi.plantWoodC + envi.plantLeafC) * fracRA +
+                               (envi.fineRootC + envi.coarseRootC) * fracRB;
+        fluxes.eventOutputC /= climLen;
+        if (ctx.nitrogenCycle) {
+          fluxes.eventOutputN += (envi.plantWoodC / params.woodCN +
+                                  envi.plantLeafC / params.leafCN) *
+                                     fracRA +
+                                 (envi.fineRootC / params.fineRootCN +
+                                  envi.coarseRootC / params.woodCN) *
+                                     fracRB;
+          fluxes.eventOutputN /= climLen;
+        }
+
         writeEventOut(gEvent, 5, "fluxes.eventLitterC", litterAdd / climLen,
                       "fluxes.eventLeafC", leafDelta / climLen,
                       "fluxes.eventWoodC", woodDelta / climLen,
@@ -456,6 +499,12 @@ void processEvents(void) {
           // nitrogen cycle model is off, so no update needed
         }
 
+        // MASS BALANCE: this is a system input
+        fluxes.eventInputC += orgC;
+        if (ctx.nitrogenCycle) {
+          fluxes.eventInputN += orgN + minN;
+        }
+
         writeEventOut(gEvent, 3, "fluxes.eventOrgN", orgN / climLen,
                       "fluxes.eventLitterC", orgC / climLen, "fluxes.eventMinN",
                       minN / climLen);

src/sipnet/events.h

@@ -117,14 +117,33 @@ void openEventOutFile(int printHeader);
  *
  * year day event_type \<param_name>=\<delta>[,\<param_name>=\<delta>,...]
  *
- * \param out       FILE pointer to events.out
  * \param oneEvent     Pointer to oneEvent node
- * \param numParams Number of param/value pairs to write
+ * \param numParams Number of param/value PAIRS to write
  * \param ...       Pairs of (char*, double) arguments to write, 2*numParams
  *                  values
  */
 void writeEventOut(EventNode *oneEvent, int numParams, ...);
 
+/*!
+ * \brief Write a line to events.out for a 'computed' event
+ *
+ * Same as writeEventOut, but for events that are computed internally, such
+ * as leaf on/leaf off events.
+ *
+ * Output format:
+ *
+ * year day event_type \<param_name>=\<delta>[,\<param_name>=\<delta>,...]
+ *
+ * \param year Year of event
+ * \param day  Day of event
+ * \param type Type of event
+ * \param numParams Number of param/value PAIRS to write
+ * \param ...       Pairs of (char*, double) arguments to write, 2*numParams
+ *                  values
+ */
+void writeComputedEventOut(int year, int day, const char *type, int numParams,
+                           ...);
+
 /*!
  * Close the event output file
  * @param file FILE pointer for output file