diff --git a/src/sipnet/debug_log.c b/src/sipnet/debug_log.c index 8253e0eb..501739c7 100644 --- a/src/sipnet/debug_log.c +++ b/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 @@ -185,7 +185,6 @@ void initDebugArrays() { 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}, diff --git a/src/sipnet/events.c b/src/sipnet/events.c index 3d05db57..ce56e7fb 100644 --- a/src/sipnet/events.c +++ b/src/sipnet/events.c @@ -544,19 +544,20 @@ void writeEventsOut(void) { // Last, check for calculated leaf on/off "events" double climLen = climate->length; - if (fluxes.leafOnCreation > 0.0) { + if (!hasSeenLeafOn && eventTrackers.totalLeafOnAmount > 0) { // clang-format off writeComputedEventOut( climate->year, climate->day, eventTypeToString(LEAFON), 2, - "leafOnCreation", fluxes.leafOnCreation * climLen, - "leafOnCreationFromWood", fluxes.leafOnCreationFromWood * climLen); + "leafOnCreation", pendingEvents.leafOnCreation, + "leafOnCreationFromWood", pendingEvents.leafOnCreationFromWood); // clang-format on } - if (fluxes.leafOffLitter > 0.0) { + if (!hasSeenLeafOff && eventTrackers.totalLeafOffFraction > 0) { // clang-format off + // TODO: this isn't right; what about N? writeComputedEventOut( climate->year, climate->day, eventTypeToString(LEAFOFF), 1, - "leafOffLitter", fluxes.leafOffLitter * climLen); + "leafOffLitter", pendingEvents.leafOffC * climLen); // clang-format on } } @@ -630,10 +631,10 @@ void processEventsForCarbon(EventNode *event) { // Reset harvest tracking eventTrackers.harvestTrackers = (HarvestTrackers){0}; + eventTrackers.totalLeafOnAmount = 0.0; + eventTrackers.totalLeafOffFraction = 0.0; - // Make sure we don't have more than 100% leaves fall for leaf-off - double totalLeafFallFrac = 0.0; - + // If we see more than one of these, log that we are combining them int hasSeenLeafOn = 0; int hasSeenLeafOff = 0; int hasSeenHarvest = 0; @@ -730,6 +731,7 @@ void processEventsForCarbon(EventNode *event) { double totalMass = aboveMass + belowMass; HarvestTrackers *ht = &eventTrackers.harvestTrackers; if (totalMass > TINY) { + ht->harvestOccurred = 1; double massRemoved = fracRA * aboveMass + fracRB * belowMass; double massTransferred = fracTA * aboveMass + fracTB * belowMass; ht->totalFracRemoved += massRemoved / totalMass; @@ -756,40 +758,8 @@ void processEventsForCarbon(EventNode *event) { event->year, event->day); } - // Litter increase - 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); - const double woodDelta = -envi.plantWoodC * (fracRA + fracTA); - const double accountingDelta = - -envi.plantCAccountingDelta * (fracRA + fracTA); - // Below-ground changes: - const double fineDelta = -envi.fineRootC * (fracRB + fracTB); - const double coarseDelta = -envi.coarseRootC * (fracRB + fracTB); - - // Pool updates: - if (!ctx.litterPool) { - // send it all to the soil - soilAdd += litterAdd; - litterAdd = 0.0; - } - fluxes.harvestLitterC += litterAdd / climLen; - fluxes.harvestSoilC += soilAdd / climLen; - fluxes.eventLeafC += leafDelta / climLen; - fluxes.eventWoodC += woodDelta / climLen; - fluxes.eventAccountingC += accountingDelta / climLen; - fluxes.eventFineRootC += fineDelta / climLen; - fluxes.eventCoarseRootC += coarseDelta / climLen; - - // MASS BALANCE: removed fractions are system outputs - const double outputC = ((woodC + envi.plantLeafC) * fracRA + - (envi.fineRootC + envi.coarseRootC) * fracRB); - fluxes.eventOutputC += outputC / climLen; - - // No log here, as eventLeafC may change due to limitation checks + // Harvest is processed as a pending event + } break; case TILLAGE: { // BIG NOTE: this is the one event type that is NOT modeled as a flux; @@ -831,20 +801,10 @@ void processEventsForCarbon(EventNode *event) { event->year, event->day); } hasSeenLeafOn = 1; + eventTrackers.totalLeafOnAmount += params.leafGrowth; - double leafOnFlux = params.leafGrowth / climLen; - checkLeafOnLimitation(&leafOnFlux); - fluxes.eventLeafOnCreation += leafOnFlux; - double totalSourceC = envi.plantWoodC + envi.coarseRootC; - if (totalSourceC > TINY) { - double leafOnFluxFromWood = - leafOnFlux * envi.plantWoodC / totalSourceC; - fluxes.eventLeafOnCreationFromWood += leafOnFluxFromWood; - } - // Unlike planting, this is NOT a system input, so no adjustments to - // eventInputC + // Leaf on is processed as a pending event - // No log here, in case we have more than one } break; case LEAFOFF: { if (hasSeenLeafOff) { @@ -854,26 +814,24 @@ void processEventsForCarbon(EventNode *event) { event->year, event->day); } hasSeenLeafOff = 1; - - double frac = params.fracLeafFall; - totalLeafFallFrac += frac; - if (totalLeafFallFrac > 1.0) { + eventTrackers.totalLeafOffFraction += params.fracLeafFall; + if (eventTrackers.totalLeafOffFraction > 1.0) { + // TODO: [consider] just trim to 100%? logError("Total leaf fall for leaf-off event(s) is greater than 100 " "percent (%.3f) for year %d day %d\n", - totalLeafFallFrac, event->year, event->day); + eventTrackers.totalLeafOffFraction, event->year, event->day); exit(EXIT_CODE_BAD_PARAMETER_VALUE); } - double leafOff = envi.plantLeafC * params.fracLeafFall; - fluxes.eventLeafOffLitterC += leafOff / climLen; - // No log here, as this value may change due to limitation checks + // Leaf off is processed as a pending event + } break; - case PLANTDEATH: + case PLANTDEATH: { // There should be no way to get here, but covering our bases... logWarning("PLANTDEATH event found for year %d day %d, but not " "implemented as an input event; ignoring\n", event->year, event->day); - break; + } break; default: logError("Unknown event type (%d) in processEvents()\n", event->type); exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM); @@ -900,11 +858,6 @@ void processEventsForNitrogen(EventNode *event) { // Checks performed in processEventsForCarbon are not repeated here unless // necessary for nitrogen handling. - // Make sure not to double-process these event types, in case there is more - // than one - int hasSeenLeafOff = 0; - int hasSeenHarvest = 0; - while (event != NULL && event->year <= climYear && event->day <= climDay) { switch (event->type) { case IRRIGATION: { @@ -928,55 +881,9 @@ void processEventsForNitrogen(EventNode *event) { appendLog(event, 1, "eventInputN", inputN); } break; - case HARVEST: { - if (hasSeenHarvest) { - break; // already processed the combined harvest events - } - hasSeenHarvest = 1; - - // Harvest can both remove biomass and move biomass to the soil/litter - // pools - HarvestTrackers *ht = &eventTrackers.harvestTrackers; - const double fracRA = ht->totalFracRemovedAbove; - const double fracRB = ht->totalFracRemovedBelow; - const double fracTA = ht->totalFracTransferredAbove; - const double fracTB = ht->totalFracTransferredBelow; - - // No need to allocate to biomass N pools, we don't track that N - // explicitly. We do need to handle soil and litter N, though. - // Note: ctx.nitrogenCycle implies ctx.litterPool - // Litter N increase - const double totalAbove = (envi.plantLeafC / params.leafCN) + - (envi.plantWoodC / params.woodCN); - const double totalBelow = (envi.fineRootC / params.fineRootCN) + - (envi.coarseRootC / params.woodCN); - double litterNAdd = fracTA * totalAbove; - double soilNAdd = fracTB * totalBelow; - fluxes.eventSoilOrgN += soilNAdd / climLen; - fluxes.harvestLitterN += litterNAdd / climLen; - - // MASS BALANCE: removed fractions are system outputs - // just plantWoodC here, not woodC - double outputN = (envi.plantWoodC / params.woodCN + - envi.plantLeafC / params.leafCN) * - fracRA + - (envi.fineRootC / params.fineRootCN + - envi.coarseRootC / params.woodCN) * - fracRB; - fluxes.eventOutputN += outputN / climLen; - - // NOTE: leaf/litter/output C may have been adjusted by - // limitation checks; make sure to do the same for N - double leafDiffFlux = ht->leafCAdjustment; - - // Figure out how to distribute - double litterFlux, outputFlux; - partitionLeafAdjustment(leafDiffFlux, &litterFlux, &outputFlux); - fluxes.harvestLitterN -= litterFlux / params.leafCN; - fluxes.eventOutputN -= outputFlux / params.leafCN; - } break; + case HARVEST: case TILLAGE: { - // Nothing to do for tillage events + // Nothing to do for harvest (pending event) or tillage } break; case FERTILIZATION: { const FertilizationParams *fertParams = event->eventParams; @@ -998,30 +905,9 @@ void processEventsForNitrogen(EventNode *event) { "eventInputN", (orgN + minN)); // clang-format on } break; - case LEAFON: { - // Nitrogen is handled implicitly by relative CN ratios. Missing N - // from low-N wood to higher-N leaves is accounted for in - // calcNFixationAndUptakeFluxes() via calcPlantNDemandFlux() - - // Unlike planting, this is NOT a system input, so no adjustments to - // eventInputN - } break; + case LEAFON: case LEAFOFF: { - if (hasSeenLeafOff) { - break; // already processed the combined leaf-off events - } - hasSeenLeafOff = 1; - - // We need to use fluxes.eventLeafOffLitter here instead of - // recalculating it, as it may have been reduced. - - // Nitrogen - need to account for leaf N moving to litter, as with - // harvests - double leafOff = fluxes.eventLeafOffLitterC * climLen; - calcLeafOffNEffects(leafOff, &fluxes.eventLeafOffNResorption, - &fluxes.eventLeafOffLitterN); - - // No log here, as these values may change + // Leaf on/off are processed as a pending event } break; case PLANTDEATH: // Nothing to do here @@ -1035,55 +921,185 @@ void processEventsForNitrogen(EventNode *event) { } } -void updatePoolsForEvents(void) { - // CARBON - // Harvest and planting events - envi.plantWoodC += fluxes.eventWoodC * climate->length; - envi.plantCAccountingDelta += fluxes.eventAccountingC * climate->length; - envi.plantLeafC += fluxes.eventLeafC * climate->length; +void processPendingEvents(void) { + // Certain events have the capacity to drive pools to zero in one step, + // which causes significant complexity when trying to handle them as other + // events. So, instead of introducing a lot of corrective code to handle + // overshoots, we process them here after pools have been updated due to + // other modeling effects. This way, we can more easily guarantee that these + // events don't do more than they should. + // + // Order of processing: + // 1. leaf on + // 2. harvest + // 3. leaf off + + double len = climate->length; + + ///////////////// + // LEAF ON + if (eventTrackers.totalLeafOnAmount > 0.0) { + double leafOnAmount = eventTrackers.totalLeafOnAmount; + checkLeafOnLimitation(&leafOnAmount); + pendingEvents.leafOnCreation = leafOnAmount; + + double totalSourceC = envi.plantWoodC + envi.coarseRootC; + if (totalSourceC > TINY) { + double leafOnCreationFromWood = + leafOnAmount * envi.plantWoodC / totalSourceC; + double leafOnCreationFromRoot = + leafOnAmount * envi.coarseRootC / totalSourceC; + pendingEvents.leafOnCreation = + leafOnCreationFromWood + leafOnCreationFromRoot; + pendingEvents.leafOnCreationFromWood = leafOnCreationFromWood; + // Update pools + envi.plantLeafC += leafOnAmount; + envi.plantWoodC -= leafOnCreationFromWood; + envi.coarseRootC -= leafOnCreationFromRoot; + // Unlike planting, this is NOT a system input, so no adjustments to + // eventInputC or eventInputN + // Nitrogen is handled implicitly by relative CN ratios. Missing N + // from low-N wood to higher-N leaves is accounted for in + // calcNFixationAndUptakeFluxes() via calcPlantNDemandFlux() + } + } - // Harvest and fertilization events - envi.soilC += (fluxes.eventSoilC + fluxes.harvestSoilC) * climate->length; - if (ctx.litterPool) { - envi.litterC += - (fluxes.eventLitterC + fluxes.harvestLitterC) * climate->length; + ///////////////// + // HARVEST + HarvestTrackers *ht = &eventTrackers.harvestTrackers; + const double fracRA = ht->totalFracRemovedAbove; + const double fracRB = ht->totalFracRemovedBelow; + const double fracTA = ht->totalFracTransferredAbove; + const double fracTB = ht->totalFracTransferredBelow; + const double woodC = envi.plantWoodC + envi.plantCAccountingDelta; + + // Litter increase + 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); + const double woodDelta = -envi.plantWoodC * (fracRA + fracTA); + const double accountingDelta = + -envi.plantCAccountingDelta * (fracRA + fracTA); + // Below-ground changes: + const double fineDelta = -envi.fineRootC * (fracRB + fracTB); + const double coarseDelta = -envi.coarseRootC * (fracRB + fracTB); + + if (!ctx.litterPool) { + // send it all to the soil + soilAdd += litterAdd; + litterAdd = 0.0; } + pendingEvents.harvestLitterC += litterAdd; + pendingEvents.harvestSoilC += soilAdd; + pendingEvents.harvestLeafC += leafDelta; + pendingEvents.harvestWoodC += woodDelta; + pendingEvents.harvestAccountingC += accountingDelta; + pendingEvents.harvestFineRootC += fineDelta; + pendingEvents.harvestCoarseRootC += coarseDelta; + + // MASS BALANCE: removed fractions are system outputs + const double outputC = ((woodC + envi.plantLeafC) * fracRA + + (envi.fineRootC + envi.coarseRootC) * fracRB); + pendingEvents.harvestOutputC += outputC; + + // Update pools + envi.plantWoodC += woodDelta; + envi.plantCAccountingDelta += accountingDelta; + envi.plantLeafC += leafDelta; + envi.fineRootC += fineDelta; + envi.coarseRootC += coarseDelta; + envi.litterC += litterAdd; + envi.soilC += soilAdd; + // The balance checker is expecting this one as a flux + fluxes.eventOutputC += outputC / len; - // Leaf on and off events - // Leaf on draws from wood and coarse root pools in proportion to their sizes - envi.plantWoodC -= fluxes.eventLeafOnCreationFromWood * climate->length; - double eventLeafOnCreationFromRoot = - fluxes.eventLeafOnCreation - fluxes.eventLeafOnCreationFromWood; - envi.coarseRootC -= eventLeafOnCreationFromRoot * climate->length; - envi.plantLeafC += (fluxes.eventLeafOnCreation - fluxes.eventLeafOffLitterC) * - climate->length; - if (ctx.litterPool) { - envi.litterC += fluxes.eventLeafOffLitterC * climate->length; - } else { - envi.soilC += fluxes.eventLeafOffLitterC * climate->length; + if (ctx.nitrogenCycle) { + // No need to allocate to biomass N pools, we don't track that N + // explicitly. We do need to handle soil and litter N, though. + // Note: ctx.nitrogenCycle implies ctx.litterPool + // Litter/soil N increase + const double totalAbove = + (envi.plantLeafC / params.leafCN) + (envi.plantWoodC / params.woodCN); + const double totalBelow = (envi.fineRootC / params.fineRootCN) + + (envi.coarseRootC / params.woodCN); + double litterNAdd = fracTA * totalAbove; + double soilNAdd = fracTB * totalBelow; + pendingEvents.harvestSoilN += soilNAdd; + pendingEvents.harvestLitterN += litterNAdd; + + // MASS BALANCE: removed fractions are system outputs + // just plantWoodC here, not woodC + double outputN = + (envi.plantWoodC / params.woodCN + envi.plantLeafC / params.leafCN) * + fracRA + + (envi.fineRootC / params.fineRootCN + + envi.coarseRootC / params.woodCN) * + fracRB; + pendingEvents.harvestOutputN += outputN; + // Update pools + envi.litterN += litterNAdd; + envi.soilOrgN += soilNAdd; } - // Harvest and planting events + ///////////////// + // LEAF OFF + if (eventTrackers.totalLeafOffFraction > 0.0) { + double leafOffFraction = eventTrackers.totalLeafOffFraction; + double leafOff = envi.plantLeafC * leafOffFraction; + + // Update pools + envi.plantLeafC -= leafOff; + pendingEvents.leafOffC -= leafOff; + if (ctx.litterPool) { + envi.litterC += leafOff; + } else { + envi.soilC += leafOff; + } + + if (ctx.nitrogenCycle) { + // Nitrogen - need to account for leaf N moving to litter, as with + // harvests + calcLeafOffNEffects(leafOff, &pendingEvents.leafOffNResorption, + &pendingEvents.leafOffLitterN); + envi.plantStorageN += pendingEvents.leafOffNResorption; + envi.litterN += pendingEvents.leafOffLitterN; + } + } +} + +void updatePoolsForEvents(void) { + double len = climate->length; + + // Reminder: harvest, leaf on, and leaf off events are handled in + // processPendingEvents + + // CARBON + // Planting events + envi.plantWoodC += fluxes.eventWoodC * len; + envi.plantLeafC += fluxes.eventLeafC * len; envi.coarseRootC += fluxes.eventCoarseRootC * climate->length; envi.fineRootC += fluxes.eventFineRootC * climate->length; + // Fertilization events + envi.soilC += fluxes.eventSoilC * len; + if (ctx.litterPool) { + envi.litterC += fluxes.eventLitterC * len; + } + // WATER // Irrigation events envi.soilWater += fluxes.eventSoilWater * climate->length; // NITROGEN - // Harvest, fertilization, and leaf-off events + // 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 + fluxes.harvestLitterN + - fluxes.eventLeafOffLitterN) * - climate->length; - double leafOnNFlux = calcLeafOnNFromC(fluxes.eventLeafOnCreation); - envi.plantStorageN += - (fluxes.eventLeafOffNResorption - leafOnNFlux) * climate->length; + envi.minN += fluxes.eventMinN * len; + envi.litterN += fluxes.eventLitterN * len; } } diff --git a/src/sipnet/events.h b/src/sipnet/events.h index 88e96984..c3096e4e 100644 --- a/src/sipnet/events.h +++ b/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; diff --git a/src/sipnet/limitations.c b/src/sipnet/limitations.c index 5feb058d..713095f9 100644 --- a/src/sipnet/limitations.c +++ b/src/sipnet/limitations.c @@ -11,19 +11,13 @@ #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; } @@ -31,7 +25,7 @@ void checkLeafOnLimitation(double *leafOnFlux) { // 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); } diff --git a/src/sipnet/limitations.h b/src/sipnet/limitations.h index 0e88f81d..4678c889 100644 --- a/src/sipnet/limitations.h +++ b/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 * diff --git a/src/sipnet/nitrogen.c b/src/sipnet/nitrogen.c index a60fa3fb..ec2e6433 100644 --- a/src/sipnet/nitrogen.c +++ b/src/sipnet/nitrogen.c @@ -69,7 +69,7 @@ static void calcNPoolFluxes(void) { // 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); @@ -188,11 +188,9 @@ void calcNResorptionFluxes(void) { 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 @@ void updateNitrogenPools(void) { // 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 @@ void updateNitrogenPools(void) { 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; } diff --git a/src/sipnet/nitrogen.h b/src/sipnet/nitrogen.h index 8a30830a..5cbbeb00 100644 --- a/src/sipnet/nitrogen.h +++ b/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 diff --git a/src/sipnet/restart.c b/src/sipnet/restart.c index 7541a991..f61ff7b4 100644 --- a/src/sipnet/restart.c +++ b/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 @@ -252,7 +252,6 @@ void initResetState(RestartState *state, MeanTracker *npp) { 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}; diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 22e50510..5876b79b 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -808,23 +808,16 @@ void calcLeafOnOffFluxes(void) { // 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 @@ void updateTrackers(double oldSoilWater) { 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 @@ void updateMainPools(void) { 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 @@ void updatePoolsForSoil(void) { 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 @@ void updatePoolsForSoil(void) { // :: 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 @@ void updateState(void) { updatePoolsAndBalance(); /////////////////////// - // 3. Update trackers + // 3. Process Pending Events + + // Harvest and leaf on/off events have their own handling + processPendingEvents(); + + /////////////////////// + // 4. Update trackers updateTrackers(oldSoilWater); diff --git a/src/sipnet/state.c b/src/sipnet/state.c index b082cbcd..ca7db50b 100644 --- a/src/sipnet/state.c +++ b/src/sipnet/state.c @@ -10,6 +10,7 @@ ClimateNode *climate; // current climate Params params; Envi envi; // state variables Fluxes fluxes; // fluxes +PendingEvents pendingEvents; Trackers trackers; PhenologyTrackers phenologyTrackers; PlantSurvivalTracker plantSurvivalTracker; @@ -17,7 +18,3 @@ PlantSurvivalTracker plantSurvivalTracker; double getTotalWoodC(void) { return envi.plantWoodC + envi.plantCAccountingDelta; } - -double getLeafLitterFlux(void) { - return fluxes.leafLitter + fluxes.leafOffLitter; -} diff --git a/src/sipnet/state.h b/src/sipnet/state.h index d57eb92b..7c1d3f72 100644 --- a/src/sipnet/state.h +++ b/src/sipnet/state.h @@ -546,10 +546,10 @@ typedef struct FluxVars { double woodCreation; // creation from leaf-on event; separated from leafCreation for N balance (as // we take this from the wood pool) (g C * m^-2 ground area * day^-1) - double leafOnCreation; + // double leafOnCreation; TODO: remove // Portion of leaf-on creation C that comes from wood C (the rest comes from // coarse root C) - double leafOnCreationFromWood; + // double leafOnCreationFromWood; TODO: remove // leaf fall from calculated leaf off 'events' (gdd, soil temp, day of year) // (g C * m^-2 grount area * day^-1) double leafOffLitter; @@ -663,6 +663,33 @@ typedef struct FluxVars { // Global var extern Fluxes fluxes; +// Mostly for debugging, though these do make logging the pending events easier +typedef struct PendingEventVars { + // leaf on + double leafOnCreation; + double leafOnCreationFromWood; + + // harvest + double harvestLeafC; + double harvestWoodC; + double harvestAccountingC; + double harvestFineRootC; + double harvestCoarseRootC; + double harvestLitterC; + double harvestSoilC; + double harvestLitterN; + double harvestSoilN; + double harvestOutputC; + double harvestOutputN; + + // leaf off + double leafOffC; + double leafOffNResorption; + double leafOffLitterN; +} PendingEvents; + +extern PendingEvents pendingEvents; + typedef struct TrackerVars { // variables to track various things // g C * m^-2 taken up in this time interval; GROSS photosynthesis double gpp; @@ -707,8 +734,6 @@ typedef struct TrackerVars { // variables to track various things double yearlyNpp; // g C * m^-2 given off, year to date double yearlyNee; - // g C * m^-2 litterfall, year to date: SUM litter - double yearlyLitter; // g C * m^-2 taken up, to date: GROSS photosynthesis double totGpp; // g C * m^-2 respired, to date @@ -771,6 +796,4 @@ extern PlantSurvivalTracker plantSurvivalTracker; double getTotalWoodC(void); -double getLeafLitterFlux(void); - #endif // SIPNET_STATE_H