diff --git a/src/sipnet/limitations.c b/src/sipnet/limitations.c index 918c109d..2971f755 100644 --- a/src/sipnet/limitations.c +++ b/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(); diff --git a/src/sipnet/nitrogen.c b/src/sipnet/nitrogen.c index 6a4f2a96..b50d3985 100644 --- a/src/sipnet/nitrogen.c +++ b/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 diff --git a/src/sipnet/nitrogen.h b/src/sipnet/nitrogen.h index 4d549074..aff91a6f 100644 --- a/src/sipnet/nitrogen.h +++ b/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 *