SIP217 Update harvest and fert event handling

.github/workflows/ci.yml

@@ -10,6 +10,7 @@ on:
       - 'docs/**'
 
 jobs:
+
   # Build and Integration Test
   build:
     strategy:

.github/workflows/debug.yml

@@ -0,0 +1,24 @@
+name: Setup Linux Debug Session
+
+on:
+  workflow_dispatch:
+    inputs:
+      debug_enabled:
+        type: boolean
+        description: 'Run the build with tmate debugging enabled (https://github.com/marketplace/actions/debugging-with-tmate)'
+        required: false
+        default: false
+
+jobs:
+  # Debug environment
+  debug:
+    runs-on: ubuntu-latest
+    environment: debug
+    steps:
+      - uses: actions/checkout@v4
+      - name: Install gdb
+        if: ${{ github.event_name == 'workflow_dispatch' && inputs.debug_enabled }}
+        run: sudo apt install gdb
+      - name: Setup tmate session
+        if: ${{ github.event_name == 'workflow_dispatch' && inputs.debug_enabled }}
+        uses: mxschmitt/action-tmate@v3

src/sipnet/events.c

@@ -303,13 +303,19 @@ void doWriteEventOut(int year, int day, const char *type, int numParams,
   // Standard prefix for all
   fprintf(eventOutFile, "%4d  %3d  %-7s  ", year, day, type);
 
+  // For debugging on linux
+  char *param;
+  double val;
   // Variable output per oneEvent type
   for (ind = 0; ind < numParams - 1; ind++) {
-    fprintf(eventOutFile, "%s=%-.2f,", va_arg(args, char *),
-            va_arg(args, double));
+    // For debugging on linux
+    param = va_arg(args, char *);
+    val = va_arg(args, double);
+    fprintf(eventOutFile, "%s=%-.2f,", param, val);
   }
-  fprintf(eventOutFile, "%s=%-.2f\n", va_arg(args, char *),
-          va_arg(args, double));
+  param = va_arg(args, char *);
+  val = va_arg(args, double);
+  fprintf(eventOutFile, "%s=%-.2f\n", param, val);
 }
 
 void writeEventOut(EventNode *oneEvent, int numParams, ...) {
@@ -350,9 +356,12 @@ void resetEventFluxes(void) {
   fluxes.eventCoarseRootC = 0.0;
   fluxes.eventEvap = 0.0;
   fluxes.eventSoilWater = 0.0;
+  fluxes.eventSoilC = 0.0;
   fluxes.eventLitterC = 0.0;
   fluxes.eventMinN = 0.0;
+  fluxes.eventSoilOrgN = 0.0;
   fluxes.eventLitterN = 0.0;
+
   // mass balance
   fluxes.eventInputC = 0.0;
   fluxes.eventOutputC = 0.0;
@@ -448,7 +457,8 @@ void processEvents(void) {
                       "fluxes.eventInputN", inputN / climLen);
       } break;
       case HARVEST: {
-        // Harvest can both remove biomass and move biomass to the litter pool
+        // Harvest can both remove biomass and move biomass to the soil/litter
+        // pools
         const HarvestParams *harvParams = gEvent->eventParams;
         const double fracRA = harvParams->fractionRemovedAbove;
         const double fracTA = harvParams->fractionTransferredAbove;
@@ -457,8 +467,9 @@ void processEvents(void) {
 
         const double woodC = envi.plantWoodC + envi.plantWoodCStorageDelta;
         // Litter increase
-        const double litterAdd = fracTA * (envi.plantLeafC + woodC) +
-                                 fracTB * (envi.fineRootC + envi.coarseRootC);
+        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);
@@ -468,24 +479,35 @@ void processEvents(void) {
         const double coarseDelta = -envi.coarseRootC * (fracRB + fracTB);
 
         // Pool updates:
+        if (!ctx.litterPool) {
+          // send it all to the soil
+          soilAdd += litterAdd;
+          litterAdd = 0.0;
+        }
         fluxes.eventLitterC += litterAdd / climLen;
+        fluxes.eventSoilC += soilAdd / climLen;
         fluxes.eventLeafC += leafDelta / climLen;
         fluxes.eventWoodC += woodDelta / climLen;
         fluxes.eventFineRootC += fineDelta / climLen;
         fluxes.eventCoarseRootC += coarseDelta / climLen;
 
         // No need to allocate to biomass N pools, we don't track that N
-        // explicitly. We do need to handle litter N, though.
+        // explicitly. We do need to handle soil and litter N, though.
+        // Note: ctx.nitrogenCycle implies ctx.litterPool
         // Litter N increase
         double litterNAdd = 0.0;
+        double soilNAdd = 0.0;
         if (ctx.nitrogenCycle) {
           const double totalAbove = (envi.plantLeafC / params.leafCN) +
                                     (envi.plantWoodC / params.woodCN);
           const double totalBelow = (envi.fineRootC / params.fineRootCN) +
                                     (envi.coarseRootC / params.woodCN);
-          litterNAdd = (fracTA * totalAbove) + (fracTB * totalBelow);
+          litterNAdd = fracTA * totalAbove;
+          soilNAdd = fracTB * totalBelow;
+          fluxes.eventSoilOrgN += soilNAdd / climLen;
           fluxes.eventLitterN += litterNAdd / climLen;
         }
+
         // MASS BALANCE: removed fractions are system outputs
         const double outputC = ((woodC + envi.plantLeafC) * fracRA +
                                 (envi.fineRootC + envi.coarseRootC) * fracRB);
@@ -501,14 +523,15 @@ void processEvents(void) {
                          fracRB);
           fluxes.eventOutputN += outputN / climLen;
         }
-
         writeEventOut(
-            gEvent, 8, "fluxes.eventLitterC", litterAdd / climLen,
-            "fluxes.eventLeafC", leafDelta / climLen, "fluxes.eventWoodC",
-            woodDelta / climLen, "fluxes.eventFineRootC", fineDelta / climLen,
+            gEvent, 10, "fluxes.eventSoilC", soilAdd / climLen,
+            "fluxes.eventLitterC", litterAdd / climLen, "fluxes.eventLeafC",
+            leafDelta / climLen, "fluxes.eventWoodC", woodDelta / climLen,
+            "fluxes.eventFineRootC", fineDelta / climLen,
             "fluxes.eventCoarseRootC", coarseDelta / climLen,
-            "fluxes.eventLitterN", litterNAdd / climLen, "fluxes.eventOutputC",
-            outputC / climLen, "fluxes.eventOutputN", outputN / climLen);
+            "fluxes.eventSoilOrgN", soilNAdd / climLen, "fluxes.eventLitterN",
+            litterNAdd / climLen, "fluxes.eventOutputC", outputC / climLen,
+            "fluxes.eventOutputN", outputN / climLen);
       } break;
       case TILLAGE: {
         // BIG NOTE: this is the one event type that is NOT modeled as a flux;
@@ -526,19 +549,25 @@ void processEvents(void) {
         const double orgC = fertParams->orgC;
         double orgN = fertParams->orgN;
         double minN = fertParams->minN;
-        if (!ctx.nitrogenCycle) {
+        if (ctx.litterPool) {
+          fluxes.eventLitterC += orgC / climLen;
+        } else {
+          fluxes.eventSoilC += orgC / climLen;
+        }
+
+        if (ctx.nitrogenCycle) {
           // As the warning says in readEventData(), we ignore N when the
           // nitrogen cycle model is off
-          orgN = 0.0;
-          minN = 0.0;
+          // Implies ctx.litterPool
+          fluxes.eventLitterN += orgN / climLen;
+          fluxes.eventMinN += minN / climLen;
         }
-        fluxes.eventLitterC += orgC / climLen;
-        fluxes.eventLitterN += orgN / climLen;
-        fluxes.eventMinN += minN / climLen;
 
         // MASS BALANCE: this is a system input
         fluxes.eventInputC += orgC / climLen;
-        fluxes.eventInputN += (orgN + minN) / climLen;
+        if (ctx.nitrogenCycle) {
+          fluxes.eventInputN += (orgN + minN) / climLen;
+        }
 
         writeEventOut(gEvent, 5, "fluxes.eventOrgN", orgN / climLen,
                       "fluxes.eventLitterC", orgC / climLen, "fluxes.eventMinN",
@@ -561,10 +590,9 @@ void updatePoolsForEvents(void) {
   envi.plantLeafC += fluxes.eventLeafC * climate->length;
 
   // Harvest and fertilization events
+  envi.soilC += fluxes.eventSoilC * climate->length;
   if (ctx.litterPool) {
     envi.litterC += fluxes.eventLitterC * climate->length;
-  } else {
-    envi.soilC += fluxes.eventLitterC * climate->length;
   }
 
   // Harvest and planting events
@@ -576,11 +604,12 @@ void updatePoolsForEvents(void) {
   envi.soilWater += fluxes.eventSoilWater * climate->length;
 
   // NITROGEN
-  // Fertilization events
-  // (Harvest and planting events TBD)
-  // Note: litter_pool is required for nitrogen_cycle
+  // Harvest and fertilization events
+  // (Planting events don't explicitly handle N)
+  // Note: nitrogen_cycle implies litter_pool
   if (ctx.nitrogenCycle) {
     envi.minN += fluxes.eventMinN * climate->length;
+    envi.soilOrgN += fluxes.eventSoilOrgN * climate->length;
     envi.litterN += fluxes.eventLitterN * climate->length;
   }
 }

src/sipnet/state.h

@@ -565,10 +565,14 @@ typedef struct FluxVars {
   double eventEvap;
   // irrigation water that goes to the soil
   double eventSoilWater;
-  // carbon added to litter pool (if used) or soil pool (if not)
+  // carbon added to soil pool
+  double eventSoilC;
+  // carbon added to litter pool
   double eventLitterC;
   // nitrogen added to soil mineral N pool
   double eventMinN;
+  // nitrogen added to soil org N pool
+  double eventSoilOrgN;
   // nitrogen added to litter N pool
   double eventLitterN;
   // MASS BALANCE HELPERS

tests/sipnet/test_events_infrastructure/events_output_header.out

@@ -1,7 +1,7 @@
 year  day  type     param_name=delta[,param_name=delta,...]
 2023   65  plant    fluxes.eventLeafC=6.00,fluxes.eventWoodC=8.00,fluxes.eventFineRootC=10.00,fluxes.eventCoarseRootC=12.00,fluxes.eventInputC=36.00,fluxes.eventInputN=0.00
 2023   70  irrig    fluxes.eventSoilWater=10.00,fluxes.eventEvap=0.00
-2023  200  harv     fluxes.eventLitterC=10.93,fluxes.eventLeafC=-11.86,fluxes.eventWoodC=-9.50,fluxes.eventFineRootC=-7.46,fluxes.eventCoarseRootC=-7.78,fluxes.eventLitterN=0.00,fluxes.eventOutputC=25.66,fluxes.eventOutputN=0.00
+2023  200  harv     fluxes.eventSoilC=3.81,fluxes.eventLitterC=7.12,fluxes.eventLeafC=-11.86,fluxes.eventWoodC=-9.50,fluxes.eventFineRootC=-7.46,fluxes.eventCoarseRootC=-7.78,fluxes.eventSoilOrgN=0.00,fluxes.eventLitterN=0.00,fluxes.eventOutputC=25.66,fluxes.eventOutputN=0.00
 2024   65  plant    fluxes.eventLeafC=6.00,fluxes.eventWoodC=10.00,fluxes.eventFineRootC=14.00,fluxes.eventCoarseRootC=18.00,fluxes.eventInputC=48.00,fluxes.eventInputN=0.00
 2024   70  irrig    fluxes.eventSoilWater=5.00,fluxes.eventEvap=5.00
-2024  200  harv     fluxes.eventLitterC=8.51,fluxes.eventLeafC=-2.78,fluxes.eventWoodC=-3.27,fluxes.eventFineRootC=-5.04,fluxes.eventCoarseRootC=-5.93,fluxes.eventLitterN=0.00,fluxes.eventOutputC=8.51,fluxes.eventOutputN=0.00
+2024  200  harv     fluxes.eventSoilC=5.48,fluxes.eventLitterC=3.02,fluxes.eventLeafC=-2.78,fluxes.eventWoodC=-3.27,fluxes.eventFineRootC=-5.04,fluxes.eventCoarseRootC=-5.93,fluxes.eventSoilOrgN=0.00,fluxes.eventLitterN=0.00,fluxes.eventOutputC=8.51,fluxes.eventOutputN=0.00

tests/sipnet/test_events_infrastructure/events_output_no_header.out

@@ -1,6 +1,6 @@
 2023   65  plant    fluxes.eventLeafC=20.00,fluxes.eventWoodC=10.00,fluxes.eventFineRootC=8.00,fluxes.eventCoarseRootC=6.00,fluxes.eventInputC=44.00,fluxes.eventInputN=0.00
 2023   70  irrig    fluxes.eventSoilWater=10.00,fluxes.eventEvap=0.00
-2023  200  harv     fluxes.eventLitterC=24.80,fluxes.eventLeafC=-9.60,fluxes.eventWoodC=-6.40,fluxes.eventFineRootC=-9.60,fluxes.eventCoarseRootC=-9.60,fluxes.eventLitterN=0.00,fluxes.eventOutputC=10.40,fluxes.eventOutputN=0.00
+2023  200  harv     fluxes.eventSoilC=24.80,fluxes.eventLitterC=0.00,fluxes.eventLeafC=-9.60,fluxes.eventWoodC=-6.40,fluxes.eventFineRootC=-9.60,fluxes.eventCoarseRootC=-9.60,fluxes.eventSoilOrgN=0.00,fluxes.eventLitterN=0.00,fluxes.eventOutputC=10.40,fluxes.eventOutputN=0.00
 2024   65  plant    fluxes.eventLeafC=20.00,fluxes.eventWoodC=10.00,fluxes.eventFineRootC=8.00,fluxes.eventCoarseRootC=6.00,fluxes.eventInputC=44.00,fluxes.eventInputN=0.00
 2024   70  irrig    fluxes.eventSoilWater=5.00,fluxes.eventEvap=5.00
-2024  200  harv     fluxes.eventLitterC=24.28,fluxes.eventLeafC=-20.64,fluxes.eventWoodC=-11.76,fluxes.eventFineRootC=-5.76,fluxes.eventCoarseRootC=-4.96,fluxes.eventLitterN=0.00,fluxes.eventOutputC=18.84,fluxes.eventOutputN=0.00
+2024  200  harv     fluxes.eventSoilC=24.28,fluxes.eventLitterC=0.00,fluxes.eventLeafC=-20.64,fluxes.eventWoodC=-11.76,fluxes.eventFineRootC=-5.76,fluxes.eventCoarseRootC=-4.96,fluxes.eventSoilOrgN=0.00,fluxes.eventLitterN=0.00,fluxes.eventOutputC=18.84,fluxes.eventOutputN=0.00

tests/sipnet/test_events_types/events_two_harvest.in

@@ -1,2 +1,2 @@
 2024  70 harv   0.1 0.2 0.3 0.4
-2024  70 harv   0.25 0.25 0.25 0.25
+2024  70 harv   0.2 0.1 0.2 0.1

tests/sipnet/test_events_types/testEventFertilization.c

@@ -61,6 +61,7 @@ int run(void) {
   initEvents("events_one_fert.in", 0);
   setupEvents();
   procEvents();
+  closeEventOutFile();
 
   // First fert: (15-5-10)
   expOrgN = 0;  // nitrogen cycle is off
@@ -78,6 +79,7 @@ int run(void) {
   initEvents("events_two_fert.in", 1);
   setupEvents();
   procEvents();
+  closeEventOutFile();
 
   // First event same as above (15-5-10)
   expOrgN = 3 + 15;