SIP208, SIP209 Add carbon balance check (#248)

* Add balance.c|h

* Improve do-single-output help

* Carbon and nitrogen balance checks

* Add event fluxes to track system inputs and outputs; add nppStorage pool to envi

* Add external writeEvent function, refactor original

* Add event for leaf on/off; split woodC into two pools; add balance checks

* Updates for balance checks and leaf on/off events

* Add ability to compare to base

* Remove tolerance from balance tracker

* Update col list for niwot

* Avoid negative zero

* Remove negative zeros

* Add test for balance check

* PR feedback

* Turn off microbes, as it is now broken

* Test updates after balance check updates

* Split leaf-on carbon creation into own flux

* Fix leafOn flux

* Updates for leafOn split

* Force skipping russell_4

* Update for new event fluxes

Author: Alomir

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,124 @@
+#include "balance.h"
+
+#include <math.h>
+
+#include "common/context.h"
+#include "common/exitCodes.h"
+#include "common/logging.h"
+#include "state.h"
+
+// Definition of global balance tracker struct
+BalanceTracker balanceTracker;
+
+void getMassTotals(double *carbon, double *nitrogen) {
+  *carbon = (envi.plantWoodC + envi.plantWoodCStorageDelta) + envi.plantLeafC +
+            envi.fineRootC + envi.coarseRootC + envi.soilC;
+  if (ctx.litterPool) {
+    *carbon += envi.litterC;
+  }
+
+  if (ctx.nitrogenCycle) {
+    // Note: this is the one place where we ust plantWoodC by itself; it's the
+    // reason plantNSCWoodCDelta was created, so that we can ignore it here.
+    *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 (fabs(balanceTracker.deltaC) < 1e-8) {
+    balanceTracker.deltaC = 0.0;
+  }
+  if (fabs(balanceTracker.deltaN) < 1e-8) {
+    balanceTracker.deltaN = 0.0;
+  }
+
+  int err = 0;
+  if (fabs(balanceTracker.deltaC) > 0.0) {
+    err = 1;
+    logInternalError(
+        "Carbon balance check failed (delta=%8.4f, Y: %d D: %d T: %4.2f)\n",
+        balanceTracker.deltaC, climate->year, climate->day, climate->time);
+  }
+  // RE-ENABLE WHEN N IS FIXED
+  // if (fabs(balanceTracker.deltaN) > 0.0) {
+  //  err = 1;
+  //  logInternalError("Nitrogen balance check failed (delta=%8.4f)\n",
+  //                   balanceTracker.deltaN);
+  //}
+  if (err) {
+    logInternalError("Exiting\n");
+    //  exit(EXIT_CODE_INTERNAL_ERROR);
+  }
+}

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");