SIP336 Add C and N limits to leaf on

src/sipnet/limitations.c

@@ -9,6 +9,73 @@
 #include "nitrogen.h"
 #include "state.h"
 
+/**
+ * Check for carbon and nitrogen limitations for leaf-on events
+ */
+static void checkLeafOnLimitation(void) {
+  // Leaf on events are limited by:
+  // * leafGrowth parameter (this is already calculated)
+  // * available carbon
+  // * available nitrogen
+  double leafOnFlux = fluxes.leafOnCreation + fluxes.eventLeafOnCreation;
+  double leafOnCDemand = leafOnFlux * climate->length;
+  double leafOnNDemand = 0.0;
+  double nLimitation = 1.0;
+  double availableN = 0.0;
+
+  if (leafOnCDemand < TINY) {
+    // Nothing to check
+    return;
+  }
+
+  // First up, carbon. We do not draw from the storage pool for this.
+  double availableC = envi.plantWoodC + envi.coarseRootC;
+  double cLimitation = (availableC > TINY) ? (availableC / leafOnCDemand) : 0;
+
+  // Next, nitrogen. This one is a little trickier, as 'available nitrogen' is
+  // harder to calculate (but we defer that to the nitrogen module).
+  if (ctx.nitrogenCycle) {
+    // Needed N for this transfer is (what leaves need) - (what wood provides)
+    // Reminder: both wood and coarseRoot use params.woodCN
+    leafOnNDemand =
+        leafOnCDemand / params.leafCN - leafOnCDemand / params.woodCN;
+    double minNFluxes = calcMinNNonUptakeFluxes() + fluxes.eventMinN;
+    availableN = envi.minN + minNFluxes * climate->length;
+    nLimitation = (availableN > TINY) ? (availableN / leafOnNDemand) : 0;
+  }
+  double limitation = fmin(cLimitation, nLimitation);
+  limitation = fmax(fmin(limitation, 1.0), 0.0);
+
+  if (limitation < 1) {
+    fluxes.leafOnCreation *= limitation;
+    fluxes.eventLeafOnCreation *= limitation;
+    if (ctx.nitrogenCycle) {
+      // We need to recalculate fixation and uptake since the demand has now
+      // changed
+      calcNFixationAndUptakeFluxes();
+
+      logInfo("Leaf on creation %.4f C / %.4f N exceeds available C %.4f / N "
+              "%.4f, "
+              "reducing leaf on by %.2f%% on year %d day %d time %.3f\n",
+              leafOnCDemand, leafOnNDemand, availableC, availableN,
+              (1 - limitation) * 100, climate->year, climate->day,
+              climate->time);
+      logInfo("envi.minN %.3f fluxes.nMin %.3f fluxes.nVol %.3f fluxes.nLeach "
+              "%.3f fluxes.nFix %.3f fluxes.nUptake %.3f\n",
+              envi.minN, fluxes.nMin * climate->length,
+              fluxes.nVolatilization * climate->length,
+              fluxes.nLeaching * climate->length,
+              fluxes.nFixation * climate->length,
+              fluxes.nUptake * climate->length);
+    } else {
+      logInfo("Leaf on creation %.4f exceeds available C %.4f, "
+              "reducing leaf on by %.2f%% on year %d day %d time %.3f\n",
+              leafOnCDemand, availableC, (1 - limitation) * 100, climate->year,
+              climate->day, climate->time);
+    }
+  }
+}
+
 /**
  * Check for nitrogen limitation, and reduce growth if needed
  */
@@ -46,7 +113,10 @@ static void checkNitrogenLimitation(void) {
 
 // See limitations.h
 void checkLimitations(void) {
-  // EXPECTED: calcLeafOnLimitation()
+
+  // First up - leaf on. This should be before the nitrogen limitation check,
+  // as the two are not independent.
+  checkLeafOnLimitation();
 
   if (ctx.nitrogenCycle) {
     checkNitrogenLimitation();

src/sipnet/nitrogen.c

@@ -39,18 +39,6 @@ static void calcNLeachingFlux(void) {
   fluxes.nLeaching = envi.minN * phi * params.nLeachingFrac;
 }
 
-/*!
- * Calculate plant N fixation and uptake fluxes
- */
-static void calcNFixationAndUptakeFluxes(void) {
-  // These values may change later if we are under nitrogen limitation
-  double nDemandFlux = calcPlantNDemand();
-  double nFixationFrac = calcNFixationFrac();
-
-  fluxes.nFixation = nFixationFrac * nDemandFlux;
-  fluxes.nUptake = (1 - nFixationFrac) * nDemandFlux;
-}
-
 /**
  * Calculate nitrogen fluxes for soil and litter pools
  */
@@ -84,6 +72,16 @@ static void calcNPoolFluxes(void) {
   fluxes.nMin = litterMin + soilMin;
 }
 
+// see nitrogen.h
+double calcAvailableNitrogen(void) {
+  // Return "available nitrogen", which is mineral N +/- relevant fluxes
+  // from this time step.
+  double availableN = envi.minN;
+  double minNDelta = calcMinNNonUptakeFluxes() * climate->length;
+
+  return availableN + minNDelta;
+}
+
 // see nitrogen.h
 double calcPlantNDemand(void) {
   if (!ctx.nitrogenCycle) {
@@ -130,10 +128,14 @@ double calcNFixationFrac(void) {
   return params.nFixationFracMax * nFixationInhibition;
 }
 
-// see nitrogen.h
-double calcAvailableNitrogen(void) {
-  // NOT USED YET
-  return 0.0;
+// See nitrogen.h
+void calcNFixationAndUptakeFluxes(void) {
+  // These values may change later if we are under nitrogen limitation
+  double nDemandFlux = calcPlantNDemand();
+  double nFixationFrac = calcNFixationFrac();
+
+  fluxes.nFixation = nFixationFrac * nDemandFlux;
+  fluxes.nUptake = (1 - nFixationFrac) * nDemandFlux;
 }
 
 // see nitrogen.h

src/sipnet/nitrogen.h

@@ -3,12 +3,15 @@
 
 // Nitrogen cycle related functions
 
-// TBD: will be used for leaf-on limit check; might be folded into
-// nitrogen use too, if appropriate
 /*!
  * Calculate available nitrogen for this step
  *
- * Takes into account demand, fixation, and...
+ * Calculate available nitrogen as mineral N pool plus relevant fluxes:
+ * - mineralization
+ * - leaching
+ * - volatilization
+ *
+ * This function should only be called AFTER fluxes are calculated
  *
  * @return available nitrogen
  */
@@ -37,6 +40,13 @@ double calcMinNNonUptakeFluxes(void);
  */
 double calcNFixationFrac(void);
 
+/*!
+ * Calculate plant N fixation and uptake fluxes.
+ *
+ * This function is called directly by the leaf-on limitation check
+ */
+void calcNFixationAndUptakeFluxes(void);
+
 /*!
  * Calculate all nitrogen fluxes
  *