Revamp harvest and leaf event handling

src/sipnet/debug_log.c

@@ -13,7 +13,7 @@
 
 #define NUM_LOGGED_ENVI_FIELDS 13
 #define NUM_LOGGED_FLUX_FIELDS 62
-#define NUM_LOGGED_TRACKER_FIELDS 33
+#define NUM_LOGGED_TRACKER_FIELDS 32
 #define NUM_LOGGED_PHEN_TRACKER_FIELDS 3
 #define NUM_LOGGED_SURVIVAL_FIELDS 1
 #define NUM_LOGGED_EVENT_TRACKER_FIELDS 8
@@ -28,7 +28,7 @@
 
 _Static_assert(sizeof(Envi) == DEBUG_LAYOUT_ENVI_SIZE,
                "Debug log schema drift: Envi changed; update debug_log.c");
 _Static_assert(sizeof(Fluxes) == DEBUG_LAYOUT_FLUX_SIZE,
                "Debug log schema drift: Fluxes changed; update debug_log.c");
 _Static_assert(sizeof(Trackers) == DEBUG_LAYOUT_TRACKER_SIZE,
                "Debug log schema drift: Trackers changed; update debug_log.c");
@@ -39,7 +39,7 @@
     sizeof(PlantSurvivalTracker) == DEBUG_LAYOUT_SURVIVAL_SIZE,
     "Debug log schema drift: PlantSurvival changed; update debug_log.c");
 _Static_assert(
     sizeof(EventTrackers) == DEBUG_LAYOUT_EVENT_SIZE,
     "Debug log schema drift: EventTrackers changed; update debug_log.c");
 
 typedef enum DebugFieldType {
@@ -120,8 +120,8 @@
   debugFields->fluxDF[ind++] = (DebugField){"rFineRoot", DEBUG_FIELD_DOUBLE, &fluxes.rFineRoot},
   debugFields->fluxDF[ind++] = (DebugField){"leafCreation", DEBUG_FIELD_DOUBLE, &fluxes.leafCreation},
   debugFields->fluxDF[ind++] = (DebugField){"woodCreation", DEBUG_FIELD_DOUBLE, &fluxes.woodCreation},
   debugFields->fluxDF[ind++] = (DebugField){"leafOnCreation", DEBUG_FIELD_DOUBLE, &fluxes.leafOnCreation},
   debugFields->fluxDF[ind++] = (DebugField){"leafOnCreationFromWood", DEBUG_FIELD_DOUBLE, &fluxes.leafOnCreationFromWood},
   debugFields->fluxDF[ind++] = (DebugField){"leafOffLitter", DEBUG_FIELD_DOUBLE, &fluxes.leafOffLitter},
   debugFields->fluxDF[ind++] = (DebugField){"nVolatilization", DEBUG_FIELD_DOUBLE, &fluxes.nVolatilization},
   debugFields->fluxDF[ind++] = (DebugField){"nLeaching", DEBUG_FIELD_DOUBLE, &fluxes.nLeaching},
@@ -185,7 +185,6 @@
   debugFields->trackerDF[ind++] = (DebugField){"yearlyRh", DEBUG_FIELD_DOUBLE, &trackers.yearlyRh},
   debugFields->trackerDF[ind++] = (DebugField){"yearlyNpp", DEBUG_FIELD_DOUBLE, &trackers.yearlyNpp},
   debugFields->trackerDF[ind++] = (DebugField){"yearlyNee", DEBUG_FIELD_DOUBLE, &trackers.yearlyNee},
-  debugFields->trackerDF[ind++] = (DebugField){"yearlyLitter", DEBUG_FIELD_DOUBLE, &trackers.yearlyLitter},
   debugFields->trackerDF[ind++] = (DebugField){"totGpp", DEBUG_FIELD_DOUBLE, &trackers.totGpp},
   debugFields->trackerDF[ind++] = (DebugField){"totRtot", DEBUG_FIELD_DOUBLE, &trackers.totRtot},
   debugFields->trackerDF[ind++] = (DebugField){"totRa", DEBUG_FIELD_DOUBLE, &trackers.totRa},

src/sipnet/events.h

@@ -207,6 +207,11 @@ void processEventsForCarbon(EventNode *event);
  */
 void processEventsForNitrogen(EventNode *event);
 
+/*!
+ * \brief Process harvest and leaf on/off events
+ */
+void processPendingEvents(void);
+
 /*!
  * Update relevant environment pools after event fluxes have been calculated
  */
@@ -225,6 +230,7 @@ typedef struct HarvestTrackersStruct {
   double totalFracTransferredAbove;
   double totalFracTransferredBelow;
   double leafCAdjustment;
+  int harvestOccurred;
 } HarvestTrackers;
 
 // Variables to track events with lingering effects
@@ -233,7 +239,13 @@ typedef struct EventTrackerStruct {
   // equal to exp(-delta_t / 30)
   double d_till_mod;
 
-  // Fraction removed and transferred from harvest event this time step
+  // Total leaf on amount
+  double totalLeafOnAmount;
+
+  // Total leaf off fraction
+  double totalLeafOffFraction;
+
+  // Fraction removed and transferred from harvest events this time step
   HarvestTrackers harvestTrackers;
 } EventTrackers;
 

src/sipnet/limitations.c

@@ -11,27 +11,21 @@
 #include "state.h"
 
 // See limitations.h
-void checkLeafOnLimitation(double *leafOnFlux) {
+void checkLeafOnLimitation(double *leafOnCDemand) {
   // Leaf on events are limited by:
   // * leafGrowth parameter (input leafOn value)
   // * available carbon
   // * available nitrogen
 
-  // Make sure to include current leaf on flux, in case there is more than
-  // one leaf on event
-  double currentFlux = fluxes.leafOnCreation + fluxes.eventLeafOnCreation;
-  double totalDemandFlux = *leafOnFlux + currentFlux;
-  double leafOnCDemand = totalDemandFlux * climate->length;
-
-  if (leafOnCDemand < TINY) {
+  if (*leafOnCDemand < TINY) {
     // Nothing to check
     return;
   }
 
   // First up, carbon. We do not draw from the C storage pool for this.
   double availableC =
       (envi.plantWoodC + envi.coarseRootC) * params.leafOnReallocFrac;
-  double cLimiter = availableC / leafOnCDemand;
+  double cLimiter = availableC / *leafOnCDemand;
 
   double leafOnNDemand = 0.0;
   double nLimiter = 1.0;
@@ -41,7 +35,7 @@ void checkLeafOnLimitation(double *leafOnFlux) {
     // Needed N for this transfer is (what leaves need) - (what wood provides)
     // Reminder: both wood and coarseRoot use params.woodCN, so no need to
     // treat those C demands separately
-    leafOnNDemand = calcLeafOnNFromC(leafOnCDemand);
+    leafOnNDemand = calcLeafOnNFromC(*leafOnCDemand);
     availableN = envi.plantStorageN;
     if (leafOnNDemand > TINY) {
       nLimiter = availableN / leafOnNDemand;
@@ -50,20 +44,19 @@ void checkLeafOnLimitation(double *leafOnFlux) {
   double limitation = unitClip(fmin(cLimiter, nLimiter));
 
   if (limitation < 1) {
-    totalDemandFlux *= limitation;
-    *leafOnFlux = totalDemandFlux - currentFlux;
+    *leafOnCDemand *= limitation;
     if (ctx.nitrogenCycle) {
       logInfo("Leaf on creation %.4f C / %.4f N exceeds available C %.4f / N "
               "%.4f (C ratio: %.4f, N ratio: %.4f), "
               "reducing leaf-on growth by %.2f%% on year %d day %d time %.3f\n",
-              leafOnCDemand, leafOnNDemand, availableC, availableN, cLimiter,
+              *leafOnCDemand, leafOnNDemand, availableC, availableN, cLimiter,
               nLimiter, (1 - limitation) * 100, climate->year, climate->day,
               climate->time);
     } else {
       logInfo("Leaf on creation %.4f exceeds available C %.4f "
               "(C ratio: %.4f, N ratio: %.4f), "
               "reducing leaf-on growth by %.2f%% on year %d day %d time %.3f\n",
-              leafOnCDemand, availableC, cLimiter, nLimiter,
+              *leafOnCDemand, availableC, cLimiter, nLimiter,
               (1 - limitation) * 100, climate->year, climate->day,
               climate->time);
     }

src/sipnet/limitations.h

@@ -12,12 +12,12 @@ void checkLimitations(void);
 /**
  * Check for carbon and nitrogen limitations for leaf-on events
  *
- * This function is called where leaf-on flux is calculated
+ * This function is called where leaf-on events are processed
  *
- * @param[inout] leafOnFlux pointer to leaf-on flux variable, which may be
- *               modified if limitation is in effect
+ * @param[inout] leafOnCDemand pointer to leaf-on demand variable, which may be
+ *               modified if limitation is in effect (g C m^-2)
  */
-void checkLeafOnLimitation(double *leafOnFlux);
+void checkLeafOnLimitation(double *leafOnCDemand);
 
 /** Check for carbon-specific limitations
  *

src/sipnet/nitrogen.c

@@ -69,7 +69,7 @@
   // The litter org N flux is determined by the carbon fluxes from wood and leaf
   // litter (modified by leaf N resorption), and N loss due to mineralization.
   // N added via fertilization is handled elsewhere.
-  fluxes.nOrgLitter = getLeafLitterFlux() / params.leafCN -
+  fluxes.nOrgLitter = fluxes.leafLitter / params.leafCN -
                       fluxes.leafOffNResorption +
                       fluxes.woodLitter / params.woodCN - litterMin -
                       litterNToSoil + (soilNInputs * saturationFraction);
@@ -129,7 +129,7 @@
 
 // see nitrogen.h
 double calcUnclaimedStorageN(void) {
   double leafOnCFlux = fluxes.leafOnCreation + fluxes.eventLeafOnCreation;
   double leafOnNFlux = calcLeafOnNFromC(leafOnCFlux);
   double unclaimedStorage = envi.plantStorageN - leafOnNFlux * climate->length;
   // The fmax here should be unnecessary, as the leaf-on demand has been capped
@@ -188,11 +188,9 @@
          fluxes.fineRootCreation / params.fineRootCN);
   }
 
-  // Leaf litter resorption; at this point. Note that event leaf off is handled
-  // in events.c
-  double nResorp =
-      params.leafNResorptionFrac * getLeafLitterFlux() / params.leafCN;
-  fluxes.leafOffNResorption = nResorp;
+  // Leaf litter resorption
+  fluxes.leafOffNResorption =
+      params.leafNResorptionFrac * fluxes.leafLitter / params.leafCN;
 
   // TODO: [consider] Should we resorb N from wood litter?
 }
@@ -215,17 +213,14 @@
   // TBD: add equation numbers once published
 
   // Storage N changes
-  // fluxes nUptake  and nFixation handle part of the demand flux (see
+  // fluxes nUptake and nFixation handle part of the demand flux (see
   // calcNFixationAndUptakeFluxes), but we expect the rest to come from the
   // storage pool
   double nDemandFlux = calcPlantNDemandFlux();
   double storageDemandFlux = nDemandFlux - fluxes.nUptake - fluxes.nFixation;
   // leaf-on; fluxes.eventLeafOnCreation is handled in events.c
-  double leafOnNFlux = calcLeafOnNFromC(fluxes.leafOnCreation);
   envi.plantStorageN +=
-      (fluxes.leafOffNResorption + fluxes.reductionNResorption -
-       storageDemandFlux - leafOnNFlux) *
-      climate->length;
+      (fluxes.reductionNResorption - storageDemandFlux) * climate->length;
 
   // Unmet uptake plus other fluxes go to soil mineral N (note we have one
   // mineral pool for soil+litter).
@@ -240,12 +235,13 @@
   envi.litterN += fluxes.nOrgLitter * climate->length;
 }
 
-void calcLeafOffNEffects(double leafOffC, double *resorptionFlux,
-                         double *litterFlux) {
+void calcLeafOffNEffects(double leafOffC, double *resorptionAmount,
+                         double *litterAmount) {
+  // TODO: move to processPendingEvents?
   double climLen = climate->length;
   double leafN = leafOffC / params.leafCN;
   double leafNResorption = leafN * params.leafNResorptionFrac;
   double litterNAdd = leafN - leafNResorption;
-  *resorptionFlux += leafNResorption / climLen;
-  *litterFlux += litterNAdd / climLen;
+  *resorptionAmount += leafNResorption / climLen;
+  *litterAmount += litterNAdd / climLen;
 }

src/sipnet/nitrogen.h

@@ -81,10 +81,9 @@ void calcNResorptionFluxes(void);
 /**
  * Helper function for handling leaf-off events
  *
- * This is used by both processEventsForNitrogen and the nitrogen limitation
- * check in limitations.
+ * This is used by processPendingEvents.
  */
-void calcLeafOffNEffects(double leafOffC, double *resorptionFlux,
-                         double *litterFlux);
+void calcLeafOffNEffects(double leafOffC, double *resorptionAmount,
+                         double *litterAmount);
 
 #endif  // NITROGEN_H

src/sipnet/restart.c

@@ -34,7 +34,7 @@
 #define NUM_SCHEMA_FIELDS 5
 #define NUM_MEAN_META_FIELDS 5
 #define NUM_ENVI_FIELDS 13
-#define NUM_TRACKER_FIELDS 33
+#define NUM_TRACKER_FIELDS 32
 #define NUM_PHENOLOGY_TRACKERS_FIELDS 3
 #define NUM_SURVIVAL_TRACKERS_FIELDS 1
 #define NUM_EVENT_TRACKERS_FIELDS 8
@@ -65,7 +65,7 @@
                    RESTART_SCHEMA_LAYOUT_SURVIVAL_TRACKERS_SIZE,
                "Restart schema drift: PlantSurvival changed; update "
                "schema_layout.* checks");
 _Static_assert(sizeof(EventTrackers) ==
                    RESTART_SCHEMA_LAYOUT_EVENT_TRACKERS_SIZE,
                "Restart schema drift: EventTrackers changed; update "
                "schema_layout.* checks");
@@ -252,7 +252,6 @@
   state->trackersPF[ind++] = (StateField){"trackers.yearlyRh",           FT_DOUBLE, &trackers.yearlyRh,           0};
   state->trackersPF[ind++] = (StateField){"trackers.yearlyNpp",          FT_DOUBLE, &trackers.yearlyNpp,          0};
   state->trackersPF[ind++] = (StateField){"trackers.yearlyNee",          FT_DOUBLE, &trackers.yearlyNee,          0};
-  state->trackersPF[ind++] = (StateField){"trackers.yearlyLitter",       FT_DOUBLE, &trackers.yearlyLitter,       0};
   state->trackersPF[ind++] = (StateField){"trackers.totGpp",             FT_DOUBLE, &trackers.totGpp,             0};
   state->trackersPF[ind++] = (StateField){"trackers.totRtot",            FT_DOUBLE, &trackers.totRtot,            0};
   state->trackersPF[ind++] = (StateField){"trackers.totRa",              FT_DOUBLE, &trackers.totRa,              0};

src/sipnet/sipnet.c

@@ -795,7 +795,7 @@
 void calcLeafOnOffFluxes(void) {
   // Calc additional fluxes at start/end of growing season
   // Note that these are basically events, and we will track them as such
   double len = climate->length;
 
   // first check for new year; if new year, reset trackers (since we haven't
   // done leaf growth or fall yet in this new year):
@@ -808,23 +808,16 @@
   // check for start of growing season:
   if (!phenologyTrackers.didLeafGrowth && pastLeafGrowth()) {
     // we just reached the start of the growing season
-    double leafOn = params.leafGrowth / len;
-    double leafOnFromWood = 0;
-    checkLeafOnLimitation(&leafOn);
-    fluxes.leafOnCreation += leafOn;
-    double totalSourceC = envi.plantWoodC + envi.coarseRootC;
-    if (totalSourceC > TINY) {
-      leafOnFromWood = leafOn * envi.plantWoodC / totalSourceC;
-      fluxes.leafOnCreationFromWood += leafOnFromWood;
-    }
+    // register this for handling as an event
+    eventTrackers.totalLeafOnAmount += params.leafGrowth;
     phenologyTrackers.didLeafGrowth = 1;
   }
 
   // check for end of growing season:
   if (!phenologyTrackers.didLeafFall && pastLeafFall()) {
     // we just reached the end of the growing season
-    double leafOff = (envi.plantLeafC * params.fracLeafFall) / len;
-    fluxes.leafOffLitter += leafOff;
+    // register this for handling as an event
+    eventTrackers.totalLeafOffFraction += params.fracLeafFall;
     phenologyTrackers.didLeafFall = 1;
   }
 }
@@ -1441,10 +1434,6 @@
   trackers.soilWetnessFrac =
       (oldSoilWater + envi.soilWater) / (2.0 * params.soilWHC);
 
-  // If we get another event flux in this function, we should create an
-  // updateTrackersForEvents() function in events.c|h
-  trackers.yearlyLitter += getLeafLitterFlux() + fluxes.eventLeafOffLitterC;
-
   if (ctx.gdd) {
     trackers.gdd += climate->gdd;
   } else {
@@ -1572,17 +1561,15 @@
                           fluxes.fineRootCreation + fluxes.coarseRootCreation;
   envi.plantCAccountingDelta +=
       ((fluxes.photosynthesis - r_a) - nppAllocations) * climate->length;
-  envi.plantWoodC += (fluxes.woodCreation - fluxes.woodLitter -
-                      fluxes.leafOnCreationFromWood) *
-                     climate->length;
+  envi.plantWoodC +=
+      (fluxes.woodCreation - fluxes.woodLitter) * climate->length;
 
   // :: from [1], eq (A2), where:
   //     L   = fluxes.leafCreation
   //     L_L = fluxes.leafLitter
   // Note: we have split leafCreation into two parts
   envi.plantLeafC +=
-      (fluxes.leafCreation + fluxes.leafOnCreation - getLeafLitterFlux()) *
-      climate->length;
+      (fluxes.leafCreation + -fluxes.leafLitter) * climate->length;
 
   // :: from [1], eq (A4), where:
   //     P = (fluxes.rain - fluxes.immedEvap)
@@ -1622,8 +1609,8 @@
         ctx.carbonSaturation ? unitClip(envi.soilC / params.soilCSaturation)
                              : 0.0;
     // :: from [2], litter model description
-    envi.litterC += (fluxes.woodLitter + getLeafLitterFlux() +
-                     (soilInputs * saturationFraction) - fluxes.litterToSoil -
+    envi.litterC += (fluxes.woodLitter + fluxes.leafLitter +
+                     soilInputs * saturationFraction - fluxes.litterToSoil -
                      fluxes.rLitter - fluxes.litterMethane) *
                     climate->length;
     // from [2] and [3], litter and root terms respectively
@@ -1639,18 +1626,15 @@
     // :: from [3], root terms
     envi.soilC +=
         (fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.woodLitter +
-         getLeafLitterFlux() - fluxes.rSoil - fluxes.soilMethane) *
+         fluxes.leafLitter - fluxes.rSoil - fluxes.soilMethane) *
         climate->length;
   }
 
   // :: from [3], root model description, except that we deal with
   //    respiration in the woodC calculations (or, put another way, it's
   //    implicit in the creation term)
-  double leafOnCreationFromRoot =
-      fluxes.leafOnCreation - fluxes.leafOnCreationFromWood;
-  envi.coarseRootC += (fluxes.coarseRootCreation - fluxes.coarseRootLoss -
-                       leafOnCreationFromRoot) *
-                      climate->length;
+  envi.coarseRootC +=
+      (fluxes.coarseRootCreation - fluxes.coarseRootLoss) * climate->length;
   envi.fineRootC +=
       (fluxes.fineRootCreation - fluxes.fineRootLoss) * climate->length;
 }
@@ -1830,7 +1814,13 @@
   updatePoolsAndBalance();
 
   ///////////////////////
-  // 3. Update trackers
+  // 3. Process Pending Events
+
+  // Harvest and leaf on/off events have their own handling
+  processPendingEvents();
+
+  ///////////////////////
+  // 4. Update trackers
 
   updateTrackers(oldSoilWater);
 

src/sipnet/state.c

@@ -10,14 +10,11 @@ ClimateNode *climate;  // current climate
 Params params;
 Envi envi;  // state variables
 Fluxes fluxes;  // fluxes
+PendingEvents pendingEvents;
 Trackers trackers;
 PhenologyTrackers phenologyTrackers;
 PlantSurvivalTracker plantSurvivalTracker;
 
 double getTotalWoodC(void) {
   return envi.plantWoodC + envi.plantCAccountingDelta;
 }
-
-double getLeafLitterFlux(void) {
-  return fluxes.leafLitter + fluxes.leafOffLitter;
-}