From f37cd77f3ea58cbb9d42f11f2e66eb6d3268b204 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Wed, 16 Sep 2026 10:07:27 -0400 Subject: [PATCH 01/14] Interim commit --- src/sipnet/debug_log.c | 3 +- src/sipnet/events.c | 9 ++++-- src/sipnet/sipnet.c | 8 ++++- src/sipnet/state.h | 2 ++ .../events_output_header.out | 4 +-- .../events_output_no_header.out | 4 +-- .../test_events_types/testEventHarvest.c | 32 ++++++++++++++++++- 7 files changed, 53 insertions(+), 9 deletions(-) diff --git a/src/sipnet/debug_log.c b/src/sipnet/debug_log.c index fa8e99ba..483ff3fe 100644 --- a/src/sipnet/debug_log.c +++ b/src/sipnet/debug_log.c @@ -21,7 +21,7 @@ typedef struct DebugField { } DebugField; #define NUM_LOGGED_ENVI_FIELDS 13 -#define NUM_LOGGED_FLUX_FIELDS 56 +#define NUM_LOGGED_FLUX_FIELDS 57 #define NUM_LOGGED_TRACKER_FIELDS 33 #define NUM_LOGGED_PHEN_TRACKER_FIELDS 3 #define NUM_LOGGED_SURVIVAL_FIELDS 1 @@ -101,6 +101,7 @@ void initDebugArrays() { debugFields->fluxDF[ind++] = (DebugField){"reductionNResorption", DEBUG_FIELD_DOUBLE, &fluxes.reductionNResorption}, debugFields->fluxDF[ind++] = (DebugField){"eventLeafC", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafC}, debugFields->fluxDF[ind++] = (DebugField){"eventWoodC", DEBUG_FIELD_DOUBLE, &fluxes.eventWoodC}, + debugFields->fluxDF[ind++] = (DebugField){"eventAccountingC", DEBUG_FIELD_DOUBLE, &fluxes.eventAccountingC}, debugFields->fluxDF[ind++] = (DebugField){"eventFineRootC", DEBUG_FIELD_DOUBLE, &fluxes.eventFineRootC}, debugFields->fluxDF[ind++] = (DebugField){"eventCoarseRootC", DEBUG_FIELD_DOUBLE, &fluxes.eventCoarseRootC}, debugFields->fluxDF[ind++] = (DebugField){"eventEvap", DEBUG_FIELD_DOUBLE, &fluxes.eventEvap}, diff --git a/src/sipnet/events.c b/src/sipnet/events.c index cc99ce0b..3d2250f0 100644 --- a/src/sipnet/events.c +++ b/src/sipnet/events.c @@ -568,7 +568,9 @@ void processEvents(void) { // 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 = -woodC * (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); @@ -583,6 +585,7 @@ void processEvents(void) { fluxes.eventSoilC += soilAdd / climLen; fluxes.eventLeafC += leafDelta / climLen; fluxes.eventWoodC += woodDelta / climLen; + fluxes.eventAccountingC += accountingDelta / climLen; fluxes.eventFineRootC += fineDelta / climLen; fluxes.eventCoarseRootC += coarseDelta / climLen; @@ -620,11 +623,12 @@ void processEvents(void) { } // clang-format off writeEventOut( - gEvent, 10, + gEvent, 11, "eventSoilC", soilAdd, "eventLitterC", litterAdd, "eventLeafC", leafDelta, "eventWoodC", woodDelta, + "eventAccountingC", accountingDelta, "eventFineRootC", fineDelta, "eventCoarseRootC", coarseDelta, "eventSoilOrgN", soilNAdd, @@ -745,6 +749,7 @@ 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; // Harvest and fertilization events diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 5153c552..1de3f7cc 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -1529,8 +1529,14 @@ void initPhenologyTrackers(void) { // this year } -// Check that woodC and total root C are both positive +// Check that woodC and total root C are both positive and that there was no +// terminating harvest int hasSufficientBiomass(void) { + if (eventTrackers.harvestFracRemoved + eventTrackers.harvestFracTransferred > + 1.0 - TINY) { + return 0; + } + double totalWoodC = getTotalWoodC(); double totalRootC = envi.fineRootC + envi.coarseRootC; // We want to check that both plantWoodC AND totalWoodC are positive, as well diff --git a/src/sipnet/state.h b/src/sipnet/state.h index 658e106f..f77bc791 100644 --- a/src/sipnet/state.h +++ b/src/sipnet/state.h @@ -595,6 +595,8 @@ typedef struct FluxVars { double eventLeafC; // plantWoodC addition double eventWoodC; + // N-free accounting carbon removed by harvest + double eventAccountingC; // plantFineRootC addition double eventFineRootC; // plantCoarseRootC addition diff --git a/tests/sipnet/test_events_infrastructure/events_output_header.out b/tests/sipnet/test_events_infrastructure/events_output_header.out index ec34fc93..d00de444 100644 --- a/tests/sipnet/test_events_infrastructure/events_output_header.out +++ b/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 eventLeafC=3.00,eventWoodC=4.00,eventFineRootC=5.00,eventCoarseRootC=6.00,eventInputC=18.00,eventInputN=0.00 2023 70 irrig eventSoilWater=5.00,eventEvap=0.00 -2023 200 harv eventSoilC=1.90,eventLitterC=3.56,eventLeafC=-5.93,eventWoodC=-4.75,eventFineRootC=-3.73,eventCoarseRootC=-3.89,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=12.83,eventOutputN=0.00 +2023 200 harv eventSoilC=1.90,eventLitterC=3.56,eventLeafC=-5.93,eventWoodC=-4.75,eventAccountingC=-0.00,eventFineRootC=-3.73,eventCoarseRootC=-3.89,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=12.83,eventOutputN=0.00 2024 65 plant eventLeafC=3.00,eventWoodC=5.00,eventFineRootC=7.00,eventCoarseRootC=9.00,eventInputC=24.00,eventInputN=0.00 2024 70 irrig eventSoilWater=2.50,eventEvap=2.50 -2024 200 harv eventSoilC=2.74,eventLitterC=1.51,eventLeafC=-1.39,eventWoodC=-1.63,eventFineRootC=-2.52,eventCoarseRootC=-2.97,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=4.25,eventOutputN=0.00 +2024 200 harv eventSoilC=2.74,eventLitterC=1.51,eventLeafC=-1.39,eventWoodC=-1.63,eventAccountingC=-0.00,eventFineRootC=-2.52,eventCoarseRootC=-2.97,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=4.25,eventOutputN=0.00 diff --git a/tests/sipnet/test_events_infrastructure/events_output_no_header.out b/tests/sipnet/test_events_infrastructure/events_output_no_header.out index e8465063..3c727e29 100644 --- a/tests/sipnet/test_events_infrastructure/events_output_no_header.out +++ b/tests/sipnet/test_events_infrastructure/events_output_no_header.out @@ -1,6 +1,6 @@ 2023 65 plant eventLeafC=10.00,eventWoodC=5.00,eventFineRootC=4.00,eventCoarseRootC=3.00,eventInputC=22.00,eventInputN=0.00 2023 70 irrig eventSoilWater=5.00,eventEvap=0.00 -2023 200 harv eventSoilC=12.40,eventLitterC=0.00,eventLeafC=-4.80,eventWoodC=-3.20,eventFineRootC=-4.80,eventCoarseRootC=-4.80,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=5.20,eventOutputN=0.00 +2023 200 harv eventSoilC=12.40,eventLitterC=0.00,eventLeafC=-4.80,eventWoodC=-3.20,eventAccountingC=-0.00,eventFineRootC=-4.80,eventCoarseRootC=-4.80,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=5.20,eventOutputN=0.00 2024 65 plant eventLeafC=10.00,eventWoodC=5.00,eventFineRootC=4.00,eventCoarseRootC=3.00,eventInputC=22.00,eventInputN=0.00 2024 70 irrig eventSoilWater=2.50,eventEvap=2.50 -2024 200 harv eventSoilC=12.14,eventLitterC=0.00,eventLeafC=-10.32,eventWoodC=-5.88,eventFineRootC=-2.88,eventCoarseRootC=-2.48,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=9.42,eventOutputN=0.00 +2024 200 harv eventSoilC=12.14,eventLitterC=0.00,eventLeafC=-10.32,eventWoodC=-5.88,eventAccountingC=-0.00,eventFineRootC=-2.88,eventCoarseRootC=-2.48,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=9.42,eventOutputN=0.00 diff --git a/tests/sipnet/test_events_types/testEventHarvest.c b/tests/sipnet/test_events_types/testEventHarvest.c index 9dbc98f4..1cb089fd 100644 --- a/tests/sipnet/test_events_types/testEventHarvest.c +++ b/tests/sipnet/test_events_types/testEventHarvest.c @@ -140,9 +140,39 @@ int run(void) { return status; } +int checkAccountingHarvest(void) { + int status = 0; + for (int i = -1; i <= 1; i++) { + prepTypesTest(); + updateIntContext("litterPool", 1, CTX_TEST); + updateIntContext("nitrogenCycle", 1, CTX_TEST); + initEnv(); + envi.plantCAccountingDelta = i; + updateBalanceTrackerPreUpdate(); + initEvents("events_one_harvest.in", "events.out", 0); + setupEvents(); + procEvents(); + closeEventOutFile(); + updateBalanceTrackerPostUpdate(); + const double dc = balanceTracker.postTotalC - balanceTracker.preTotalC + + fluxes.eventOutputC * climate->length; + const double dn = balanceTracker.postTotalN - balanceTracker.preTotalN + + fluxes.eventOutputN * climate->length; + if (fabs(dc) > 1e-12 || fabs(dn) > 1e-12 || + fabs(envi.plantWoodC - 1.8) > 1e-12 || + fabs(envi.plantCAccountingDelta - 0.6 * i) > 1e-12) { + logTest("Accounting harvest failed: A=%d, C residual=%.15g, N " + "residual=%.15g\n", + i, dc, dn); + status = 1; + } + } + return status; +} + int main(void) { logTest("Starting run()\n"); - int status = run(); + int status = run() | checkAccountingHarvest(); if (status) { logTest("FAILED testEventHarvest with status %d\n", status); exit(status); From ef504f7394a9364880cc0ded7c4d91a5d14f0744 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Mon, 21 Sep 2026 14:54:08 -0400 Subject: [PATCH 02/14] Interim commit --- src/sipnet/balance.c | 6 +- src/sipnet/events.c | 38 +- src/sipnet/events.h | 12 +- src/sipnet/restart.c | 12 +- src/sipnet/sipnet.c | 69 ++- tests/sipnet/test_modeling/Makefile | 11 +- tests/sipnet/test_modeling/termination.clim | 3 + tests/sipnet/test_modeling/termination.param | 78 +++ .../test_modeling/testCompleteHarvest.c | 494 ++++++++++++++++++ 9 files changed, 681 insertions(+), 42 deletions(-) create mode 100644 tests/sipnet/test_modeling/termination.clim create mode 100644 tests/sipnet/test_modeling/termination.param create mode 100644 tests/sipnet/test_modeling/testCompleteHarvest.c diff --git a/src/sipnet/balance.c b/src/sipnet/balance.c index 1c43608c..2efe55ce 100644 --- a/src/sipnet/balance.c +++ b/src/sipnet/balance.c @@ -47,7 +47,8 @@ void updateBalanceTrackerPostClamp(void) { balanceTracker.clampedC = balanceTracker.finalC - balanceTracker.postTotalC; if (balanceTracker.clampedC < -EPS) { // This shouldn't happen, by construction - logInternalError("Non-negative clamping has cause carbon loss\n"); + logInternalError("Non-negative clamping has cause carbon loss %f\n", + balanceTracker.clampedC); } if (balanceTracker.clampedC < EPS) { balanceTracker.clampedC = 0; @@ -56,7 +57,8 @@ void updateBalanceTrackerPostClamp(void) { balanceTracker.clampedN = balanceTracker.finalN - balanceTracker.postTotalN; if (balanceTracker.clampedN < -EPS) { // This shouldn't happen, by construction - logInternalError("Non-negative clamping has cause nitrogen loss\n"); + logInternalError("Non-negative clamping has cause nitrogen loss %f\n", + balanceTracker.clampedN); } if (balanceTracker.clampedN < EPS) { balanceTracker.clampedN = 0; diff --git a/src/sipnet/events.c b/src/sipnet/events.c index 3d2250f0..3b6d18b3 100644 --- a/src/sipnet/events.c +++ b/src/sipnet/events.c @@ -58,7 +58,13 @@ EventNode *createEventNode(int year, int day, int eventType, // Validate the params if ((fracRA + fracTA > 1) || (fracRB + fracTB > 1)) { logError("invalid harvest newEvent for year %d day %d; above and below " - "must each add to 1 or less", + "must each add to 1 or less\n", + year, day); + exit(EXIT_CODE_BAD_PARAMETER_VALUE); + } + if (fracRA < 0.0 || fracRB < 0.0 || fracTA < 0.0 || fracTB < 0.0) { + logError("invalid harvest newEvent for year %d day %d; fractions must " + "be non-negative\n", year, day); exit(EXIT_CODE_BAD_PARAMETER_VALUE); } @@ -465,8 +471,7 @@ void processEvents(void) { } // Reset harvest tracking - eventTrackers.harvestFracRemoved = 0; - eventTrackers.harvestFracTransferred = 0; + eventTrackers.harvestTrackers = (HarvestTrackers){0}; while (gEvent != NULL && gEvent->year <= climYear && gEvent->day <= climDay) { // The events file has been tested on read, so we know this event list @@ -545,8 +550,8 @@ void processEvents(void) { // pools const HarvestParams *harvParams = gEvent->eventParams; const double fracRA = harvParams->fractionRemovedAbove; - const double fracTA = harvParams->fractionTransferredAbove; const double fracRB = harvParams->fractionRemovedBelow; + const double fracTA = harvParams->fractionTransferredAbove; const double fracTB = harvParams->fractionTransferredBelow; const double woodC = envi.plantWoodC + envi.plantCAccountingDelta; @@ -554,11 +559,32 @@ void processEvents(void) { double aboveMass = woodC + envi.plantLeafC; double belowMass = envi.fineRootC + envi.coarseRootC; double totalMass = aboveMass + belowMass; + HarvestTrackers *ht = &eventTrackers.harvestTrackers; if (totalMass > TINY) { double massRemoved = fracRA * aboveMass + fracRB * belowMass; double massTransferred = fracTA * aboveMass + fracTB * belowMass; - eventTrackers.harvestFracRemoved += massRemoved / totalMass; - eventTrackers.harvestFracTransferred += massTransferred / totalMass; + ht->totalFracRemoved += massRemoved / totalMass; + ht->totalFracTransferred += massTransferred / totalMass; + ht->totalFracRemovedAbove += fracRA; + ht->totalFracRemovedBelow += fracRB; + ht->totalFracTransferredAbove += fracTA; + ht->totalFracTransferredBelow += fracTB; + if (ht->totalFracRemovedAbove + ht->totalFracTransferredAbove > + 1.0 + TINY || + ht->totalFracRemovedBelow + ht->totalFracTransferredBelow > + 1.0 + TINY) { + logError("Harvest event(s) at year %d day %d has total above-ground" + " or below-ground removal + transfer fraction > 1.0" + " (above %.3f, below %.3f)\n", + gEvent->year, gEvent->day, + ht->totalFracRemovedAbove + ht->totalFracTransferredAbove, + ht->totalFracRemovedBelow + ht->totalFracTransferredBelow); + exit(EXIT_CODE_BAD_PARAMETER_VALUE); + } + } else { + logWarning("Harvest event at year %d day %d has no biomass to remove " + "or transfer\n", + gEvent->year, gEvent->day); } // Litter increase diff --git a/src/sipnet/events.h b/src/sipnet/events.h index 22a43fcb..d829f6bb 100644 --- a/src/sipnet/events.h +++ b/src/sipnet/events.h @@ -209,6 +209,15 @@ void updatePoolsForEvents(void); */ void freeEventList(void); +typedef struct HarvestTrackersStruct { + double totalFracRemoved; + double totalFracTransferred; + double totalFracRemovedAbove; + double totalFracRemovedBelow; + double totalFracTransferredAbove; + double totalFracTransferredBelow; +} HarvestTrackers; + // Variables to track events with lingering effects typedef struct EventTrackerStruct { // Tillage effect on Rh; exponentially decays at each time step by a factor @@ -216,8 +225,7 @@ typedef struct EventTrackerStruct { double d_till_mod; // Fraction removed and transferred from harvest event this time step - double harvestFracRemoved; - double harvestFracTransferred; + HarvestTrackers harvestTrackers; } EventTrackers; extern EventTrackers eventTrackers; diff --git a/src/sipnet/restart.c b/src/sipnet/restart.c index e8c2a4b6..78689e22 100644 --- a/src/sipnet/restart.c +++ b/src/sipnet/restart.c @@ -37,7 +37,7 @@ #define NUM_TRACKER_FIELDS 33 #define NUM_PHENOLOGY_TRACKERS_FIELDS 3 #define NUM_SURVIVAL_TRACKERS_FIELDS 1 -#define NUM_EVENT_TRACKERS_FIELDS 3 +#define NUM_EVENT_TRACKERS_FIELDS 7 // This one shouldn't change #define NUM_END_FIELDS 1 @@ -291,9 +291,13 @@ void initResetState(RestartState *state, MeanTracker *npp) { } ind = 0; - state->eventPF[ind++] = (StateField){"event_trackers.d_till_mod", FT_DOUBLE, &eventTrackers.d_till_mod, 0}; - state->eventPF[ind++] = (StateField){"event_trackers.harvestFracRemoved", FT_DOUBLE, &eventTrackers.harvestFracRemoved, 0}; - state->eventPF[ind++] = (StateField){"event_trackers.harvestFracTransferred", FT_DOUBLE, &eventTrackers.harvestFracTransferred, 0}; + state->eventPF[ind++] = (StateField){"event_trackers.d_till_mod", FT_DOUBLE, &eventTrackers.d_till_mod, 0}; + state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracRemoved", FT_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemoved, 0}; + state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracTransferred", FT_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferred, 0}; + state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracRemovedAbove", FT_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemovedAbove, 0}; + state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracRemovedBelow", FT_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemovedBelow, 0}; + state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracTransferredAbove", FT_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferredAbove, 0}; + state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracTransferredBelow", FT_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferredBelow, 0}; state->eventPF[ind++] = (StateField){"event_trackers.invalid", FT_INVALID, NULL, FIELD_INVALID}; if (ind != NUM_EVENT_TRACKERS_FIELDS + 1) { logInternalError("Restart array size mismatch: eventPF\n"); diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 1de3f7cc..e54aaacb 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -1532,7 +1532,9 @@ void initPhenologyTrackers(void) { // Check that woodC and total root C are both positive and that there was no // terminating harvest int hasSufficientBiomass(void) { - if (eventTrackers.harvestFracRemoved + eventTrackers.harvestFracTransferred > + // If there was a harvest termination event, the answer is no + if (eventTrackers.harvestTrackers.totalFracRemoved + + eventTrackers.harvestTrackers.totalFracTransferred > 1.0 - TINY) { return 0; } @@ -1712,18 +1714,15 @@ void checkForMortality(void) { plantSurvivalTracker.isAlive = 0; double totalWoodC = getTotalWoodC(); double totalRootC = envi.fineRootC + envi.coarseRootC; - - if (eventTrackers.harvestFracRemoved + - eventTrackers.harvestFracTransferred >= - TINY) { + HarvestTrackers *ht = &eventTrackers.harvestTrackers; + int harvestOccurred = + ht->totalFracRemoved + ht->totalFracTransferred >= TINY; + if (harvestOccurred) { logInfo("Plant mortality detected after harvest event: total fraction " - "removed %.3f total fraction transferred %.3f; woodC %f " - "totalWoodC %f coarseRootC %f fineRootC %f year %d day %d " + "removed %.3f total fraction transferred %.3f on year %d day %d " "time %6.3f; zeroing out biomass pools\n", - eventTrackers.harvestFracRemoved, - eventTrackers.harvestFracTransferred, envi.plantWoodC, totalWoodC, - envi.coarseRootC, envi.fineRootC, climate->year, climate->day, - climate->time); + ht->totalFracRemoved, ht->totalFracTransferred, climate->year, + climate->day, climate->time); } else { logWarning( "Plant mortality detected as wood or total root carbon is zero " @@ -1738,19 +1737,34 @@ void checkForMortality(void) { // multiple processes being modeled, it's believable that there may be a bit // of overshoot when a plant dies - for example, a 100% harvest event with // any overall loss (respiration, turnover). - envi.soilC += totalRootC; + + // Also, if there was a harvest, reduce by the appropriate removal fraction + // (For no harvest, these will be 1) + double aboveRemovalReduction = (1 - ht->totalFracRemovedAbove); + double belowRemovalReduction = (1 - ht->totalFracRemovedBelow); + + envi.soilC += totalRootC * belowRemovalReduction; + double aboveC = + envi.plantWoodC + envi.plantLeafC + envi.plantCAccountingDelta; if (ctx.litterPool) { - envi.litterC += - envi.plantWoodC + envi.plantLeafC + envi.plantCAccountingDelta; + envi.litterC += aboveC * aboveRemovalReduction; } else { - envi.soilC += - envi.plantWoodC + envi.plantLeafC + envi.plantCAccountingDelta; + envi.soilC += aboveC * aboveRemovalReduction; } + fluxes.eventOutputC += (aboveC * ht->totalFracRemovedAbove + + totalRootC * ht->totalFracRemovedBelow) / + climate->length; + if (ctx.nitrogenCycle) { // litter pool implied - envi.soilOrgN += + double aboveN = + envi.plantWoodC / params.woodCN + envi.plantLeafC / params.leafCN; + double belowN = envi.fineRootC / params.fineRootCN + envi.coarseRootC / params.woodCN; - envi.litterN += envi.plantWoodC / params.woodCN + - envi.plantLeafC / params.leafCN + envi.plantStorageN; + envi.soilOrgN += belowN * belowRemovalReduction; + envi.litterN += aboveN * aboveRemovalReduction + envi.plantStorageN; + fluxes.eventOutputN += (aboveN * ht->totalFracRemovedAbove + + belowN * ht->totalFracRemovedBelow) / + climate->length; } // Force pools to zero @@ -1766,11 +1780,11 @@ void checkForMortality(void) { resetMeanTracker(meanNPP, 0.0); if (ctx.events) { - writeComputedEventOut( - climate->year, climate->day, eventTypeToString(PLANTDEATH), 4, - "harvestFracRemoved", eventTrackers.harvestFracRemoved, - "harvestFracTransferred", eventTrackers.harvestFracTransferred, - "totalWoodC", totalWoodC, "totalRootC", totalRootC); + writeComputedEventOut(climate->year, climate->day, + eventTypeToString(PLANTDEATH), 4, + "harvestFracRemoved", ht->totalFracRemoved, + "harvestFracTransferred", ht->totalFracTransferred, + "totalWoodC", totalWoodC, "totalRootC", totalRootC); } } } @@ -1797,12 +1811,13 @@ void updatePoolsAndBalance() { updateNitrogenPools(); } + // Check for obvious plant death (wood and/or roots at zero); also + // correct for harvest termination + checkForMortality(); + // Calc total C and N after pool updates updateBalanceTrackerPostUpdate(); - // Check for obvious plant death (wood and/or roots at zero) - checkForMortality(); - // Verify none of our stocks have gone negative (set any that are to zero). ensureNonNegativeStocks(); diff --git a/tests/sipnet/test_modeling/Makefile b/tests/sipnet/test_modeling/Makefile index 5cb4a25b..5ab3bfb5 100644 --- a/tests/sipnet/test_modeling/Makefile +++ b/tests/sipnet/test_modeling/Makefile @@ -8,7 +8,16 @@ LDFLAGS=-L$(ROOT_DIR)/libs LDLIBS=-lsipnet -lsipnet_common -lm # List test files in this directory here -TEST_CFILES=testNitrogenCycle.c testDependencyFunctions.c testBalance.c testMethane.c testSoilMoisture.c testCarbonSaturation.c testPlantMortality.c testFluxCalculations.c +TEST_CFILES= \ + testBalance.c \ + testCarbonSaturation.c \ + testCompleteHarvest.c \ + testDependencyFunctions.c \ + testFluxCalculations.c \ + testMethane.c \ + testNitrogenCycle.c \ + testPlantMortality.c \ + testSoilMoisture.c \ # The rest is boilerplate, likely copyable as is to a new test directory TEST_OBJ_FILES=$(TEST_CFILES:%.c=%.o) diff --git a/tests/sipnet/test_modeling/termination.clim b/tests/sipnet/test_modeling/termination.clim new file mode 100644 index 00000000..79baf35a --- /dev/null +++ b/tests/sipnet/test_modeling/termination.clim @@ -0,0 +1,3 @@ +2017 19 21 0.125 10.574335 8.893 6.502 3.059 256.369 119.4744 1025.2 3.23371 +2017 20 0 0.125 10.85354 8.879 5.898 2.163 281.998 119.9049 1023.9 2.411573 +2017 20 3 0.125 8.795038 8.862 0.4043 1.941 167.848 173.0103 967.9 2.896282 diff --git a/tests/sipnet/test_modeling/termination.param b/tests/sipnet/test_modeling/termination.param new file mode 100644 index 00000000..fe8fdd45 --- /dev/null +++ b/tests/sipnet/test_modeling/termination.param @@ -0,0 +1,78 @@ +plantWoodInit 70.0 +laiInit 2.3481421625603751 +leafCSpWt 29.409172953079818 +litterInit 525.0 +soilInit 5550.0 +soilWFracInit 1.9282707581831169 +snowInit 0 +plantStorageNInit 0.001 +fineRootFrac 0.25 +coarseRootFrac 0.055 +aMax 603.51460591037255 +aMaxFrac 0.75 +psnTMin 2 +psnTOpt 17.796068268198724 +dVpdSlope 0.015793854979347089 +dVpdExp 2 +halfSatPar 15.300351381651099 +attenuation 0.57999999999999996 +baseVegResp 0.0060000000000000001 +baseFolRespFrac 0.016122588146105607 +baseSoilResp 0.030052329942810149 +baseFineRootResp 0.2695450954857716 +baseCoarseRootResp 0.0060000000000000001 +vegRespQ10 1.4019202256757124 +fineRootQ10 2.6000000000000001 +coarseRootQ10 3.20392084014602 +soilRespQ10 1.6759340808079286 +growthRespFrac 0.20000000000000001 +frozenSoilFolREff 0 +frozenSoilThreshold 0 +soilRespMoistEffect 1 +leafOnDay 0 +gddLeafOn 500 +soilTempLeafOn 12 +leafOnReallocFrac 0.10337817506913401 +leafOffDay 0 +leafGrowth 126 +fracLeafFall 1 +woodTurnoverRate 0.014 +leafTurnoverRate 1.09744012449142 +fineRootTurnoverRate 0.2682821481714458 +coarseRootTurnoverRate 0.056000000000000001 +litterBreakdownRate 0.56007424802511896 +fracLitterRespired 0.76768911631115766 +fineRootAllocation 0.10609482754666355 +woodAllocation 0.33470604289674222 +leafAllocation 0.53456185029321135 +waterRemoveFrac 0.087999999999999995 +frozenSoilEff 1 +wueConst 10.9 +soilWHC 12 +immedEvapFrac 0.10000000000000001 +leafPoolDepth 0.10000000000000001 +fastFlowFrac 0 +snowMelt 0.14999999999999999 +rdConst 1149.9961760992301 +rSoilConst1 8.1999999999999993 +rSoilConst2 4.2999999999999998 +cFracLeaf 0.46999999999999997 +mineralNInit 0 +soilOrgNInit 528.09892365915437 +litterOrgNInit 18.860086734898815 +nVolatilizationFrac 6.8860489049188205e-05 +nLeachingFrac 0.248596695361815 +leafNResorptionFrac 0.53599062060008895 +leafCN 15.0 +woodCN 77.0 +fineRootCN 40.0 +kCN 190.0 +nFixationFracMax 0 +halfNFixationMax 0.99630850231037704 +fAnoxia 0.29999999999999999 +anaerobicDecompRate 0.26666853510696198 +anaerobicTransExp 9.9741824548147182 +soilMethaneRate 1.0029297483187848e-05 +litterMethaneRate 4.9804017375032515e-06 +waterDrainFrac 0.10000000000000001 +soilCSaturation 1.0 diff --git a/tests/sipnet/test_modeling/testCompleteHarvest.c b/tests/sipnet/test_modeling/testCompleteHarvest.c new file mode 100644 index 00000000..e0e22e35 --- /dev/null +++ b/tests/sipnet/test_modeling/testCompleteHarvest.c @@ -0,0 +1,494 @@ +// Native three-step fixture from sipnet-termination-repro.tar.gz, derived from +// wet_d0.1_m1, point 332581. The first step primes mean NPP; positive PAR is +// intentional even at the midnight harvest timestamp. This is a constructed +// model regression fixture, not a field observation or restart replay. +#include "utils/tUtils.h" +#include "sipnet/events.c" +#include "sipnet/sipnet.c" +#include + +static int failures; +static ModelParams *testParams; + +static void near(double actual, double expected, const char *name) { + if (!isfinite(actual) || fabs(actual - expected) > 1e-8) { + logTest("%s: %.15g, expected %.15g\n", name, actual, expected); + failures++; + } +} + +static void balanced(void) { + near(balanceTracker.deltaC, 0, "C balance"); + near(balanceTracker.deltaN, 0, "N balance"); + near(balanceTracker.clampedC, 0, "C clamping"); + near(balanceTracker.clampedN, 0, "N clamping"); +} + +static void start(int mode, const char *harvest) { + initContext(); + ctx.events = 1; + ctx.litterPool = mode > 0; + ctx.nitrogenCycle = mode == 2; + ctx.anaerobic = 1; + ctx.waterHResp = 1; + ctx.growthResp = 0; + ctx.leafWater = 0; + ctx.flooding = 1; + ctx.gdd = 0; + ctx.soilPhenol = 0; + envi = (Envi){0}; + fluxes = (Fluxes){0}; + eventTrackers = (EventTrackers){0}; + initModel(&testParams, "termination.param", "termination.clim"); + setupModel(); + gEvents = harvest ? createEventNode(2017, 20, HARVEST, harvest) : NULL; + setupEvents(); + openEventOutFile("events.out", 1); + updateState(); // Prime mean NPP with the native first timestep. + balanced(); + climate = climate->nextClim; +} + +static void finish(void) { + cleanupModel(); + deleteModelParams(testParams); + gEvents = gEvent = NULL; +} + +static void clearPlant(void) { + near(envi.plantLeafC, 0, "leaf"); + near(envi.plantWoodC, 0, "wood"); + near(envi.fineRootC, 0, "fine roots"); + near(envi.coarseRootC, 0, "coarse roots"); + near(envi.plantCAccountingDelta, 0, "accounting C"); + near(envi.plantStorageN, 0, "storage N"); + near(plantSurvivalTracker.isAlive, 0, "dead"); + near(getMeanTrackerMean(meanNPP), 0, "mean NPP"); +} + +static void fullCase(int mode, int accounting, int dark, int limited, + const char *harvest, double aboveExport, + double belowExport) { + // A no-harvest run supplies the independently calculated end-of-step pools. + logTest("Running full case with mode %d accounting %d dark %d limited %d" + " harvest %s aboveExport %.2f belowExport %.2f\n", + mode, accounting, dark, limited, harvest ? harvest : "none", + aboveExport, belowExport); + Envi end = {0}; + for (int h = 0; h < 2; h++) { + start(mode, h ? harvest : NULL); + envi.plantCAccountingDelta = accounting; + if (dark) + climate->par = 0; + if (limited && mode == 2) { + envi.minN = 0; + envi.plantStorageN = 0; + } + updateState(); + balanced(); + if (!h) { + // First time through, capture envi state + end = envi; + // logTest("(end) nLeach %f nVol %f\n", + // fluxes.nLeaching * climate->length, fluxes.nVolatilization * + // climate->length); + // logTest("(end) fluxes.eventOutputN %f\n", fluxes.eventOutputN * + // climate->length); + } else { + // logTest("(envi) fluxes.eventOutputN %f\n", fluxes.eventOutputN * + // climate->length); + const double aboveC = + end.plantLeafC + end.plantWoodC + end.plantCAccountingDelta; + const double belowC = end.fineRootC + end.coarseRootC; + const double aboveN = mode == 2 ? end.plantLeafC / params.leafCN + + end.plantWoodC / params.woodCN + : 0; + const double belowN = mode == 2 ? end.fineRootC / params.fineRootCN + + end.coarseRootC / params.woodCN + : 0; + near(envi.soilC, + end.soilC + belowC * (1 - belowExport) + + (mode == 0 ? aboveC * (1 - aboveExport) : 0), + "soil routing"); + near(envi.litterC, end.litterC + (mode ? aboveC * (1 - aboveExport) : 0), + "litter routing"); + near(fluxes.eventOutputC * climate->length, + aboveC * aboveExport + belowC * belowExport, "C export"); + near(fluxes.eventOutputN * climate->length, + aboveN * aboveExport + belowN * belowExport, "N export"); + // logTest("N export: above %f fracAbove %f below %f fracbelow %f " + // "eventOutputN %f\n", + // aboveN, aboveExport, belowN, belowExport, fluxes.eventOutputN * + // climate->length); + + if (mode == 2) { + near(envi.litterN, + end.litterN + aboveN * (1 - aboveExport) + end.plantStorageN, + "litter N including storage"); + near(envi.soilOrgN, end.soilOrgN + belowN * (1 - belowExport), + "soil N"); + } + clearPlant(); + climate = climate->nextClim; + updateState(); + balanced(); + clearPlant(); + near(fluxes.photosynthesis, 0, "no fallow photosynthesis"); + // A later planting can restore living tissue. + EventNode *plant = createEventNode(2017, 21, PLANTING, "2 3 4 5"); + gEvents->nextEvent = plant; + gEvent = plant; + climate->day = 21; + updateState(); + balanced(); + near(plantSurvivalTracker.isAlive, 1, "alive after planting"); + } + finish(); + } + // if (mode == 2) { + // logTest("Press return to continue\n"); + // getchar(); + // } +} + +static void partialCase(int accounting, const char *harvest, double fraction) { + start(2, harvest); + envi.plantCAccountingDelta = accounting; + Envi before = envi; + resetFluxes(); + processEvents(); + updateBalanceTrackerPreUpdate(); + updatePoolsForEvents(); + // near(updatePoolsForFullHarvest(), 0, "partial harvest not deferred"); + near(envi.plantWoodC, before.plantWoodC * (1 - fraction), "partial wood"); + near(envi.plantCAccountingDelta, accounting * (1 - fraction), + "partial accounting"); + updateBalanceTrackerPostUpdate(); + near(balanceTracker.postTotalC + fluxes.eventOutputC * climate->length, + balanceTracker.preTotalC, "partial C conservation"); + near(balanceTracker.postTotalN + fluxes.eventOutputN * climate->length, + balanceTracker.preTotalN, "partial N conservation"); + finish(); +} + +static void partialTimestepCase(int accounting) { + start(2, ".1 .2 .3 .4"); + envi.plantCAccountingDelta = accounting; + updateState(); + balanced(); + near(plantSurvivalTracker.isAlive, 1, "partial harvest remains alive"); + climate = climate->nextClim; + updateState(); + balanced(); + near(plantSurvivalTracker.isAlive, 1, "partial harvest can continue growth"); + finish(); +} + +static void invalidCase(int which) { + pid_t child = fork(); + if (child < 0) { + failures++; + return; + } + if (child == 0) { + start(2, "0 0 1 1"); + if (which < 4) { // 0, 1, 2, 3 + const char *bad[] = {"0 0 -1 1", "-0.1 0 1 1", "0 0 1.1 1", "1 0.1 0 1"}; + createEventNode(2017, 20, HARVEST, bad[which]); + } else { // 4, 5 + EventNode *extra = createEventNode(2017, 20, HARVEST, "0 0 .2 .2"); + if (which < 5) + gEvents->nextEvent = extra; + else { + extra->nextEvent = gEvents; + gEvents = extra; + setupEvents(); + } + processEvents(); + } + _exit(99); + } + int status; + waitpid(child, &status, 0); + int expected = + which < 6 ? EXIT_CODE_BAD_PARAMETER_VALUE : EXIT_CODE_INTERNAL_ERROR; + if (!WIFEXITED(status) || WEXITSTATUS(status) != expected) { + logTest("Invalid case %d did not exit with code %d\n", which, expected); + failures++; + } +} + +static void coincidentFertilizerCase(void) { + Envi expected = {0}; + for (int h = 0; h < 2; h++) { + start(2, h ? "0 0 1 1" : NULL); + EventNode *fert = createEventNode(2017, 20, FERTILIZATION, "1 2 3"); + fert->nextEvent = gEvents; + gEvents = fert; + setupEvents(); + updateState(); + balanced(); + if (!h) + expected = envi; + else { + near(envi.soilC, + expected.soilC + expected.fineRootC + expected.coarseRootC, + "fertilizer plus harvest soil C"); + near(envi.litterC, + expected.litterC + expected.plantWoodC + expected.plantLeafC + + expected.plantCAccountingDelta, + "fertilizer plus harvest litter C"); + near(envi.minN, expected.minN, "fertilizer applied once"); + clearPlant(); + } + finish(); + } +} + +static void sameEnvi(const Envi *a, const Envi *b) { +#define E(f) near(a->f, b->f, #f) + E(plantWoodC); + E(plantLeafC); + E(soilC); + E(soilWater); + E(litterC); + E(snow); + E(coarseRootC); + E(fineRootC); + E(minN); + E(soilOrgN); + E(litterN); + E(plantStorageN); + E(plantCAccountingDelta); +#undef E +} +static void coincidentEventCase(int type, const char *arguments) { + Envi end; + Fluxes ordinary; + for (int order = 0; order < 3; order++) { + start(2, order ? ".25 .75 .75 .25" : NULL); + params.leafGrowth = 1; + params.fracLeafFall = .25; + EventNode *extra = createEventNode(2017, 20, type, arguments); + if (order == 2) + gEvents->nextEvent = extra; + else { + extra->nextEvent = gEvents; + gEvents = extra; + } + setupEvents(); + updateState(); + balanced(); + if (!order) { + end = envi; + ordinary = fluxes; + if ((type == LEAFON && fluxes.eventLeafOnCreation <= 0) || + (type == LEAFOFF && fluxes.eventLeafOffLitter <= 0) || + (type == IRRIGATION && fluxes.eventSoilWater <= 0)) + failures++; + } else { + double dt = climate->length; + double ac = end.plantWoodC + end.plantCAccountingDelta + end.plantLeafC; + double bc = end.fineRootC + end.coarseRootC; + double an = + end.plantWoodC / params.woodCN + end.plantLeafC / params.leafCN; + double bn = + end.fineRootC / params.fineRootCN + end.coarseRootC / params.woodCN; + Envi expected = end; + expected.soilC += .25 * bc; + expected.litterC += .75 * ac; + expected.soilOrgN += .25 * bn; + expected.litterN += .75 * an + end.plantStorageN; + expected.plantWoodC = expected.plantLeafC = expected.fineRootC = 0; + expected.coarseRootC = expected.plantCAccountingDelta = + expected.plantStorageN = 0; + sameEnvi(&envi, &expected); + clearPlant(); + near(fluxes.eventInputC, ordinary.eventInputC, "C input once"); + near(fluxes.eventInputN, ordinary.eventInputN, "N input once"); + near(fluxes.eventOutputC, + ordinary.eventOutputC + (.25 * ac + .75 * bc) / dt, "C export"); + near(fluxes.eventOutputN, + ordinary.eventOutputN + (.25 * an + .75 * bn) / dt, "N export"); +#define F(f) near(fluxes.f, ordinary.f, #f) + F(eventSoilWater); + F(eventEvap); + F(eventLeafOnCreation); + F(eventLeafOnCreationFromWood); + F(eventLeafOffLitter); + F(eventLeafOffNResorption); +#undef F + climate = climate->nextClim; + updateState(); + balanced(); + clearPlant(); + near(fluxes.photosynthesis, 0, "no fallow photosynthesis"); + } + finish(); + } +} + +static void restartCase(void) { + Envi expected; + for (int resumed = 0; resumed < 2; resumed++) { + start(2, "0 0 1 1"); + if (resumed) { + restartNoteProcessedClimateStep(firstClimate); + restartWriteCheckpoint("termination.restart", meanNPP); + envi = (Envi){0}; + resetMeanTracker(meanNPP, 0); + restartLoadCheckpoint("termination.restart", meanNPP); + } + updateState(); + balanced(); + clearPlant(); + if (resumed) { + restartNoteProcessedClimateStep(climate); + restartWriteCheckpoint("termination.restart", meanNPP); + } + climate = climate->nextClim; + if (resumed) { + envi = (Envi){0}; + restartLoadCheckpoint("termination.restart", meanNPP); + } + updateState(); + balanced(); + clearPlant(); + if (!resumed) + expected = envi; + else { + sameEnvi(&envi, &expected); + near(envi.soilC, expected.soilC, "restart soil C"); + near(envi.litterC, expected.litterC, "restart litter C"); + near(envi.minN, expected.minN, "restart mineral N"); + near(envi.soilOrgN, expected.soilOrgN, "restart soil N"); + near(envi.litterN, expected.litterN, "restart litter N"); + } + finish(); + } +} + +static void leafBudgetCases(void) { + for (int kind = 0; kind < 4; kind++) + for (int dark = 0; dark < 2; dark++) + for (int sign = -1; sign <= 1; sign++) + for (int account = -1; account <= 1; account++) + for (int limited = 0; limited < 2; limited++) + for (int resorb = 0; resorb < 3; resorb++) + for (int harvest = 0; harvest < 2; harvest++) { + start(2, harvest ? "0 0 1 1" : NULL); + envi.plantLeafC = 1; + envi.plantWoodC = envi.fineRootC = envi.coarseRootC = 100; + envi.plantCAccountingDelta = account; + envi.plantStorageN = 0; + envi.minN = limited ? 0 : 1000; + params.leafAllocation = params.woodAllocation = .25; + params.fineRootAllocation = params.coarseRootAllocation = .25; + params.leafTurnoverRate = 2; + params.woodTurnoverRate = params.fineRootTurnoverRate = 0; + params.coarseRootTurnoverRate = params.litterBreakdownRate = 0; + params.baseSoilResp = params.soilMethaneRate = + params.litterMethaneRate = 0; + params.nVolatilizationFrac = params.nLeachingFrac = 0; + params.nFixationFracMax = 0; + params.leafNResorptionFrac = .5 * resorb; + params.fracLeafFall = kind == 0 ? 1 : kind == 1 ? .25 : .75; + if (dark) + climate->par = 0; + resetMeanTracker(meanNPP, 20 * sign); + if (kind != 3) { + EventNode *off = createEventNode(2017, 20, LEAFOFF, ""); + off->nextEvent = gEvents; + gEvents = off; + if (kind == 2) { + EventNode *second = createEventNode(2017, 20, LEAFOFF, ""); + second->nextEvent = off->nextEvent; + off->nextEvent = second; + } + setupEvents(); + } + double beforeLitter = envi.litterC; + double beforeLitterN = envi.litterN; + updateState(); + balanced(); + double shed = (fluxes.leafLitter + fluxes.eventLeafOffLitter) * + climate->length; + if (shed < -1e-10 || shed > 1 + 1e-10 || + envi.plantLeafC < -1e-10) { + logTest("Leaf budget overdraw: %.17g, leaf %.17g\n", shed, + envi.plantLeafC); + failures++; + } + double eventExpected = kind == 3 ? 0 : kind == 1 ? .25 : 1; + near(fluxes.eventLeafOffLitter * climate->length, eventExpected, + "event leaves conserved"); + if (!harvest) { + near(envi.litterC - beforeLitter, shed, + "leaf litter transfer"); + near(envi.litterN - beforeLitterN, + shed * (1 - params.leafNResorptionFrac) / params.leafCN, + "leaf litter N transfer"); + } else { + clearPlant(); + climate = climate->nextClim; + updateState(); + balanced(); + clearPlant(); + near(fluxes.photosynthesis, 0, + "no regrowth after complete harvest"); + } + finish(); + } +} + +// Full export during net loss must not borrow carbon from an empty soil pool. +static void exportLossCase(void) { + start(1, "1 1 0 0"); + envi.litterC = envi.soilC = envi.litterN = envi.soilOrgN = 0; + params.baseSoilResp = params.litterBreakdownRate = 0; + params.soilMethaneRate = params.litterMethaneRate = 0; + climate->par = 0; + resetMeanTracker(meanNPP, -20); + updateState(); + balanced(); + clearPlant(); + finish(); +} + +int main(void) { + const char *harvest[] = {"0 0 1 1", "1 1 0 0", ".25 .75 .75 .25"}; + const double above[] = {0, 1, .25}, below[] = {0, 1, .75}; + // for (int mode = 0; mode < 3; mode++) + // for (int account = -1; account <= 1; account++) + // for (int dark = 0; dark < 2; dark++) + // for (int limited = 0; limited <= (mode == 2); limited++) + // for (int route = 0; route < 3; route++) + // fullCase(mode, account, dark, limited, harvest[route], + // above[route], + // below[route]); + + fullCase(0, -1, 0, 0, harvest[0], above[0], below[0]); + + for (int account = -1; account <= 1; account++) { + partialTimestepCase(account); + partialCase(account, ".1 .2 .3 .4", .4); + partialCase(account, "0 0 .999999999 .999999999", .999999999); + } + + for (int which = 0; which < 6; which++) { + invalidCase(which); + } + + logTest("Press return to continue\n"); + getchar(); + + coincidentFertilizerCase(); + coincidentEventCase(LEAFON, ""); + coincidentEventCase(LEAFOFF, ""); + coincidentEventCase(IRRIGATION, "2 0"); + restartCase(); + leafBudgetCases(); + exportLossCase(); + logTest("Complete harvest failures: %d\n", failures); + return failures != 0; +} From 3d13355b8c1825c6951fa8ade286a2b9c40b45fe Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Fri, 25 Sep 2026 13:52:44 -0400 Subject: [PATCH 03/14] Interim commit --- src/sipnet/events.c | 322 +++++++++++++++++- src/sipnet/events.h | 31 +- src/sipnet/limitations.c | 53 ++- src/sipnet/nitrogen.c | 9 +- src/sipnet/sipnet.c | 32 +- src/sipnet/state.c | 4 + src/sipnet/state.h | 7 +- .../test_modeling/testCompleteHarvest.c | 59 +++- 8 files changed, 459 insertions(+), 58 deletions(-) diff --git a/src/sipnet/events.c b/src/sipnet/events.c index 3b6d18b3..66401ede 100644 --- a/src/sipnet/events.c +++ b/src/sipnet/events.c @@ -43,6 +43,8 @@ EventNode *createEventNode(int year, int day, int eventType, newEvent->year = year; newEvent->day = day; newEvent->type = eventType; + newEvent->numLogParamPairs = 0; + newEvent->logLine = NULL; switch (eventType) { case HARVEST: { @@ -452,7 +454,325 @@ int isFirstEventBefore(int year, int day) { return firstEvent->day < day; } -void processEvents(void) { +void processEventsForCarbon(EventNode *event) { + // Event fluxes have all been reset to zero at the start of the time step, + // so we can just add to them as needed + + // If event starts off NULL, this function will just fall through, as it + // should. + const int climYear = climate->year; + const int climDay = climate->day; + const double climLen = climate->length; + + // As this is used as a divisor in many places, let's make sure it's >0 + if (climLen <= 0) { + logError("climate length (%f) on year %d day %d is non-positive; please " + "fix and re-run", + climLen, climYear, climDay); + exit(EXIT_CODE_BAD_PARAMETER_VALUE); + } + + // Reset harvest tracking + eventTrackers.harvestTrackers = (HarvestTrackers){0}; + + while (event != NULL && event->year <= climYear && event->day <= climDay) { + // The events file has been tested on read, so we know this event list + // should be in chrono order. However, we need to check to make sure the + // current event is not in the past, as that would indicate an event that + // did not have a corresponding climate file record. + if (event->year < climYear || event->day < climDay) { + logError("Agronomic event found for year: %d day: %d that does not " + "have a corresponding record in the climate file\n", + event->year, event->day); + exit(EXIT_CODE_INPUT_FILE_ERROR); + } + + switch (event->type) { + case IRRIGATION: { + const IrrigationParams *irrParams = event->eventParams; + const double amount = irrParams->amountAdded; + double soilAmount, evapAmount; + if (irrParams->method == CANOPY) { + // Part of the irrigation evaporates, and the rest makes it to the + // soil. Evaporated fraction: + evapAmount = params.immedEvapFrac * amount; + // Soil fraction: + soilAmount = amount - evapAmount; + } else if (irrParams->method == SOIL) { + // All goes to the soil + evapAmount = 0.0; + soilAmount = amount; + } else { + logError("Unknown irrigation method type: %d\n", irrParams->method); + exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM); + } + fluxes.eventEvap += evapAmount / climLen; + fluxes.eventSoilWater += soilAmount / climLen; + } break; + case PLANTING: { + const PlantingParams *plantParams = event->eventParams; + const double leafC = plantParams->leafC; + const double woodC = plantParams->woodC; + const double fineRootC = plantParams->fineRootC; + const double coarseRootC = plantParams->coarseRootC; + + // Update the fluxes + fluxes.eventLeafC += leafC / climLen; + fluxes.eventWoodC += woodC / climLen; + fluxes.eventFineRootC += fineRootC / climLen; + fluxes.eventCoarseRootC += coarseRootC / climLen; + + // No need to allocate to biomass N pools, we don't track that N + // explicitly + + // MASS BALANCE: this is a system input + const double inputC = leafC + woodC + fineRootC + coarseRootC; + fluxes.eventInputC += inputC / climLen; + } break; + case HARVEST: { + // Harvest can both remove biomass and move biomass to the soil/litter + // pools + const HarvestParams *harvParams = event->eventParams; + const double fracRA = harvParams->fractionRemovedAbove; + const double fracRB = harvParams->fractionRemovedBelow; + const double fracTA = harvParams->fractionTransferredAbove; + const double fracTB = harvParams->fractionTransferredBelow; + const double woodC = envi.plantWoodC + envi.plantCAccountingDelta; + + // Record fraction of total biomass removed and transferred + double aboveMass = woodC + envi.plantLeafC; + double belowMass = envi.fineRootC + envi.coarseRootC; + double totalMass = aboveMass + belowMass; + HarvestTrackers *ht = &eventTrackers.harvestTrackers; + if (totalMass > TINY) { + double massRemoved = fracRA * aboveMass + fracRB * belowMass; + double massTransferred = fracTA * aboveMass + fracTB * belowMass; + ht->totalFracRemoved += massRemoved / totalMass; + ht->totalFracTransferred += massTransferred / totalMass; + ht->totalFracRemovedAbove += fracRA; + ht->totalFracRemovedBelow += fracRB; + ht->totalFracTransferredAbove += fracTA; + ht->totalFracTransferredBelow += fracTB; + if (ht->totalFracRemovedAbove + ht->totalFracTransferredAbove > + 1.0 + TINY || + ht->totalFracRemovedBelow + ht->totalFracTransferredBelow > + 1.0 + TINY) { + logError("Harvest event(s) at year %d day %d has total above-ground" + " or below-ground removal + transfer fraction > 1.0" + " (above %.3f, below %.3f)\n", + event->year, event->day, + ht->totalFracRemovedAbove + ht->totalFracTransferredAbove, + ht->totalFracRemovedBelow + ht->totalFracTransferredBelow); + exit(EXIT_CODE_BAD_PARAMETER_VALUE); + } + } else { + logWarning("Harvest event at year %d day %d has no biomass to remove " + "or transfer\n", + 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.eventLitterC += litterAdd / climLen; + fluxes.eventSoilC += 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; + } break; + case TILLAGE: { + // BIG NOTE: this is the one event type that is NOT modeled as a flux; + // see updateEventTrackers() for more + const TillageParams *tillParams = event->eventParams; + // Update the tillage mod for R_H calculations; this will be slowly + // reduced by an exponential decay function. Note we add here, not set, + // as there may be lingering effects from a prior tillage. + eventTrackers.d_till_mod += tillParams->tillageEffect; + } break; + case FERTILIZATION: { + const FertilizationParams *fertParams = event->eventParams; + const double orgC = fertParams->orgC; + if (ctx.litterPool) { + fluxes.eventLitterC += orgC / climLen; + } else { + fluxes.eventSoilC += orgC / climLen; + } + + // MASS BALANCE: this is a system input + fluxes.eventInputC += orgC / climLen; + } break; + case LEAFON: { + double leafOnFlux = params.leafGrowth / climLen; + checkLeafOnLimitation(&leafOnFlux); + fluxes.eventLeafOnCreation += leafOnFlux; + double totalSourceC = envi.plantWoodC + envi.coarseRootC; + if (totalSourceC > TINY) { + fluxes.eventLeafOnCreationFromWood += + leafOnFlux * envi.plantWoodC / totalSourceC; + } + // Unlike planting, this is NOT a system input, so no adjustments to + // eventInputC + } break; + case LEAFOFF: { + double leafOff = envi.plantLeafC * params.fracLeafFall; + fluxes.eventLeafOffLitter += leafOff / climLen; + } break; + 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; + default: + logError("Unknown event type (%d) in processEvents()\n", event->type); + exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM); + } + + event = event->nextEvent; + } +} + +void processEventsForNitrogen(EventNode *event) { + if (!ctx.nitrogenCycle) { + return; + } + + // Event fluxes have all been reset to zero at the start of the time step, + // so we can just add to them as needed + + // If event starts off NULL, this function will just fall through, as it + // should. + const int climYear = climate->year; + const int climDay = climate->day; + const double climLen = climate->length; + + // Checks performed in processEventsForCarbon are not repeated here unless + // necessary for nitrogen handling. + + while (event != NULL && event->year <= climYear && event->day <= climDay) { + switch (event->type) { + case IRRIGATION: { + // Nothing to do for irrigation events + } break; + case PLANTING: { + const PlantingParams *plantParams = event->eventParams; + const double leafC = plantParams->leafC; + const double woodC = plantParams->woodC; + const double fineRootC = plantParams->fineRootC; + const double coarseRootC = plantParams->coarseRootC; + + // No need to allocate to biomass N pools, we don't track that N + // explicitly + + // MASS BALANCE: this is a system input + double inputN = leafC / params.leafCN + woodC / params.woodCN + + fineRootC / params.fineRootCN + + coarseRootC / params.woodCN; + fluxes.eventInputN += inputN / climLen; + } break; + case HARVEST: { + // Harvest can both remove biomass and move biomass to the soil/litter + // pools + const HarvestParams *harvParams = event->eventParams; + const double fracRA = harvParams->fractionRemovedAbove; + const double fracRB = harvParams->fractionRemovedBelow; + const double fracTA = harvParams->fractionTransferredAbove; + const double fracTB = harvParams->fractionTransferredBelow; + + // 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.eventLitterN += 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; + } break; + case TILLAGE: { + // Nothing to do for tillage events + } break; + case FERTILIZATION: { + const FertilizationParams *fertParams = event->eventParams; + double orgN = fertParams->orgN; + double minN = fertParams->minN; + // As the warning says in readEventData(), we ignore N when the + // nitrogen cycle model is off + // Implies ctx.litterPool + fluxes.eventLitterN += orgN / climLen; + fluxes.eventMinN += minN / climLen; + + // MASS BALANCE: this is a system input + fluxes.eventInputN += (orgN + minN) / climLen; + } 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 LEAFOFF: { + // Nitrogen - need to account for leaf N moving to litter, as with + // harvests + double leafOff = envi.plantLeafC * params.fracLeafFall; + double leafN = leafOff / params.leafCN; + double leafNResorption = leafN * params.leafNResorptionFrac; + double litterNAdd = leafN - leafNResorption; + fluxes.eventLeafOffNResorption += leafNResorption / climLen; + fluxes.eventLitterN += litterNAdd / climLen; + } break; + case PLANTDEATH: + // Nothing to do here + break; + default: + logError("Unknown event type (%d) in processEvents()\n", event->type); + exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM); + } + + event = event->nextEvent; + } +} + +void writeAllEvents(void) { // Event fluxes have all been reset to zero at the start of the time step, // so we can just add to them as needed diff --git a/src/sipnet/events.h b/src/sipnet/events.h index d829f6bb..de757fe2 100644 --- a/src/sipnet/events.h +++ b/src/sipnet/events.h @@ -81,6 +81,8 @@ struct EventNode { event_type_t type; int year, day; void *eventParams; + int numLogParamPairs; + char *logLine; EventNode *nextEvent; }; @@ -187,17 +189,30 @@ void setupEvents(void); int isFirstEventBefore(int year, int day); /*! - * \brief Process events for current location/year/day + * \brief Process carbon effects from events for current day * - * For a given year and day (as determined by the global `climate` - * pointer), process all events listed in the global `events` pointer for the - * referenced location. + * Process all events for the current day, calculating all carbon + * effects. For each event, modify flux variables according to the model for + * that event type. + */ +void processEventsForCarbon(EventNode *event); + +/*! + * \brief Process nitrogen effects from events for current day + * + * Process all events for the current day, calculating all nitrogen + * effects. For each event, modify flux variables according to the model for + * that event type. * - * For each event, modify flux variables according to the model for that event, - * and write a row to the configured event output file listing the modified - * variables and the delta applied. + * Carbon and nitrogen effects are calculated separately to allow carbon + * limitation checks to be run before any nitrogen calculations are made. + */ +void processEventsForNitrogen(EventNode *event); + +/*! + * \brief Write events to the events output file for the current date */ -void processEvents(void); +void writeAllEvents(void); /*! * Update relevant environment pools after event fluxes have been calculated diff --git a/src/sipnet/limitations.c b/src/sipnet/limitations.c index 726cedc2..1473b3dc 100644 --- a/src/sipnet/limitations.c +++ b/src/sipnet/limitations.c @@ -150,18 +150,55 @@ static void checkNegativeCreation(void) { // appropriately. double len = climate->length; + + logInfo("BEFORE: leafLitter %f eventLOLitter %f leafC %f leafCreation %f\n", + fluxes.leafLitter * len, fluxes.eventLeafOffLitter * len, + envi.plantLeafC, fluxes.leafCreation * len); + // Above ground // If leafCreation is too negative, we need to deduct from wood instead - // Use only the continuous turnover term to match previous logic - but see - // SIPNET issue #372. - double leafLitterTurnover = envi.plantLeafC * params.leafTurnoverRate; - double leafDeficit = - envi.plantLeafC / len + fluxes.leafCreation - leafLitterTurnover; - if (leafDeficit < 0) { - fluxes.woodCreation += leafDeficit; - fluxes.leafCreation -= leafDeficit; + // Need to make sure we don't go too far the other way. Leaf off litter + // (either fluxes.leafLitter or fluxes.eventLeafOffLitter) might also + // incorrectly drive the pool negative, sp adjust for that too. + + // First we handle leaf litter. Assuming params.leafTurnoverRate is valid + // (ie, <=1), availableLeafRate should be non-negative. + double availableLeafRate = envi.plantLeafC / len - fluxes.leafLitter; + + // Reminder: litter and creation fluxes "point" in the opposite direction, so + // their signs are opposite here + double deficit = fluxes.leafOffLitter + fluxes.eventLeafOffLitter - + fluxes.leafCreation - availableLeafRate; + if (deficit > 0) { + // If negative growth + leaf off is too much, let's first reduce leaf off + if (fluxes.eventLeafOffLitter > 0) { + double adjust = fmin(fluxes.eventLeafOffLitter, deficit); + fluxes.eventLeafOffLitter -= adjust; + deficit -= adjust; + // If we change EVENT leaf off litter, we need to adjust eventLitterN, as + // it has already been calculated + fluxes.eventLitterN + } + if (deficit > 0) { + if (fluxes.leafOffLitter > 0) { + double adjust = fmin(fluxes.leafOffLitter, deficit); + fluxes.leafOffLitter -= adjust; + deficit -= adjust; + } + } + // Next, move some negative growth to the wood pool + if (deficit > 0) { + // If this is too much for the wood pool to handle, we have an error that + // will be caught in ensureNonNegative + fluxes.woodCreation -= deficit; + fluxes.leafCreation += deficit; + } } + logInfo("AFTER: leafLitter %f eventLOLitter %f leafC %f leafCreation %f\n", + fluxes.leafLitter * len, fluxes.eventLeafOffLitter * len, + envi.plantLeafC, fluxes.leafCreation * len); + // Below ground double fineRootDeficit = envi.fineRootC / len + fluxes.fineRootCreation - fluxes.fineRootLoss; diff --git a/src/sipnet/nitrogen.c b/src/sipnet/nitrogen.c index 739eb290..8c9ffbf1 100644 --- a/src/sipnet/nitrogen.c +++ b/src/sipnet/nitrogen.c @@ -68,7 +68,7 @@ static void calcNPoolFluxes(void) { // litter (modified by leaf N resorption), and N loss due to mineralization. // N added via fertilization is handled elsewhere. fluxes.nOrgLitter = - fluxes.leafLitter / params.leafCN - fluxes.leafOffNResorption + + getLeafLitterFlux() / params.leafCN - fluxes.leafOffNResorption + fluxes.woodLitter / params.woodCN - litterMin - fluxes.litterToSoil / litterCN + (soilNInputs * saturationFraction); @@ -186,11 +186,10 @@ void calcNResorptionFluxes(void) { fluxes.fineRootCreation / params.fineRootCN); } - // Leaf litter resorption; at this point, fluxes.leafLitter counts both normal - // turnover and leaf-off calcs. Note that event leaf off is handled in - // events.c + // Leaf litter resorption; at this point. Note that event leaf off is handled + // in events.c double nResorp = - params.leafNResorptionFrac * fluxes.leafLitter / params.leafCN; + params.leafNResorptionFrac * getLeafLitterFlux() / params.leafCN; fluxes.leafOffNResorption += nResorp; // TODO: Should we resorb N from wood litter? diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index e54aaacb..5014e4e3 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -791,17 +791,8 @@ void calcWoodAndLeafFluxes(void) { * the start or end of the growing season. These transition fluxes are tracked * separately from the continuous leaf creation and litter fluxes calculated in * calcLeafFluxes(). - * - * @param[out] leafOnCreation Additional leaf creation flux at leaf-on - * (g C/m^2 ground/day) - * @param[out] leafOnFromWood Carbon transferred from wood to leaves at - * leaf-on (g C/m^2 ground/day) - * @param[out] leafLitter Additional leaf litter flux at leaf-off - * (g C/m^2 ground/day) - * @param[in] plantLeafC Leaf carbon pool size (g C/m^2 ground area) */ -void calcLeafOnOffFluxes(double *leafOnCreation, double *leafOnFromWood, - double *leafLitter, double plantLeafC) { +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 @@ -818,10 +809,10 @@ void calcLeafOnOffFluxes(double *leafOnCreation, double *leafOnFromWood, // we just reached the start of the growing season double leafOn = params.leafGrowth / climate->length; checkLeafOnLimitation(&leafOn); - *leafOnCreation += leafOn; + fluxes.leafOnCreation += leafOn; double totalSourceC = envi.plantWoodC + envi.coarseRootC; if (totalSourceC > TINY) { - *leafOnFromWood += leafOn * envi.plantWoodC / totalSourceC; + fluxes.leafOnCreationFromWood += leafOn * envi.plantWoodC / totalSourceC; } phenologyTrackers.didLeafGrowth = 1; // This is a computed event - however, the value may get reduced by @@ -833,8 +824,8 @@ void calcLeafOnOffFluxes(double *leafOnCreation, double *leafOnFromWood, if (!phenologyTrackers.didLeafFall && pastLeafFall()) { // we just reached the end of the growing season double len = climate->length; - double leafOff = (plantLeafC * params.fracLeafFall) / len; - *leafLitter += leafOff; + double leafOff = (envi.plantLeafC * params.fracLeafFall) / len; + fluxes.leafOffLitter += leafOff; phenologyTrackers.didLeafFall = 1; if (leafOff > TINY && ctx.events) { writeComputedEventOut(climate->year, climate->day, @@ -1301,8 +1292,7 @@ void calculateFluxes(void) { calcWoodAndLeafFluxes(); // Leaf on/off - calcLeafOnOffFluxes(&fluxes.leafOnCreation, &fluxes.leafOnCreationFromWood, - &fluxes.leafLitter, envi.plantLeafC); + calcLeafOnOffFluxes(); // Litter pool, if LITTER is on calcLitterFluxes(); @@ -1478,7 +1468,7 @@ void updateTrackers(double oldSoilWater) { // If we get another event flux in this function, we should create an // updateTrackersForEvents() function in events.c|h - trackers.yearlyLitter += fluxes.leafLitter + fluxes.eventLeafOffLitter; + trackers.yearlyLitter += getLeafLitterFlux() + fluxes.eventLeafOffLitter; if (ctx.gdd) { trackers.gdd += climate->gdd; @@ -1616,7 +1606,7 @@ void updateMainPools(void) { // L_L = fluxes.leafLitter // Note: we have split leafCreation into two parts envi.plantLeafC += - (fluxes.leafCreation + fluxes.leafOnCreation - fluxes.leafLitter) * + (fluxes.leafCreation + fluxes.leafOnCreation - getLeafLitterFlux()) * climate->length; // :: from [1], eq (A4), where: @@ -1657,7 +1647,7 @@ void updatePoolsForSoil(void) { ctx.carbonSaturation ? unitClip(envi.soilC / params.soilCSaturation) : 0.0; // :: from [2], litter model description - envi.litterC += (fluxes.woodLitter + fluxes.leafLitter + + envi.litterC += (fluxes.woodLitter + getLeafLitterFlux() + (soilInputs * saturationFraction) - fluxes.litterToSoil - fluxes.rLitter - fluxes.litterMethane) * climate->length; @@ -1669,12 +1659,12 @@ void updatePoolsForSoil(void) { // Normal pool (single pool, no microbes) // :: from [1] (and others, TBD), eq (A3), where: // L_w = fluxes.woodLitter - // L_l = fluxes.leafLitter + // L_l = fluxes.leafLitter + fluxes.leafOffLitter // R_h = fluxes.rSoil // :: from [3], root terms envi.soilC += (fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.woodLitter + - fluxes.leafLitter - fluxes.rSoil - fluxes.soilMethane) * + getLeafLitterFlux() - fluxes.rSoil - fluxes.soilMethane) * climate->length; } diff --git a/src/sipnet/state.c b/src/sipnet/state.c index fdaa4391..b082cbcd 100644 --- a/src/sipnet/state.c +++ b/src/sipnet/state.c @@ -17,3 +17,7 @@ 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 f77bc791..d341eade 100644 --- a/src/sipnet/state.h +++ b/src/sipnet/state.h @@ -473,7 +473,7 @@ typedef struct FluxVars { // GROSS photosynthesis (g C * m^-2 ground area * day^-1) double photosynthesis; - // Leaf fall (g C * m^-2 ground area * day^-1) + // Leaf fall (g C * m^-2 ground area * day^-1); excludes leaf-off events double leafLitter; // Wood flux to litter (g C * m^-2 ground area * day^-1) double woodLitter; @@ -550,6 +550,9 @@ typedef struct FluxVars { // Portion of leaf-on creation C that comes from wood C (the rest comes from // coarse root C) double leafOnCreationFromWood; + // leaf fall from calculated leaf off 'events' (gdd, soil temp, day of year) + // (g C * m^-2 grount area * day^-1) + double leafOffLitter; // **************************************** // Fluxes for nitrogen cycle @@ -757,4 +760,6 @@ extern PlantSurvivalTracker plantSurvivalTracker; double getTotalWoodC(void); +double getLeafLitterFlux(void); + #endif // SIPNET_STATE_H diff --git a/tests/sipnet/test_modeling/testCompleteHarvest.c b/tests/sipnet/test_modeling/testCompleteHarvest.c index e0e22e35..d2d48ed6 100644 --- a/tests/sipnet/test_modeling/testCompleteHarvest.c +++ b/tests/sipnet/test_modeling/testCompleteHarvest.c @@ -70,8 +70,8 @@ static void fullCase(int mode, int accounting, int dark, int limited, const char *harvest, double aboveExport, double belowExport) { // A no-harvest run supplies the independently calculated end-of-step pools. - logTest("Running full case with mode %d accounting %d dark %d limited %d" - " harvest %s aboveExport %.2f belowExport %.2f\n", + logTest("Running full case with mode: %d accounting: %d dark: %d limited: %d" + " harvest: %s aboveExport: %.2f belowExport: %.2f\n", mode, accounting, dark, limited, harvest ? harvest : "none", aboveExport, belowExport); Envi end = {0}; @@ -152,6 +152,8 @@ static void fullCase(int mode, int accounting, int dark, int limited, } static void partialCase(int accounting, const char *harvest, double fraction) { + logTest("Running partial case with accounting: %d harvest: %s fraction: %f\n", + accounting, harvest, fraction); start(2, harvest); envi.plantCAccountingDelta = accounting; Envi before = envi; @@ -172,6 +174,7 @@ static void partialCase(int accounting, const char *harvest, double fraction) { } static void partialTimestepCase(int accounting) { + logTest("Running partial timestep case\n"); start(2, ".1 .2 .3 .4"); envi.plantCAccountingDelta = accounting; updateState(); @@ -219,6 +222,7 @@ static void invalidCase(int which) { } static void coincidentFertilizerCase(void) { + logTest("Running coincident fertilizer case\n"); Envi expected = {0}; for (int h = 0; h < 2; h++) { start(2, h ? "0 0 1 1" : NULL); @@ -263,6 +267,8 @@ static void sameEnvi(const Envi *a, const Envi *b) { #undef E } static void coincidentEventCase(int type, const char *arguments) { + logTest("Running coincident event case with type: %s params: %s\n", + eventTypeToString(type), *arguments ? arguments : ""); Envi end; Fluxes ordinary; for (int order = 0; order < 3; order++) { @@ -329,6 +335,7 @@ static void coincidentEventCase(int type, const char *arguments) { } static void restartCase(void) { + logTest("Running restart case\n"); Envi expected; for (int resumed = 0; resumed < 2; resumed++) { start(2, "0 0 1 1"); @@ -369,13 +376,18 @@ static void restartCase(void) { } static void leafBudgetCases(void) { - for (int kind = 0; kind < 4; kind++) - for (int dark = 0; dark < 2; dark++) - for (int sign = -1; sign <= 1; sign++) - for (int account = -1; account <= 1; account++) - for (int limited = 0; limited < 2; limited++) - for (int resorb = 0; resorb < 3; resorb++) + for (int kind = 0; kind < 4; kind++) { + for (int dark = 0; dark < 2; dark++) { + for (int sign = -1; sign <= 1; sign++) { + for (int account = -1; account <= 1; account++) { + for (int limited = 0; limited < 2; limited++) { + for (int resorb = 0; resorb < 3; resorb++) { for (int harvest = 0; harvest < 2; harvest++) { + logTest("*** Running leaf budget case with kind: %d dark: %d " + "sign %d" + " account %d limited %d resorb %d harvest %s\n", + kind, dark, sign, account, limited, resorb, + harvest ? "0 0 1 1" : "none"); start(2, harvest ? "0 0 1 1" : NULL); envi.plantLeafC = 1; envi.plantWoodC = envi.fineRootC = envi.coarseRootC = 100; @@ -411,7 +423,8 @@ static void leafBudgetCases(void) { double beforeLitterN = envi.litterN; updateState(); balanced(); - double shed = (fluxes.leafLitter + fluxes.eventLeafOffLitter) * + double shed = (fluxes.leafLitter + fluxes.leafOffLitter + + fluxes.eventLeafOffLitter) * climate->length; if (shed < -1e-10 || shed > 1 + 1e-10 || envi.plantLeafC < -1e-10) { @@ -419,7 +432,9 @@ static void leafBudgetCases(void) { envi.plantLeafC); failures++; } - double eventExpected = kind == 3 ? 0 : kind == 1 ? .25 : 1; + + double expEventLOL[] = {1, .25, 1, 0}; + double eventExpected = expEventLOL[kind]; near(fluxes.eventLeafOffLitter * climate->length, eventExpected, "event leaves conserved"); if (!harvest) { @@ -438,7 +453,18 @@ static void leafBudgetCases(void) { "no regrowth after complete harvest"); } finish(); - } + } // harvest loop + } // resorb loop + + logTest("Complete harvest failures: %d\n", failures); + logTest("Press return to continue\n"); + getchar(); + + } // limited loop + } // account loop + } // sign loop + } // dark loop + } // kind loop } // Full export during net loss must not borrow carbon from an empty soil pool. @@ -475,19 +501,24 @@ int main(void) { partialCase(account, "0 0 .999999999 .999999999", .999999999); } + logTest("Running invalid case checks; six errors expected\n"); for (int which = 0; which < 6; which++) { invalidCase(which); } - logTest("Press return to continue\n"); - getchar(); - coincidentFertilizerCase(); coincidentEventCase(LEAFON, ""); coincidentEventCase(LEAFOFF, ""); coincidentEventCase(IRRIGATION, "2 0"); restartCase(); + + logTest("\n\n\n"); + leafBudgetCases(); + + logTest("Press return to continue\n"); + getchar(); + exportLossCase(); logTest("Complete harvest failures: %d\n", failures); return failures != 0; From 2e274317879aa60d0e8fc8d6dbf7a18d5b553b34 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Mon, 28 Sep 2026 16:40:45 -0400 Subject: [PATCH 04/14] Interim commit, through leaf case 251 --- src/common/exitCodes.h | 3 +- src/common/util.c | 106 ++++ src/common/util.h | 22 + src/sipnet/debug_log.c | 2 +- src/sipnet/events.c | 508 +++++------------- src/sipnet/events.h | 41 +- src/sipnet/limitations.c | 219 +++++--- src/sipnet/nitrogen.c | 24 +- src/sipnet/nitrogen.h | 9 + src/sipnet/sipnet.c | 64 +-- src/sipnet/state.h | 4 +- .../test_modeling/testCompleteHarvest.c | 84 ++- 12 files changed, 550 insertions(+), 536 deletions(-) diff --git a/src/common/exitCodes.h b/src/common/exitCodes.h index 1f17776b..2c399490 100644 --- a/src/common/exitCodes.h +++ b/src/common/exitCodes.h @@ -23,7 +23,8 @@ typedef enum { EXIT_CODE_FILE_OPEN_OR_READ_ERROR = 6, EXIT_CODE_INTERNAL_ERROR = 7, EXIT_CODE_BAD_CLI_ARGUMENT = 8, - EXIT_CODE_BAD_RESTART_PARAMETER = 9 + EXIT_CODE_BAD_RESTART_PARAMETER = 9, + EXIT_CODE_MEMORY_ALLOCATION_FAILURE = 10 } exit_code_t; #endif diff --git a/src/common/util.c b/src/common/util.c index fb087421..632c8e75 100644 --- a/src/common/util.c +++ b/src/common/util.c @@ -76,3 +76,109 @@ double calcRatio(const double num, const double den) { // For global linkage extern inline double unitClip(double preClip); + +DynamicString *dsCreate(size_t initial_capacity) { + DynamicString *ds = malloc(sizeof(DynamicString)); + if (!ds) + return NULL; + + // Ensure we have room for at least a null terminator + ds->capacity = (initial_capacity > 0) ? initial_capacity : 16; + ds->buffer = malloc(ds->capacity * sizeof(char)); + + if (!ds->buffer) { + free(ds); + return NULL; + } + + ds->buffer[0] = '\0'; // Start with an empty string + ds->length = 0; + return ds; +} + +int dsAppend(DynamicString *ds, const char *str) { + if (!ds || !str) + return 0; + + size_t append_len = strlen(str); + // +1 is crucial to ensure room for the null terminator + size_t needed_capacity = ds->length + append_len + 1; + + // Double capacity until it's large enough for the new content + if (needed_capacity > ds->capacity) { + size_t new_capacity = ds->capacity * 2; + while (new_capacity < needed_capacity) { + new_capacity *= 2; + } + + // Safely reallocate using a temporary pointer + char *temp = realloc(ds->buffer, new_capacity * sizeof(char)); + if (!temp) { + return 0; // Reallocation failed, original data remains intact + } + + ds->buffer = temp; + ds->capacity = new_capacity; + } + + // Copy the new string over the old null terminator + memcpy(ds->buffer + ds->length, str, append_len); + ds->length += append_len; + ds->buffer[ds->length] = '\0'; // Manually place the new null terminator + + return 1; // Success +} + +int dsAppendFormatted(DynamicString *ds, const char *format, ...) { + if (!ds || !format) + return 0; + + // 1. Determine how much space the formatted text needs + va_list args; + va_start(args, format); + // Make a copy of args because vsnprintf consumes the list + va_list args_copy; + va_copy(args_copy, args); + + // Pass NULL and 0 to just count the characters needed + int formatted_len = vsnprintf(NULL, 0, format, args_copy); + va_end(args_copy); + + if (formatted_len < 0) { + va_end(args); + return 0; // Formatting error occurred + } + + // 2. Ensure the buffer is big enough + size_t needed_capacity = ds->length + (size_t)formatted_len + 1; + if (needed_capacity > ds->capacity) { + size_t new_capacity = ds->capacity * 2; + while (new_capacity < needed_capacity) { + new_capacity *= 2; + } + + char *temp = realloc(ds->buffer, new_capacity * sizeof(char)); + if (!temp) { + va_end(args); + return 0; // Reallocation failed + } + ds->buffer = temp; + ds->capacity = new_capacity; + } + + // 3. Write the formatted string directly into the builder's buffer + // Write directly to the position of the current null terminator + vsnprintf(ds->buffer + ds->length, (size_t)formatted_len + 1, format, args); + va_end(args); + + // 4. Update the length of the string builder + ds->length += (size_t)formatted_len; + return 1; +} + +void dsFree(DynamicString *ds) { + if (ds) { + free(ds->buffer); + free(ds); + } +} diff --git a/src/common/util.h b/src/common/util.h index c8a31f26..8344d7a3 100644 --- a/src/common/util.h +++ b/src/common/util.h @@ -9,7 +9,10 @@ #define UTIL_H #include +#include // Required for va_list, va_start, va_end #include +#include +#include #define TINY 0.000001 // to avoid those nasty divide-by-zero errors @@ -37,4 +40,23 @@ double calcRatio(double num, double den); */ inline double unitClip(double preClip) { return fmin(fmax(preClip, 0.0), 1.0); } +// DYNAMIC STRING +typedef struct DynamicStringStruct { + char *buffer; // Pointer to the character array + size_t length; // Number of characters currently in the string + size_t capacity; // Total allocated space (including null terminator) +} DynamicString; + +// Initialize the builder with a reasonable default capacity +DynamicString *dsCreate(size_t initial_capacity); + +// Append text, automatically doubling capacity if needed +int dsAppend(DynamicString *ds, const char *str); + +// Append formatted text using printf-style syntax +int dsAppendFormatted(DynamicString *ds, const char *format, ...); + +// Free all memory associated with the builder +void dsFree(DynamicString *ds); + #endif diff --git a/src/sipnet/debug_log.c b/src/sipnet/debug_log.c index 483ff3fe..1c62d77a 100644 --- a/src/sipnet/debug_log.c +++ b/src/sipnet/debug_log.c @@ -117,7 +117,7 @@ void initDebugArrays() { debugFields->fluxDF[ind++] = (DebugField){"eventOutputN", DEBUG_FIELD_DOUBLE, &fluxes.eventOutputN}, debugFields->fluxDF[ind++] = (DebugField){"eventLeafOnCreation", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOnCreation}, debugFields->fluxDF[ind++] = (DebugField){"eventLeafOnCreationFromWood", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOnCreationFromWood}, - debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitter", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitter}, + debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitter", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitterC}, debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffNResorption", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffNResorption}, debugFields->fluxDF[ind++] = (DebugField){"soilMethane", DEBUG_FIELD_DOUBLE, &fluxes.soilMethane}, debugFields->fluxDF[ind ] = (DebugField){"litterMethane", DEBUG_FIELD_DOUBLE, &fluxes.litterMethane}; diff --git a/src/sipnet/events.c b/src/sipnet/events.c index 66401ede..0edc556a 100644 --- a/src/sipnet/events.c +++ b/src/sipnet/events.c @@ -44,7 +44,7 @@ EventNode *createEventNode(int year, int day, int eventType, newEvent->day = day; newEvent->type = eventType; newEvent->numLogParamPairs = 0; - newEvent->logLine = NULL; + newEvent->logLine = dsCreate(0); switch (eventType) { case HARVEST: { @@ -377,51 +377,57 @@ EventNode *readEventData(const char *eventFile) { void openEventOutFile(const char *eventOutFilePath, int printHeader) { eventOutFile = openFile(eventOutFilePath, "w"); if (printHeader) { - // Use format string analogous to the one in writeEventOut for + // Use format string analogous to the one in writeEventsOut for // better alignment (won't be perfect, but definitely better) fprintf(eventOutFile, "%4s %3s %-7s %s", "year", "day", "type", "param_name=delta[,param_name=delta,...]\n"); } } -void doWriteEventOut(int year, int day, const char *type, int numParams, - va_list args) { - int ind = 0; - - // Spec: - // year day event_type [,=,...] - - // 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++) { - // For debugging on linux - param = va_arg(args, char *); - val = va_arg(args, double); - fprintf(eventOutFile, "%s=%-.2f,", param, val); - } - param = va_arg(args, char *); - val = va_arg(args, double); - fprintf(eventOutFile, "%s=%-.2f\n", param, val); -} - -void writeEventOut(EventNode *oneEvent, int numParams, ...) { +void appendLog(EventNode *event, int numParams, ...) { va_list args; va_start(args, numParams); - doWriteEventOut(oneEvent->year, oneEvent->day, - eventTypeToString(oneEvent->type), numParams, args); + for (int ind = 0; ind < numParams; ind++) { + char *param = va_arg(args, char *); + double val = va_arg(args, double); + char suffix = (ind == numParams - 1) ? '\n' : ','; + int success = + dsAppendFormatted(event->logLine, "%s=%-.2f%c", param, val, suffix); + if (!success) { + logError("appending event log line failed\n"); + exit(EXIT_CODE_INTERNAL_ERROR); + } + } va_end(args); } +void writeEventsOut(void) { + // We move the global events pointer here (instead of a copy as in the + // processEvents functions), as this is the last pass + const int climYear = climate->year; + const int climDay = climate->day; + while (gEvent != NULL && gEvent->year <= climYear && gEvent->day <= climDay) { + fprintf(eventOutFile, "%4d %3d %-7s %s", gEvent->year, gEvent->day, + eventTypeToString(gEvent->type), gEvent->logLine->buffer); + gEvent = gEvent->nextEvent; + } +} + void writeComputedEventOut(int year, int day, const char *type, int numParams, ...) { va_list args; va_start(args, numParams); - doWriteEventOut(year, day, type, numParams, args); + // Standard prefix for all + fprintf(eventOutFile, "%4d %3d %-7s ", year, day, type); + + // Variable output per oneEvent type + for (int ind = 0; ind < numParams; ind++) { + char *param = va_arg(args, char *); + double val = va_arg(args, double); + char suffix = (ind == numParams - 1) ? '\n' : ','; + fprintf(eventOutFile, "%s=%-.2f%c", param, val, suffix); + } + va_end(args); } @@ -454,6 +460,8 @@ int isFirstEventBefore(int year, int day) { return firstEvent->day < day; } +EventNode *getCurrentEvent(void) { return gEvent; } + void processEventsForCarbon(EventNode *event) { // Event fluxes have all been reset to zero at the start of the time step, // so we can just add to them as needed @@ -508,6 +516,10 @@ void processEventsForCarbon(EventNode *event) { } fluxes.eventEvap += evapAmount / climLen; fluxes.eventSoilWater += soilAmount / climLen; + + appendLog(gEvent, 2, "eventSoilWater", soilAmount, "eventEvap", + evapAmount); + } break; case PLANTING: { const PlantingParams *plantParams = event->eventParams; @@ -528,6 +540,16 @@ void processEventsForCarbon(EventNode *event) { // MASS BALANCE: this is a system input const double inputC = leafC + woodC + fineRootC + coarseRootC; fluxes.eventInputC += inputC / climLen; + + // clang-format off + appendLog(gEvent, 5, + "eventLeafC", leafC, + "eventWoodC", woodC, + "eventFineRootC", fineRootC, + "eventCoarseRootC", coarseRootC, + "eventInputC", inputC); + // clang-format on + } break; case HARVEST: { // Harvest can both remove biomass and move biomass to the soil/litter @@ -603,6 +625,20 @@ void processEventsForCarbon(EventNode *event) { const double outputC = ((woodC + envi.plantLeafC) * fracRA + (envi.fineRootC + envi.coarseRootC) * fracRB); fluxes.eventOutputC += outputC / climLen; + + // clang-format off + appendLog( + gEvent, 8, + "eventSoilC", soilAdd, + "eventLitterC", litterAdd, + "eventLeafC", leafDelta, + "eventWoodC", woodDelta, + "eventAccountingC", accountingDelta, + "eventFineRootC", fineDelta, + "eventCoarseRootC", coarseDelta, + "eventOutputC", outputC); + // clang-format on + } break; case TILLAGE: { // BIG NOTE: this is the one event type that is NOT modeled as a flux; @@ -612,6 +648,10 @@ void processEventsForCarbon(EventNode *event) { // reduced by an exponential decay function. Note we add here, not set, // as there may be lingering effects from a prior tillage. eventTrackers.d_till_mod += tillParams->tillageEffect; + + appendLog(gEvent, 1, "eventTrackers.d_till_mod", + tillParams->tillageEffect); + } break; case FERTILIZATION: { const FertilizationParams *fertParams = event->eventParams; @@ -624,22 +664,39 @@ void processEventsForCarbon(EventNode *event) { // MASS BALANCE: this is a system input fluxes.eventInputC += orgC / climLen; + + // clang-format off + appendLog(gEvent, 3, + "eventLitterC", ctx.litterPool ? orgC : 0.0, + "eventSoilC", ctx.litterPool ? 0.0 : orgC, + "eventInputC", orgC); + // clang-format on + } break; case LEAFON: { double leafOnFlux = params.leafGrowth / climLen; + double leafOnFluxFromWood = 0.0; checkLeafOnLimitation(&leafOnFlux); fluxes.eventLeafOnCreation += leafOnFlux; double totalSourceC = envi.plantWoodC + envi.coarseRootC; if (totalSourceC > TINY) { - fluxes.eventLeafOnCreationFromWood += - leafOnFlux * envi.plantWoodC / totalSourceC; + leafOnFluxFromWood = leafOnFlux * envi.plantWoodC / totalSourceC; + fluxes.eventLeafOnCreationFromWood += leafOnFluxFromWood; } // Unlike planting, this is NOT a system input, so no adjustments to // eventInputC + + // clang-format off + appendLog(gEvent, 2, + "eventLeafOnCreation", leafOnFlux * climLen, + "eventLeafOnCreationFromWood", leafOnFluxFromWood * climLen); + // clang-format on } break; case LEAFOFF: { double leafOff = envi.plantLeafC * params.fracLeafFall; - fluxes.eventLeafOffLitter += leafOff / climLen; + fluxes.eventLeafOffLitterC += leafOff / climLen; + + appendLog(gEvent, 1, "eventLeafOffLitter", leafOff); } break; case PLANTDEATH: // There should be no way to get here, but covering our bases... @@ -673,6 +730,10 @@ void processEventsForNitrogen(EventNode *event) { // Checks performed in processEventsForCarbon are not repeated here unless // necessary for nitrogen handling. + // Leaf-off events are a little special, as we only process the first one + // (see LEAFOFF below) + int firstLeafOff = 1; + while (event != NULL && event->year <= climYear && event->day <= climDay) { switch (event->type) { case IRRIGATION: { @@ -693,6 +754,8 @@ void processEventsForNitrogen(EventNode *event) { fineRootC / params.fineRootCN + coarseRootC / params.woodCN; fluxes.eventInputN += inputN / climLen; + + appendLog(gEvent, 1, "eventInputN", inputN); } break; case HARVEST: { // Harvest can both remove biomass and move biomass to the soil/litter @@ -725,6 +788,14 @@ void processEventsForNitrogen(EventNode *event) { envi.coarseRootC / params.woodCN) * fracRB; fluxes.eventOutputN += outputN / climLen; + + // clang-format off + appendLog( + gEvent, 3, + "eventSoilOrgN", soilNAdd, + "eventLitterN", litterNAdd, + "eventOutputN", outputN); + // clang-format on } break; case TILLAGE: { // Nothing to do for tillage events @@ -741,353 +812,63 @@ void processEventsForNitrogen(EventNode *event) { // MASS BALANCE: this is a system input fluxes.eventInputN += (orgN + minN) / climLen; - } 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 LEAFOFF: { - // Nitrogen - need to account for leaf N moving to litter, as with - // harvests - double leafOff = envi.plantLeafC * params.fracLeafFall; - double leafN = leafOff / params.leafCN; - double leafNResorption = leafN * params.leafNResorptionFrac; - double litterNAdd = leafN - leafNResorption; - fluxes.eventLeafOffNResorption += leafNResorption / climLen; - fluxes.eventLitterN += litterNAdd / climLen; - } break; - case PLANTDEATH: - // Nothing to do here - break; - default: - logError("Unknown event type (%d) in processEvents()\n", event->type); - exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM); - } - - event = event->nextEvent; - } -} - -void writeAllEvents(void) { - // Event fluxes have all been reset to zero at the start of the time step, - // so we can just add to them as needed - - // If event starts off NULL, this function will just fall through, as it - // should. - const int climYear = climate->year; - const int climDay = climate->day; - const double climLen = climate->length; - - // As this is used as a divisor in many places, let's make sure it's >0 - if (climLen <= 0) { - logError("climate length (%f) on year %d day %d is non-positive; please " - "fix and re-run", - climLen, climYear, climDay); - exit(EXIT_CODE_BAD_PARAMETER_VALUE); - } - - // Reset harvest tracking - eventTrackers.harvestTrackers = (HarvestTrackers){0}; - - while (gEvent != NULL && gEvent->year <= climYear && gEvent->day <= climDay) { - // The events file has been tested on read, so we know this event list - // should be in chrono order. However, we need to check to make sure the - // current event is not in the past, as that would indicate an event that - // did not have a corresponding climate file record. - if (gEvent->year < climYear || gEvent->day < climDay) { - logError("Agronomic event found for year: %d day: %d that does not " - "have a corresponding record in the climate file\n", - gEvent->year, gEvent->day); - exit(EXIT_CODE_INPUT_FILE_ERROR); - } - - switch (gEvent->type) { - case IRRIGATION: { - const IrrigationParams *irrParams = gEvent->eventParams; - const double amount = irrParams->amountAdded; - double soilAmount, evapAmount; - if (irrParams->method == CANOPY) { - // Part of the irrigation evaporates, and the rest makes it to the - // soil. Evaporated fraction: - evapAmount = params.immedEvapFrac * amount; - // Soil fraction: - soilAmount = amount - evapAmount; - } else if (irrParams->method == SOIL) { - // All goes to the soil - evapAmount = 0.0; - soilAmount = amount; - } else { - logError("Unknown irrigation method type: %d\n", irrParams->method); - exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM); - } - fluxes.eventEvap += evapAmount / climLen; - fluxes.eventSoilWater += soilAmount / climLen; - writeEventOut(gEvent, 2, "eventSoilWater", soilAmount, "eventEvap", - evapAmount); - } break; - case PLANTING: { - const PlantingParams *plantParams = gEvent->eventParams; - const double leafC = plantParams->leafC; - const double woodC = plantParams->woodC; - const double fineRootC = plantParams->fineRootC; - const double coarseRootC = plantParams->coarseRootC; - - // Update the fluxes - fluxes.eventLeafC += leafC / climLen; - fluxes.eventWoodC += woodC / climLen; - fluxes.eventFineRootC += fineRootC / climLen; - fluxes.eventCoarseRootC += coarseRootC / climLen; - - // No need to allocate to biomass N pools, we don't track that N - // explicitly - - // MASS BALANCE: this is a system input - const double inputC = leafC + woodC + fineRootC + coarseRootC; - double inputN = 0.0; - fluxes.eventInputC += inputC / climLen; - if (ctx.nitrogenCycle) { - inputN = leafC / params.leafCN + woodC / params.woodCN + - fineRootC / params.fineRootCN + coarseRootC / params.woodCN; - fluxes.eventInputN += inputN / climLen; - } // clang-format off - writeEventOut(gEvent, 6, - "eventLeafC", leafC, - "eventWoodC", woodC, - "eventFineRootC", fineRootC, - "eventCoarseRootC", coarseRootC, - "eventInputC", inputC, - "eventInputN", inputN); - // clang-format on - } break; - case HARVEST: { - // 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 fracRB = harvParams->fractionRemovedBelow; - const double fracTA = harvParams->fractionTransferredAbove; - const double fracTB = harvParams->fractionTransferredBelow; - const double woodC = envi.plantWoodC + envi.plantCAccountingDelta; - - // Record fraction of total biomass removed and transferred - double aboveMass = woodC + envi.plantLeafC; - double belowMass = envi.fineRootC + envi.coarseRootC; - double totalMass = aboveMass + belowMass; - HarvestTrackers *ht = &eventTrackers.harvestTrackers; - if (totalMass > TINY) { - double massRemoved = fracRA * aboveMass + fracRB * belowMass; - double massTransferred = fracTA * aboveMass + fracTB * belowMass; - ht->totalFracRemoved += massRemoved / totalMass; - ht->totalFracTransferred += massTransferred / totalMass; - ht->totalFracRemovedAbove += fracRA; - ht->totalFracRemovedBelow += fracRB; - ht->totalFracTransferredAbove += fracTA; - ht->totalFracTransferredBelow += fracTB; - if (ht->totalFracRemovedAbove + ht->totalFracTransferredAbove > - 1.0 + TINY || - ht->totalFracRemovedBelow + ht->totalFracTransferredBelow > - 1.0 + TINY) { - logError("Harvest event(s) at year %d day %d has total above-ground" - " or below-ground removal + transfer fraction > 1.0" - " (above %.3f, below %.3f)\n", - gEvent->year, gEvent->day, - ht->totalFracRemovedAbove + ht->totalFracTransferredAbove, - ht->totalFracRemovedBelow + ht->totalFracTransferredBelow); - exit(EXIT_CODE_BAD_PARAMETER_VALUE); - } - } else { - logWarning("Harvest event at year %d day %d has no biomass to remove " - "or transfer\n", - gEvent->year, gEvent->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.eventLitterC += litterAdd / climLen; - fluxes.eventSoilC += soilAdd / climLen; - fluxes.eventLeafC += leafDelta / climLen; - fluxes.eventWoodC += woodDelta / climLen; - fluxes.eventAccountingC += accountingDelta / 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 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; - 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); - double outputN = 0.0; - fluxes.eventOutputC += outputC / climLen; - if (ctx.nitrogenCycle) { - // just plantWoodC here, not woodC - outputN = (envi.plantWoodC / params.woodCN + - envi.plantLeafC / params.leafCN) * - fracRA + - (envi.fineRootC / params.fineRootCN + - envi.coarseRootC / params.woodCN) * - fracRB; - fluxes.eventOutputN += outputN / climLen; - } - // clang-format off - writeEventOut( - gEvent, 11, - "eventSoilC", soilAdd, - "eventLitterC", litterAdd, - "eventLeafC", leafDelta, - "eventWoodC", woodDelta, - "eventAccountingC", accountingDelta, - "eventFineRootC", fineDelta, - "eventCoarseRootC", coarseDelta, - "eventSoilOrgN", soilNAdd, - "eventLitterN", litterNAdd, - "eventOutputC", outputC, - "eventOutputN", outputN); - // clang-format on - } break; - case TILLAGE: { - // BIG NOTE: this is the one event type that is NOT modeled as a flux; - // see updateEventTrackers() for more - const TillageParams *tillParams = gEvent->eventParams; - // Update the tillage mod for R_H calculations; this will be slowly - // reduced by an exponential decay function. Note we add here, not set, - // as there may be lingering effects from a prior tillage. - eventTrackers.d_till_mod += tillParams->tillageEffect; - writeEventOut(gEvent, 1, "eventTrackers.d_till_mod", - tillParams->tillageEffect); - } break; - case FERTILIZATION: { - const FertilizationParams *fertParams = gEvent->eventParams; - const double orgC = fertParams->orgC; - double orgN = 0.0; - double minN = 0.0; - if (ctx.nitrogenCycle) { - orgN = fertParams->orgN; - minN = fertParams->minN; - } - 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 - // Implies ctx.litterPool - fluxes.eventLitterN += orgN / climLen; - fluxes.eventMinN += minN / climLen; - } - - // MASS BALANCE: this is a system input - fluxes.eventInputC += orgC / climLen; - if (ctx.nitrogenCycle) { - fluxes.eventInputN += (orgN + minN) / climLen; - } - - // clang-format off - writeEventOut(gEvent, 6, - "eventLitterC", ctx.litterPool ? orgC : 0.0, - "eventSoilC", ctx.litterPool ? 0.0 : orgC, + appendLog(gEvent, 3, "eventMinN", minN, "eventLitterN", orgN, - "eventInputC", orgC, "eventInputN", (orgN + minN)); // clang-format on } break; case LEAFON: { - double leafOnFlux = params.leafGrowth / climLen; - checkLeafOnLimitation(&leafOnFlux); - fluxes.eventLeafOnCreation += leafOnFlux; - double totalSourceC = envi.plantWoodC + envi.coarseRootC; - if (totalSourceC > TINY) { - fluxes.eventLeafOnCreationFromWood += - leafOnFlux * envi.plantWoodC / totalSourceC; - } - // 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 - // eventInputC or eventInputN - - // ALSO, unlike all other events, we don't write the event here, as N - // limitation may change the amount + // eventInputN } break; case LEAFOFF: { - double leafOff = envi.plantLeafC * params.fracLeafFall; - fluxes.eventLeafOffLitter += leafOff / climLen; - - double litterNAdd = 0.0; - double leafNResorption = 0.0; - if (ctx.nitrogenCycle) { - // Nitrogen - need to account for leaf N moving to litter, as with - // harvests - double leafN = leafOff / params.leafCN; - leafNResorption = leafN * params.leafNResorptionFrac; - litterNAdd = leafN - leafNResorption; - fluxes.eventLeafOffNResorption += leafNResorption / climLen; - fluxes.eventLitterN += litterNAdd / climLen; + // We need to use fluxes.eventLeafOffLitter here instead of + // recalculating it, as it may have been reduced. Note that this means + // we are processing ALL leaf-off events in one shot in the unlikely + // event that there are more than one in this time step. + if (!firstLeafOff) { + logInfo( + "Ignoring nitrogen effects of second (or more) leaf-off " + "event at year %d day %d; all nitrogen effects are recorded with " + "the first leaf-off event in this time step\n", + event->year, event->day); + break; } - + firstLeafOff = 0; + // Nitrogen - need to account for leaf N moving to litter, as with + // harvests + double leafOff = fluxes.eventLeafOffLitterC * climLen; + double preResorp = fluxes.eventLeafOffNResorption; + double preLitter = fluxes.eventLeafOffLitterN; + calcLeafOffNEffects(leafOff, &fluxes.eventLeafOffNResorption, + &fluxes.eventLeafOffLitterN); + + double leafNResorptionFlux = fluxes.eventLeafOffNResorption - preResorp; + double litterNAddFlux = fluxes.eventLeafOffLitterN - preLitter; + logInfo("Proc events: leaf N resorption %.4f litter N = %.4f\n", + leafNResorptionFlux * climLen, litterNAddFlux * climLen); // clang-format off - writeEventOut(gEvent, 3, - "eventLeafOffLitter", leafOff, - "eventLeafOffNResorption", leafNResorption, - "eventLitterN", litterNAdd); + appendLog(gEvent, 2, + "eventLeafOffNResorption", leafNResorptionFlux * climLen, + "eventLitterN", litterNAddFlux * climLen); // clang-format on } break; 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", - gEvent->year, gEvent->day); + // Nothing to do here break; default: - logError("Unknown event type (%d) in processEvents()\n", gEvent->type); + logError("Unknown event type (%d) in processEvents()\n", event->type); exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM); } - gEvent = gEvent->nextEvent; + event = event->nextEvent; } } @@ -1110,12 +891,12 @@ void updatePoolsForEvents(void) { double eventLeafOnCreationFromRoot = fluxes.eventLeafOnCreation - fluxes.eventLeafOnCreationFromWood; envi.coarseRootC -= eventLeafOnCreationFromRoot * climate->length; - envi.plantLeafC += (fluxes.eventLeafOnCreation - fluxes.eventLeafOffLitter) * + envi.plantLeafC += (fluxes.eventLeafOnCreation - fluxes.eventLeafOffLitterC) * climate->length; if (ctx.litterPool) { - envi.litterC += fluxes.eventLeafOffLitter * climate->length; + envi.litterC += fluxes.eventLeafOffLitterC * climate->length; } else { - envi.soilC += fluxes.eventLeafOffLitter * climate->length; + envi.soilC += fluxes.eventLeafOffLitterC * climate->length; } // Harvest and planting events @@ -1133,7 +914,11 @@ void updatePoolsForEvents(void) { if (ctx.nitrogenCycle) { envi.minN += fluxes.eventMinN * climate->length; envi.soilOrgN += fluxes.eventSoilOrgN * climate->length; - envi.litterN += fluxes.eventLitterN * climate->length; + envi.litterN += + (fluxes.eventLitterN + fluxes.eventLeafOffLitterN) * climate->length; + logInfo("Proc events: envi.litterN += %f\n", + (fluxes.eventLitterN + fluxes.eventLeafOffLitterN) * + climate->length); double leafOnNFlux = calcLeafOnNFromC(fluxes.eventLeafOnCreation); envi.plantStorageN += (fluxes.eventLeafOffNResorption - leafOnNFlux) * climate->length; @@ -1150,6 +935,9 @@ void freeEventList(void) { if (prev->eventParams != NULL) { free(prev->eventParams); } + if (prev->logLine != NULL) { + dsFree(prev->logLine); + } free(prev); } } diff --git a/src/sipnet/events.h b/src/sipnet/events.h index de757fe2..ed3aa6fd 100644 --- a/src/sipnet/events.h +++ b/src/sipnet/events.h @@ -1,6 +1,8 @@ #ifndef EVENTS_H #define EVENTS_H +#include "common/util.h" + typedef enum EventType { FERTILIZATION, HARVEST, @@ -82,7 +84,7 @@ struct EventNode { int year, day; void *eventParams; int numLogParamPairs; - char *logLine; + DynamicString *logLine; EventNode *nextEvent; }; @@ -119,29 +121,20 @@ EventNode *readEventData(const char *eventFile); void openEventOutFile(const char *eventOutFile, int printHeader); /*! - * \brief Write a line to the event output file for a single oneEvent - * - * Writes a single oneEvent to the configured event output file. This is a - * variadic function which expects to receive 2*numParams values in (char*, - * double) pairs after the - * numParams argument. - * - * Output format: - * - * year day event_type \=\[,\=\,...] - * - * \param oneEvent Pointer to oneEvent node - * \param numParams Number of param/value PAIRS to write - * \param ... Pairs of (char*, double) arguments to write, 2*numParams - * values + * Append to an event's log line */ -void writeEventOut(EventNode *oneEvent, int numParams, ...); +void appendLog(EventNode *event, int numParams, ...); + +/*! + * Write out all events for this time step + */ +void writeEventsOut(void); /*! * \brief Write a line to the event output file for a computed event * - * Same as writeEventOut, but for events that are computed internally, such - * as leaf on/leaf off events. + * Write an event that is computed internally, such as leaf on/leaf off or + * plant death events. * * Output format: * @@ -188,6 +181,11 @@ void setupEvents(void); */ int isFirstEventBefore(int year, int day); +/*! + * Return today's first event for SIPNET's multi-pass calls + */ +EventNode *getCurrentEvent(void); + /*! * \brief Process carbon effects from events for current day * @@ -209,11 +207,6 @@ void processEventsForCarbon(EventNode *event); */ void processEventsForNitrogen(EventNode *event); -/*! - * \brief Write events to the events output file for the current date - */ -void writeAllEvents(void); - /*! * Update relevant environment pools after event fluxes have been calculated */ diff --git a/src/sipnet/limitations.c b/src/sipnet/limitations.c index 1473b3dc..f83363d0 100644 --- a/src/sipnet/limitations.c +++ b/src/sipnet/limitations.c @@ -63,16 +63,110 @@ void checkLeafOnLimitation(double *leafOnFlux) { } } +/** + * Check that negative growth is not driving a pool to end negative + * + * Adjust if necessary + */ +static void checkNegativeCreation(void) { + // In the case of negative growth (mean npp < 0), we might be allocating that + // negative growth to a pool that can't handle it (e.g., leaf creation is + // negative, but leaf pool is already at 0). In those cases, adjust + // appropriately. + + double len = climate->length; + + logInfo("checkNeg BEFORE: leafLitter %f eventLOLitter %f leafC %f " + "leafCreation %f " + "eventLONResorp %f eventLOLitterN %f eventLitterN %f\n", + fluxes.leafLitter * len, fluxes.eventLeafOffLitterC * len, + envi.plantLeafC, fluxes.leafCreation * len, + fluxes.eventLeafOffNResorption * len, + fluxes.eventLeafOffLitterN * len, fluxes.eventLitterN); + + // Above ground + // If leafCreation is too negative, we need to deduct from wood instead + // Need to make sure we don't go too far the other way. Leaf off litter + // (either fluxes.leafLitter or fluxes.eventLeafOffLitter) might also + // incorrectly drive the pool negative, sp adjust for that too. + + // First we handle leaf litter. Assuming params.leafTurnoverRate is valid + // (ie, <=1), availableLeafRate should be non-negative. + double availableLeafRate = envi.plantLeafC / len - fluxes.leafLitter; + + // Reminder: litter and creation fluxes "point" in the opposite direction, so + // their signs are opposite here + double deficit = fluxes.leafOffLitter + fluxes.eventLeafOffLitterC - + fluxes.leafCreation - availableLeafRate; + if (deficit > 0) { + // If negative growth + leaf off is too much, let's first reduce leaf off + if (fluxes.eventLeafOffLitterC > 0) { + double adjust = fmin(fluxes.eventLeafOffLitterC, deficit); + fluxes.eventLeafOffLitterC -= adjust; + deficit -= adjust; + } + if (deficit > 0) { + if (fluxes.leafOffLitter > 0) { + double adjust = fmin(fluxes.leafOffLitter, deficit); + fluxes.leafOffLitter -= adjust; + deficit -= adjust; + } + } + // Next, move some negative growth to the wood pool + if (deficit > 0) { + // If this is too much for the wood pool to handle, we have an error that + // will be caught in ensureNonNegative + fluxes.woodCreation -= deficit; + fluxes.leafCreation += deficit; + } + } + + logInfo("checkNeg AFTER: leafLitter %f eventLOLitter %f leafC %f " + "leafCreation %f " + "eventLONResorp %f eventLOLitterN %f eventLitterN %f\n", + fluxes.leafLitter * len, fluxes.eventLeafOffLitterC * len, + envi.plantLeafC, fluxes.leafCreation * len, + fluxes.eventLeafOffNResorption * len, + fluxes.eventLeafOffLitterN * len, fluxes.eventLitterN); + + // Below ground + double fineRootDeficit = + envi.fineRootC / len + fluxes.fineRootCreation - fluxes.fineRootLoss; + double coarseRootDeficit = envi.coarseRootC / len + + fluxes.coarseRootCreation - fluxes.coarseRootLoss; + if ((fineRootDeficit < 0.0) != (coarseRootDeficit < 0.0)) { + // If neither are negative, nothing to do + // If both are negative, the plant will die in checkForMortality() + if (fineRootDeficit < 0.0) { + fluxes.coarseRootCreation += fineRootDeficit; + fluxes.fineRootCreation -= fineRootDeficit; + } + if (coarseRootDeficit < 0.0) { + fluxes.fineRootCreation += coarseRootDeficit; + fluxes.coarseRootCreation -= coarseRootDeficit; + } + } +} + /** * Check for nitrogen limitation, and reduce growth if needed */ static void checkNitrogenLimitation(void) { + double len = climate->length; + logInfo("checkNLimit BEFORE: leafLitter %f eventLOLitter %f leafC %f " + "leafCreation %f " + "eventLONResorp %f eventLOLitterN %f eventLitterN %f\n", + fluxes.leafLitter * len, fluxes.eventLeafOffLitterC * len, + envi.plantLeafC, fluxes.leafCreation * len, + fluxes.eventLeafOffNResorption * len, + fluxes.eventLeafOffLitterN * len, fluxes.eventLitterN); + // First, determine if we are in a nitrogen-limited situation. The uptake // flux has already taken the storage pool into account, so we just need to // see if that uptake is too much, taking into account other fluxes to the // minN pool. // Calc total delta to minN pool - double len = climate->length; + // double len = climate->length; double uptakeDemand = fluxes.nUptake * len; double nonUptakeDelta = calcMinNNonUptakeFluxes() * len; double availableMinN = envi.minN + nonUptakeDelta; @@ -108,9 +202,48 @@ static void checkNitrogenLimitation(void) { fluxes.fineRootCreation *= reduction; fluxes.coarseRootCreation *= reduction; - // Reset fixation and uptake - calcNFixationAndUptakeFluxes(); + // If there was a leaf-off event in this step, it is possible that we have + // reduced leaf growth down to the point where the pool will go negative + // Re-call that check + checkNegativeCreation(); + // double leafOffFlux = fluxes.leafOffLitter + fluxes.eventLeafOffLitter; + // if (leafOffFlux > 0) { + // logInfo("Leaf-off flux %.f, leaf pool %.f, leaf creation %.f leaf + // litter %.f\n", + // leafOffFlux * len, envi.plantLeafC, fluxes.leafCreation * len, + // fluxes.leafLitter * len); + // double leafPool = envi.plantLeafC + + // (fluxes.leafCreation + fluxes.leafLitter + leafOffFlux) * len; + // if (leafPool < 0) { + // // Reduce leaf-off litter to avoid negative pool; note that only one + // of + // // these will be > 0, meaning the other is unaffected by this calc + // fluxes.leafOffLitter = fmax(0.0, fluxes.leafOffLitter + leafPool); + // fluxes.eventLeafOffLitter = + // fmax(0.0, fluxes.eventLeafOffLitter + leafPool); + // } + // // Recalc N effects for event leaf off; will be a no-op if there was + // // no event leaf off + + // Reset and recalc event leaf-off N + fluxes.eventLeafOffNResorption = 0.0; + fluxes.eventLeafOffLitterN = 0.0; + calcLeafOffNEffects(fluxes.eventLeafOffLitterC * len, + &fluxes.eventLeafOffNResorption, + &fluxes.eventLeafOffLitterN); + // } + + // Reset N calculations + calcNitrogenFluxes(); } + + logInfo("checkNLimit AFTER: leafLitter %f eventLOLitter %f leafC %f " + "leafCreation %f " + "eventLONResorp %f eventLOLitterN %f eventLitterN %f\n", + fluxes.leafLitter * len, fluxes.eventLeafOffLitterC * len, + envi.plantLeafC, fluxes.leafCreation * len, + fluxes.eventLeafOffNResorption * len, + fluxes.eventLeafOffLitterN * len, fluxes.eventLitterN); } /** @@ -138,85 +271,5 @@ void checkLimitations(void) { } } -/** - * Check that negative growth is not driving a pool to end negative - * - * Adjust if necessary - */ -static void checkNegativeCreation(void) { - // In the case of negative growth (mean npp < 0), we might be allocating that - // negative growth to a pool that can't handle it (e.g., leaf creation is - // negative, but leaf pool is already at 0). In those cases, adjust - // appropriately. - - double len = climate->length; - - logInfo("BEFORE: leafLitter %f eventLOLitter %f leafC %f leafCreation %f\n", - fluxes.leafLitter * len, fluxes.eventLeafOffLitter * len, - envi.plantLeafC, fluxes.leafCreation * len); - - // Above ground - // If leafCreation is too negative, we need to deduct from wood instead - // Need to make sure we don't go too far the other way. Leaf off litter - // (either fluxes.leafLitter or fluxes.eventLeafOffLitter) might also - // incorrectly drive the pool negative, sp adjust for that too. - - // First we handle leaf litter. Assuming params.leafTurnoverRate is valid - // (ie, <=1), availableLeafRate should be non-negative. - double availableLeafRate = envi.plantLeafC / len - fluxes.leafLitter; - - // Reminder: litter and creation fluxes "point" in the opposite direction, so - // their signs are opposite here - double deficit = fluxes.leafOffLitter + fluxes.eventLeafOffLitter - - fluxes.leafCreation - availableLeafRate; - if (deficit > 0) { - // If negative growth + leaf off is too much, let's first reduce leaf off - if (fluxes.eventLeafOffLitter > 0) { - double adjust = fmin(fluxes.eventLeafOffLitter, deficit); - fluxes.eventLeafOffLitter -= adjust; - deficit -= adjust; - // If we change EVENT leaf off litter, we need to adjust eventLitterN, as - // it has already been calculated - fluxes.eventLitterN - } - if (deficit > 0) { - if (fluxes.leafOffLitter > 0) { - double adjust = fmin(fluxes.leafOffLitter, deficit); - fluxes.leafOffLitter -= adjust; - deficit -= adjust; - } - } - // Next, move some negative growth to the wood pool - if (deficit > 0) { - // If this is too much for the wood pool to handle, we have an error that - // will be caught in ensureNonNegative - fluxes.woodCreation -= deficit; - fluxes.leafCreation += deficit; - } - } - - logInfo("AFTER: leafLitter %f eventLOLitter %f leafC %f leafCreation %f\n", - fluxes.leafLitter * len, fluxes.eventLeafOffLitter * len, - envi.plantLeafC, fluxes.leafCreation * len); - - // Below ground - double fineRootDeficit = - envi.fineRootC / len + fluxes.fineRootCreation - fluxes.fineRootLoss; - double coarseRootDeficit = envi.coarseRootC / len + - fluxes.coarseRootCreation - fluxes.coarseRootLoss; - if ((fineRootDeficit < 0.0) != (coarseRootDeficit < 0.0)) { - // If neither are negative, nothing to do - // If both are negative, the plant will die in checkForMortality() - if (fineRootDeficit < 0.0) { - fluxes.coarseRootCreation += fineRootDeficit; - fluxes.fineRootCreation -= fineRootDeficit; - } - if (coarseRootDeficit < 0.0) { - fluxes.fineRootCreation += coarseRootDeficit; - fluxes.coarseRootCreation -= coarseRootDeficit; - } - } -} - // See limitations.h void checkCarbonLimitations(void) { checkNegativeCreation(); } diff --git a/src/sipnet/nitrogen.c b/src/sipnet/nitrogen.c index 8c9ffbf1..ccb70a85 100644 --- a/src/sipnet/nitrogen.c +++ b/src/sipnet/nitrogen.c @@ -72,6 +72,13 @@ static void calcNPoolFluxes(void) { fluxes.woodLitter / params.woodCN - litterMin - fluxes.litterToSoil / litterCN + (soilNInputs * saturationFraction); + double len = climate->length; + logInfo("NPoolFluxes: leaf litter %f, leafOffNResorp %f, litterMin %f " + "nOrgLitter %f\n", + getLeafLitterFlux() / params.leafCN * len, + fluxes.leafOffNResorption * len, litterMin * len, + fluxes.nOrgLitter * len); + // soil // The soil org N flux is determined by the carbon flux from the litter pool, // carbon fluxes from roots, and N loss due to mineralization @@ -190,7 +197,7 @@ void calcNResorptionFluxes(void) { // in events.c double nResorp = params.leafNResorptionFrac * getLeafLitterFlux() / params.leafCN; - fluxes.leafOffNResorption += nResorp; + fluxes.leafOffNResorption = nResorp; // TODO: Should we resorb N from wood litter? } @@ -236,4 +243,19 @@ void updateNitrogenPools(void) { // Litter organic N envi.litterN += fluxes.nOrgLitter * climate->length; + logInfo("Update N Pools: envi.litterN += %f\n", + fluxes.nOrgLitter * climate->length); +} + +void calcLeafOffNEffects(double leafOffC, double *resorptionFlux, + double *litterFlux) { + double climLen = climate->length; + double leafN = leafOffC / params.leafCN; + double leafNResorption = leafN * params.leafNResorptionFrac; + double litterNAdd = leafN - leafNResorption; + *resorptionFlux += leafNResorption / climLen; + *litterFlux += litterNAdd / climLen; + + logInfo("calcLONEffects: resorpN %f litterN %f\n", *resorptionFlux * climLen, + *litterFlux * climLen); } diff --git a/src/sipnet/nitrogen.h b/src/sipnet/nitrogen.h index db7a475b..8a30830a 100644 --- a/src/sipnet/nitrogen.h +++ b/src/sipnet/nitrogen.h @@ -78,4 +78,13 @@ void updateNitrogenPools(void); */ void calcNResorptionFluxes(void); +/** + * Helper function for handling leaf-off events + * + * This is used by both processEventsForNitrogen and the nitrogen limitation + * check in limitations. + */ +void calcLeafOffNEffects(double leafOffC, double *resorptionFlux, + double *litterFlux); + #endif // NITROGEN_H diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 5014e4e3..068e3259 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -795,6 +795,7 @@ void calcWoodAndLeafFluxes(void) { 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): @@ -807,23 +808,26 @@ 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 / climate->length; + double leafOn = params.leafGrowth / len; + double leafOnFromWood = 0; checkLeafOnLimitation(&leafOn); fluxes.leafOnCreation += leafOn; double totalSourceC = envi.plantWoodC + envi.coarseRootC; if (totalSourceC > TINY) { - fluxes.leafOnCreationFromWood += leafOn * envi.plantWoodC / totalSourceC; + leafOnFromWood = leafOn * envi.plantWoodC / totalSourceC; + fluxes.leafOnCreationFromWood += leafOnFromWood; } phenologyTrackers.didLeafGrowth = 1; - // This is a computed event - however, the value may get reduced by - // nitrogen limitation. The writeEvent call is in writeLeafOnEventIfNeeded, - // called after N limiting is checked. + + writeComputedEventOut(climate->year, climate->day, + eventTypeToString(LEAFON), 2, "leafOnCreation", + leafOn * len, "leafOnCreationFromWood", + leafOnFromWood * len); } // check for end of growing season: if (!phenologyTrackers.didLeafFall && pastLeafFall()) { // we just reached the end of the growing season - double len = climate->length; double leafOff = (envi.plantLeafC * params.fracLeafFall) / len; fluxes.leafOffLitter += leafOff; phenologyTrackers.didLeafFall = 1; @@ -1215,31 +1219,6 @@ void calcMethaneFlux(void) { */ void resetFluxes(void) { fluxes = (struct FluxVars){0}; } -/** - * Write out a leaf-on event if one happened - * - * Delayed event writing for leaf-on, if appropriate, since the value may - * have changed due to N limitation - */ -void writeLeafOnEventIfNeeded(void) { - const char *type = eventTypeToString(LEAFON); - const double len = climate->length; - if (fluxes.leafOnCreation > TINY && ctx.events) { - writeComputedEventOut(climate->year, climate->day, type, 2, - "leafOnCreation", fluxes.leafOnCreation * len, - "leafOnCreationFromWood", - fluxes.leafOnCreationFromWood * len); - } - if (fluxes.eventLeafOnCreation > TINY && ctx.events) { - // Not really a computed event, but we don't have the event object here, so - // we use this mechanism - writeComputedEventOut( - climate->year, climate->day, type, 2, "eventLeafOnCreation", - fluxes.eventLeafOnCreation * len, "eventLeafOnCreationFromWood", - fluxes.eventLeafOnCreationFromWood * len); - } -} - /*! * Calculate flux terms for sipnet as part of main model flow * @@ -1267,6 +1246,9 @@ void calculateFluxes(void) { // Psn, moisture and water fluxes lai = envi.plantLeafC / params.leafCSpWt; // current lai + // Today's first event + EventNode *event = getCurrentEvent(); + potPsn(&potGrossPsn, &baseFolResp, lai, climate->tair, climate->vpd, climate->par); moisture(&(fluxes.transpiration), &dWater, potGrossPsn, climate->vpd, @@ -1279,6 +1261,9 @@ void calculateFluxes(void) { fluxes.snowMelt, fluxes.transpiration); getGpp(&(fluxes.photosynthesis), potGrossPsn, dWater); + // First pass for events: carbon effects + processEventsForCarbon(event); + // Vegetation respiration if (ctx.growthResp) { vegResp2(&folResp, &woodResp, &growthResp, baseFolResp); @@ -1317,15 +1302,17 @@ void calculateFluxes(void) { // so make sure this stays at the bottom of this function (or after the // carbon and water calcs, at least). if (ctx.nitrogenCycle) { + // Second pass for events: nitrogen effects + processEventsForNitrogen(event); + // General nitrogen fluxes calcNitrogenFluxes(); } // Check limitations; must be after all flux calcs checkLimitations(); - // Delayed write needed due to N limitation check - leafOn value may have - // changed - writeLeafOnEventIfNeeded(); + // Write to events file for all of today's events + writeEventsOut(); } // /////////////////////// // @@ -1468,7 +1455,7 @@ void updateTrackers(double oldSoilWater) { // If we get another event flux in this function, we should create an // updateTrackersForEvents() function in events.c|h - trackers.yearlyLitter += getLeafLitterFlux() + fluxes.eventLeafOffLitter; + trackers.yearlyLitter += getLeafLitterFlux() + fluxes.eventLeafOffLitterC; if (ctx.gdd) { trackers.gdd += climate->gdd; @@ -1844,12 +1831,7 @@ void updateState(void) { /////////////////////// // 1. Calculate Fluxes - // All event handling, which is modeled as fluxes. Note that we have this - // before the other fluxes so that everything is in place when we consider - // N limitation at the end of calculateFluxes(). - processEvents(); - - // All non-event fluxes + // All modeled fluxes, including events calculateFluxes(); /////////////////////// diff --git a/src/sipnet/state.h b/src/sipnet/state.h index d341eade..73ef466f 100644 --- a/src/sipnet/state.h +++ b/src/sipnet/state.h @@ -636,7 +636,9 @@ typedef struct FluxVars { // coarse root C) double eventLeafOnCreationFromWood; // Transfer from leafC to soil/litter from a leaf-off event - double eventLeafOffLitter; + double eventLeafOffLitterC; + // Transfer from leafN to soil/litter from a leaf-off event + double eventLeafOffLitterN; // Resorption of leaf N from a leaf-off event double eventLeafOffNResorption; diff --git a/tests/sipnet/test_modeling/testCompleteHarvest.c b/tests/sipnet/test_modeling/testCompleteHarvest.c index d2d48ed6..b07e4438 100644 --- a/tests/sipnet/test_modeling/testCompleteHarvest.c +++ b/tests/sipnet/test_modeling/testCompleteHarvest.c @@ -24,6 +24,12 @@ static void balanced(void) { near(balanceTracker.clampedN, 0, "N clamping"); } +static void processEvents(void) { + EventNode *event = getCurrentEvent(); + processEventsForCarbon(event); + processEventsForNitrogen(event); +} + static void start(int mode, const char *harvest) { initContext(); ctx.events = 1; @@ -144,11 +150,14 @@ static void fullCase(int mode, int accounting, int dark, int limited, near(plantSurvivalTracker.isAlive, 1, "alive after planting"); } finish(); + + logInfo(":\n"); + if (failures > 0) { + logTest("Complete harvest failures: %d\n", failures); + logTest("Press return to continue\n"); + getchar(); + } } - // if (mode == 2) { - // logTest("Press return to continue\n"); - // getchar(); - // } } static void partialCase(int accounting, const char *harvest, double fraction) { @@ -289,7 +298,7 @@ static void coincidentEventCase(int type, const char *arguments) { end = envi; ordinary = fluxes; if ((type == LEAFON && fluxes.eventLeafOnCreation <= 0) || - (type == LEAFOFF && fluxes.eventLeafOffLitter <= 0) || + (type == LEAFOFF && fluxes.eventLeafOffLitterC <= 0) || (type == IRRIGATION && fluxes.eventSoilWater <= 0)) failures++; } else { @@ -321,7 +330,7 @@ static void coincidentEventCase(int type, const char *arguments) { F(eventEvap); F(eventLeafOnCreation); F(eventLeafOnCreationFromWood); - F(eventLeafOffLitter); + F(eventLeafOffLitterC); F(eventLeafOffNResorption); #undef F climate = climate->nextClim; @@ -376,18 +385,32 @@ static void restartCase(void) { } static void leafBudgetCases(void) { + int caseNum = 0; for (int kind = 0; kind < 4; kind++) { + // kind: 0 1 2 3 + // fracLeafFall 1 .25 .75 .75 + // leafOff a a b c + // a: inserted at start + // b: two events inserted at start + // c: no leaf-off events for (int dark = 0; dark < 2; dark++) { + // dark : set climate->par=0 when dark for (int sign = -1; sign <= 1; sign++) { + // sign: init npp tracker with 20*sign (-20, 0, 20) for (int account = -1; account <= 1; account++) { + // account: starting value for plantCAccountingDelta for (int limited = 0; limited < 2; limited++) { + // limited: starting minN = 0 if limited (1000 else) for (int resorb = 0; resorb < 3; resorb++) { + /// resorb: leafNResorptionFrac = 0, .5, 1 (resorb/2) for (int harvest = 0; harvest < 2; harvest++) { - logTest("*** Running leaf budget case with kind: %d dark: %d " - "sign %d" - " account %d limited %d resorb %d harvest %s\n", - kind, dark, sign, account, limited, resorb, - harvest ? "0 0 1 1" : "none"); + // TODO: Add "1 1 0 0" and "0.5 0.5 0.5 0.5" cases for harvest + logTest( + "*** Running leaf budget case [%d] with kind: %d dark: %d " + "sign %d" + " account %d limited %d resorb %.1f harvest %s\n", + caseNum++, kind, dark, sign, account, limited, resorb / 2.0, + harvest ? "0 0 1 1" : "none"); start(2, harvest ? "0 0 1 1" : NULL); envi.plantLeafC = 1; envi.plantWoodC = envi.fineRootC = envi.coarseRootC = 100; @@ -419,30 +442,37 @@ static void leafBudgetCases(void) { } setupEvents(); } + double beforeLeafC = envi.plantLeafC; double beforeLitter = envi.litterC; double beforeLitterN = envi.litterN; updateState(); balanced(); double shed = (fluxes.leafLitter + fluxes.leafOffLitter + - fluxes.eventLeafOffLitter) * + fluxes.eventLeafOffLitterC) * climate->length; - if (shed < -1e-10 || shed > 1 + 1e-10 || - envi.plantLeafC < -1e-10) { - logTest("Leaf budget overdraw: %.17g, leaf %.17g\n", shed, + double pool = + envi.plantLeafC + fluxes.leafCreation * climate->length; + if (fabs(pool - shed) < 1e-10 || envi.plantLeafC < -1e-10) { + logTest("Leaf budget overdraw (2): %.17g, leaf %.17g\n", shed, envi.plantLeafC); failures++; } - double expEventLOL[] = {1, .25, 1, 0}; - double eventExpected = expEventLOL[kind]; - near(fluxes.eventLeafOffLitter * climate->length, eventExpected, - "event leaves conserved"); + double expLeafOffLitter = + kind == 3 ? 0 + : fmin(1.0, beforeLeafC + (fluxes.leafCreation - + fluxes.leafLitter) * + climate->length); + near(fluxes.eventLeafOffLitterC * climate->length, + expLeafOffLitter, "event leaves conserved"); if (!harvest) { near(envi.litterC - beforeLitter, shed, "leaf litter transfer"); near(envi.litterN - beforeLitterN, shed * (1 - params.leafNResorptionFrac) / params.leafCN, "leaf litter N transfer"); + logTest("beforeLitterN %f afterLitterN %f\n", beforeLitterN, + envi.litterN); } else { clearPlant(); climate = climate->nextClim; @@ -453,15 +483,21 @@ static void leafBudgetCases(void) { "no regrowth after complete harvest"); } finish(); + + logInfo(":\n"); + if (failures > 0) { + logTest("Complete harvest failures: %d\n", failures); + logTest("Press return to continue\n"); + getchar(); + } } // harvest loop } // resorb loop - - logTest("Complete harvest failures: %d\n", failures); - logTest("Press return to continue\n"); - getchar(); - } // limited loop } // account loop + logTest("New sign value incoming\n"); + logTest("Complete harvest failures: %d\n", failures); + logTest("Press return to continue\n"); + getchar(); } // sign loop } // dark loop } // kind loop From 27aef2ab29488ea709eae9005209ea0f9499e72e Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Mon, 28 Sep 2026 17:14:48 -0400 Subject: [PATCH 05/14] Interim commit, through leaf case 503 --- src/sipnet/limitations.c | 6 ++-- src/sipnet/nitrogen.c | 6 ++-- .../test_modeling/testCompleteHarvest.c | 29 +++++++++++-------- 3 files changed, 24 insertions(+), 17 deletions(-) diff --git a/src/sipnet/limitations.c b/src/sipnet/limitations.c index f83363d0..ae9ce9c3 100644 --- a/src/sipnet/limitations.c +++ b/src/sipnet/limitations.c @@ -78,11 +78,13 @@ static void checkNegativeCreation(void) { logInfo("checkNeg BEFORE: leafLitter %f eventLOLitter %f leafC %f " "leafCreation %f " - "eventLONResorp %f eventLOLitterN %f eventLitterN %f\n", + "eventLONResorp %f eventLOLitterN %f eventLitterN %f " + "params.fracLeafFall %f\n", fluxes.leafLitter * len, fluxes.eventLeafOffLitterC * len, envi.plantLeafC, fluxes.leafCreation * len, fluxes.eventLeafOffNResorption * len, - fluxes.eventLeafOffLitterN * len, fluxes.eventLitterN); + fluxes.eventLeafOffLitterN * len, fluxes.eventLitterN, + params.fracLeafFall); // Above ground // If leafCreation is too negative, we need to deduct from wood instead diff --git a/src/sipnet/nitrogen.c b/src/sipnet/nitrogen.c index ccb70a85..8d681ab7 100644 --- a/src/sipnet/nitrogen.c +++ b/src/sipnet/nitrogen.c @@ -73,11 +73,11 @@ static void calcNPoolFluxes(void) { fluxes.litterToSoil / litterCN + (soilNInputs * saturationFraction); double len = climate->length; - logInfo("NPoolFluxes: leaf litter %f, leafOffNResorp %f, litterMin %f " - "nOrgLitter %f\n", + logInfo("NPoolFluxes: leafLitterN %f, leafOffNResorp %f, litterMin %f " + "nOrgLitter %f params.leafNResorptionFrac %f\n", getLeafLitterFlux() / params.leafCN * len, fluxes.leafOffNResorption * len, litterMin * len, - fluxes.nOrgLitter * len); + fluxes.nOrgLitter * len, params.leafNResorptionFrac); // soil // The soil org N flux is determined by the carbon flux from the litter pool, diff --git a/tests/sipnet/test_modeling/testCompleteHarvest.c b/tests/sipnet/test_modeling/testCompleteHarvest.c index b07e4438..4a306b60 100644 --- a/tests/sipnet/test_modeling/testCompleteHarvest.c +++ b/tests/sipnet/test_modeling/testCompleteHarvest.c @@ -412,6 +412,7 @@ static void leafBudgetCases(void) { caseNum++, kind, dark, sign, account, limited, resorb / 2.0, harvest ? "0 0 1 1" : "none"); start(2, harvest ? "0 0 1 1" : NULL); + logInfo("* Post start setup *\n"); envi.plantLeafC = 1; envi.plantWoodC = envi.fineRootC = envi.coarseRootC = 100; envi.plantCAccountingDelta = account; @@ -450,19 +451,21 @@ static void leafBudgetCases(void) { double shed = (fluxes.leafLitter + fluxes.leafOffLitter + fluxes.eventLeafOffLitterC) * climate->length; - double pool = - envi.plantLeafC + fluxes.leafCreation * climate->length; - if (fabs(pool - shed) < 1e-10 || envi.plantLeafC < -1e-10) { - logTest("Leaf budget overdraw (2): %.17g, leaf %.17g\n", shed, - envi.plantLeafC); + double pool = 1 + fluxes.leafCreation * climate->length; + if (pool - shed < -1e-10 || envi.plantLeafC < -1e-10) { + logTest("Leaf budget overdraw (2): shed %.17g pool %.17g " + "|s-p| %.17g leaf %.17g expr1 %d expr2 %d\n", + shed, pool, pool - shed, envi.plantLeafC, + pool - shed < -1e-10, envi.plantLeafC < -1e-10); failures++; } double expLeafOffLitter = kind == 3 ? 0 - : fmin(1.0, beforeLeafC + (fluxes.leafCreation - - fluxes.leafLitter) * - climate->length); + : fmin(params.fracLeafFall, + beforeLeafC + (fluxes.leafCreation - + fluxes.leafLitter) * + climate->length); near(fluxes.eventLeafOffLitterC * climate->length, expLeafOffLitter, "event leaves conserved"); if (!harvest) { @@ -494,10 +497,12 @@ static void leafBudgetCases(void) { } // resorb loop } // limited loop } // account loop - logTest("New sign value incoming\n"); - logTest("Complete harvest failures: %d\n", failures); - logTest("Press return to continue\n"); - getchar(); + if (kind > 1) { + logTest("New sign value incoming\n"); + logTest("Complete harvest failures: %d\n", failures); + logTest("Press return to continue\n"); + getchar(); + } } // sign loop } // dark loop } // kind loop From e44b2532c9714be1f6184586f9bd34e8b0741723 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Tue, 29 Sep 2026 12:32:46 -0400 Subject: [PATCH 06/14] Interim commit: all test points passing --- src/sipnet/events.c | 18 ++- src/sipnet/limitations.c | 51 ------- src/sipnet/nitrogen.c | 12 -- .../test_modeling/testCompleteHarvest.c | 134 ++++++++---------- 4 files changed, 70 insertions(+), 145 deletions(-) diff --git a/src/sipnet/events.c b/src/sipnet/events.c index 0edc556a..1e1a938b 100644 --- a/src/sipnet/events.c +++ b/src/sipnet/events.c @@ -483,6 +483,9 @@ void processEventsForCarbon(EventNode *event) { // Reset harvest tracking eventTrackers.harvestTrackers = (HarvestTrackers){0}; + // Make sure we don't have more than 100% leaves fall for leaf-off + double totalLeafFallFrac = 0.0; + while (event != NULL && event->year <= climYear && event->day <= climDay) { // The events file has been tested on read, so we know this event list // should be in chrono order. However, we need to check to make sure the @@ -693,6 +696,15 @@ void processEventsForCarbon(EventNode *event) { // clang-format on } break; case LEAFOFF: { + double frac = params.fracLeafFall; + totalLeafFallFrac += frac; + if (totalLeafFallFrac > 1.0) { + 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); + exit(EXIT_CODE_BAD_PARAMETER_VALUE); + } + double leafOff = envi.plantLeafC * params.fracLeafFall; fluxes.eventLeafOffLitterC += leafOff / climLen; @@ -852,8 +864,7 @@ void processEventsForNitrogen(EventNode *event) { double leafNResorptionFlux = fluxes.eventLeafOffNResorption - preResorp; double litterNAddFlux = fluxes.eventLeafOffLitterN - preLitter; - logInfo("Proc events: leaf N resorption %.4f litter N = %.4f\n", - leafNResorptionFlux * climLen, litterNAddFlux * climLen); + // clang-format off appendLog(gEvent, 2, "eventLeafOffNResorption", leafNResorptionFlux * climLen, @@ -916,9 +927,6 @@ void updatePoolsForEvents(void) { envi.soilOrgN += fluxes.eventSoilOrgN * climate->length; envi.litterN += (fluxes.eventLitterN + fluxes.eventLeafOffLitterN) * climate->length; - logInfo("Proc events: envi.litterN += %f\n", - (fluxes.eventLitterN + fluxes.eventLeafOffLitterN) * - climate->length); double leafOnNFlux = calcLeafOnNFromC(fluxes.eventLeafOnCreation); envi.plantStorageN += (fluxes.eventLeafOffNResorption - leafOnNFlux) * climate->length; diff --git a/src/sipnet/limitations.c b/src/sipnet/limitations.c index ae9ce9c3..a074c463 100644 --- a/src/sipnet/limitations.c +++ b/src/sipnet/limitations.c @@ -76,16 +76,6 @@ static void checkNegativeCreation(void) { double len = climate->length; - logInfo("checkNeg BEFORE: leafLitter %f eventLOLitter %f leafC %f " - "leafCreation %f " - "eventLONResorp %f eventLOLitterN %f eventLitterN %f " - "params.fracLeafFall %f\n", - fluxes.leafLitter * len, fluxes.eventLeafOffLitterC * len, - envi.plantLeafC, fluxes.leafCreation * len, - fluxes.eventLeafOffNResorption * len, - fluxes.eventLeafOffLitterN * len, fluxes.eventLitterN, - params.fracLeafFall); - // Above ground // If leafCreation is too negative, we need to deduct from wood instead // Need to make sure we don't go too far the other way. Leaf off litter @@ -123,14 +113,6 @@ static void checkNegativeCreation(void) { } } - logInfo("checkNeg AFTER: leafLitter %f eventLOLitter %f leafC %f " - "leafCreation %f " - "eventLONResorp %f eventLOLitterN %f eventLitterN %f\n", - fluxes.leafLitter * len, fluxes.eventLeafOffLitterC * len, - envi.plantLeafC, fluxes.leafCreation * len, - fluxes.eventLeafOffNResorption * len, - fluxes.eventLeafOffLitterN * len, fluxes.eventLitterN); - // Below ground double fineRootDeficit = envi.fineRootC / len + fluxes.fineRootCreation - fluxes.fineRootLoss; @@ -155,13 +137,6 @@ static void checkNegativeCreation(void) { */ static void checkNitrogenLimitation(void) { double len = climate->length; - logInfo("checkNLimit BEFORE: leafLitter %f eventLOLitter %f leafC %f " - "leafCreation %f " - "eventLONResorp %f eventLOLitterN %f eventLitterN %f\n", - fluxes.leafLitter * len, fluxes.eventLeafOffLitterC * len, - envi.plantLeafC, fluxes.leafCreation * len, - fluxes.eventLeafOffNResorption * len, - fluxes.eventLeafOffLitterN * len, fluxes.eventLitterN); // First, determine if we are in a nitrogen-limited situation. The uptake // flux has already taken the storage pool into account, so we just need to @@ -208,24 +183,6 @@ static void checkNitrogenLimitation(void) { // reduced leaf growth down to the point where the pool will go negative // Re-call that check checkNegativeCreation(); - // double leafOffFlux = fluxes.leafOffLitter + fluxes.eventLeafOffLitter; - // if (leafOffFlux > 0) { - // logInfo("Leaf-off flux %.f, leaf pool %.f, leaf creation %.f leaf - // litter %.f\n", - // leafOffFlux * len, envi.plantLeafC, fluxes.leafCreation * len, - // fluxes.leafLitter * len); - // double leafPool = envi.plantLeafC + - // (fluxes.leafCreation + fluxes.leafLitter + leafOffFlux) * len; - // if (leafPool < 0) { - // // Reduce leaf-off litter to avoid negative pool; note that only one - // of - // // these will be > 0, meaning the other is unaffected by this calc - // fluxes.leafOffLitter = fmax(0.0, fluxes.leafOffLitter + leafPool); - // fluxes.eventLeafOffLitter = - // fmax(0.0, fluxes.eventLeafOffLitter + leafPool); - // } - // // Recalc N effects for event leaf off; will be a no-op if there was - // // no event leaf off // Reset and recalc event leaf-off N fluxes.eventLeafOffNResorption = 0.0; @@ -238,14 +195,6 @@ static void checkNitrogenLimitation(void) { // Reset N calculations calcNitrogenFluxes(); } - - logInfo("checkNLimit AFTER: leafLitter %f eventLOLitter %f leafC %f " - "leafCreation %f " - "eventLONResorp %f eventLOLitterN %f eventLitterN %f\n", - fluxes.leafLitter * len, fluxes.eventLeafOffLitterC * len, - envi.plantLeafC, fluxes.leafCreation * len, - fluxes.eventLeafOffNResorption * len, - fluxes.eventLeafOffLitterN * len, fluxes.eventLitterN); } /** diff --git a/src/sipnet/nitrogen.c b/src/sipnet/nitrogen.c index 8d681ab7..70ca096e 100644 --- a/src/sipnet/nitrogen.c +++ b/src/sipnet/nitrogen.c @@ -72,13 +72,6 @@ static void calcNPoolFluxes(void) { fluxes.woodLitter / params.woodCN - litterMin - fluxes.litterToSoil / litterCN + (soilNInputs * saturationFraction); - double len = climate->length; - logInfo("NPoolFluxes: leafLitterN %f, leafOffNResorp %f, litterMin %f " - "nOrgLitter %f params.leafNResorptionFrac %f\n", - getLeafLitterFlux() / params.leafCN * len, - fluxes.leafOffNResorption * len, litterMin * len, - fluxes.nOrgLitter * len, params.leafNResorptionFrac); - // soil // The soil org N flux is determined by the carbon flux from the litter pool, // carbon fluxes from roots, and N loss due to mineralization @@ -243,8 +236,6 @@ void updateNitrogenPools(void) { // Litter organic N envi.litterN += fluxes.nOrgLitter * climate->length; - logInfo("Update N Pools: envi.litterN += %f\n", - fluxes.nOrgLitter * climate->length); } void calcLeafOffNEffects(double leafOffC, double *resorptionFlux, @@ -255,7 +246,4 @@ void calcLeafOffNEffects(double leafOffC, double *resorptionFlux, double litterNAdd = leafN - leafNResorption; *resorptionFlux += leafNResorption / climLen; *litterFlux += litterNAdd / climLen; - - logInfo("calcLONEffects: resorpN %f litterN %f\n", *resorptionFlux * climLen, - *litterFlux * climLen); } diff --git a/tests/sipnet/test_modeling/testCompleteHarvest.c b/tests/sipnet/test_modeling/testCompleteHarvest.c index 4a306b60..d3d3713a 100644 --- a/tests/sipnet/test_modeling/testCompleteHarvest.c +++ b/tests/sipnet/test_modeling/testCompleteHarvest.c @@ -76,10 +76,11 @@ static void fullCase(int mode, int accounting, int dark, int limited, const char *harvest, double aboveExport, double belowExport) { // A no-harvest run supplies the independently calculated end-of-step pools. - logTest("Running full case with mode: %d accounting: %d dark: %d limited: %d" - " harvest: %s aboveExport: %.2f belowExport: %.2f\n", - mode, accounting, dark, limited, harvest ? harvest : "none", - aboveExport, belowExport); + logTest( + "*** Running full case with mode: %d accounting: %d dark: %d limited: %d" + " harvest: %s aboveExport: %.2f belowExport: %.2f\n", + mode, accounting, dark, limited, harvest ? harvest : "none", aboveExport, + belowExport); Envi end = {0}; for (int h = 0; h < 2; h++) { start(mode, h ? harvest : NULL); @@ -95,14 +96,7 @@ static void fullCase(int mode, int accounting, int dark, int limited, if (!h) { // First time through, capture envi state end = envi; - // logTest("(end) nLeach %f nVol %f\n", - // fluxes.nLeaching * climate->length, fluxes.nVolatilization * - // climate->length); - // logTest("(end) fluxes.eventOutputN %f\n", fluxes.eventOutputN * - // climate->length); } else { - // logTest("(envi) fluxes.eventOutputN %f\n", fluxes.eventOutputN * - // climate->length); const double aboveC = end.plantLeafC + end.plantWoodC + end.plantCAccountingDelta; const double belowC = end.fineRootC + end.coarseRootC; @@ -122,10 +116,6 @@ static void fullCase(int mode, int accounting, int dark, int limited, aboveC * aboveExport + belowC * belowExport, "C export"); near(fluxes.eventOutputN * climate->length, aboveN * aboveExport + belowN * belowExport, "N export"); - // logTest("N export: above %f fracAbove %f below %f fracbelow %f " - // "eventOutputN %f\n", - // aboveN, aboveExport, belowN, belowExport, fluxes.eventOutputN * - // climate->length); if (mode == 2) { near(envi.litterN, @@ -150,19 +140,13 @@ static void fullCase(int mode, int accounting, int dark, int limited, near(plantSurvivalTracker.isAlive, 1, "alive after planting"); } finish(); - - logInfo(":\n"); - if (failures > 0) { - logTest("Complete harvest failures: %d\n", failures); - logTest("Press return to continue\n"); - getchar(); - } } } static void partialCase(int accounting, const char *harvest, double fraction) { - logTest("Running partial case with accounting: %d harvest: %s fraction: %f\n", - accounting, harvest, fraction); + logTest( + "*** Running partial case with accounting: %d harvest: %s fraction: %f\n", + accounting, harvest, fraction); start(2, harvest); envi.plantCAccountingDelta = accounting; Envi before = envi; @@ -170,7 +154,6 @@ static void partialCase(int accounting, const char *harvest, double fraction) { processEvents(); updateBalanceTrackerPreUpdate(); updatePoolsForEvents(); - // near(updatePoolsForFullHarvest(), 0, "partial harvest not deferred"); near(envi.plantWoodC, before.plantWoodC * (1 - fraction), "partial wood"); near(envi.plantCAccountingDelta, accounting * (1 - fraction), "partial accounting"); @@ -183,7 +166,7 @@ static void partialCase(int accounting, const char *harvest, double fraction) { } static void partialTimestepCase(int accounting) { - logTest("Running partial timestep case\n"); + logTest("*** Running partial timestep case\n"); start(2, ".1 .2 .3 .4"); envi.plantCAccountingDelta = accounting; updateState(); @@ -206,32 +189,43 @@ static void invalidCase(int which) { start(2, "0 0 1 1"); if (which < 4) { // 0, 1, 2, 3 const char *bad[] = {"0 0 -1 1", "-0.1 0 1 1", "0 0 1.1 1", "1 0.1 0 1"}; + // createEventNode should reject these createEventNode(2017, 20, HARVEST, bad[which]); - } else { // 4, 5 + } else { // 4, 5, 6 EventNode *extra = createEventNode(2017, 20, HARVEST, "0 0 .2 .2"); - if (which < 5) + if (which == 4) { gEvents->nextEvent = extra; - else { + // Harvest greater than 100% + processEvents(); + } else if (which == 5) { extra->nextEvent = gEvents; gEvents = extra; setupEvents(); + // Harvest greater than 100% + processEvents(); + } else { // 6 + EventNode *leafOff1 = createEventNode(2017, 20, LEAFOFF, ""); + EventNode *leafOff2 = createEventNode(2017, 20, LEAFOFF, ""); + leafOff1->nextEvent = leafOff2; + gEvents = leafOff1; + setupEvents(); + processEvents(); } - processEvents(); } _exit(99); } int status; waitpid(child, &status, 0); - int expected = - which < 6 ? EXIT_CODE_BAD_PARAMETER_VALUE : EXIT_CODE_INTERNAL_ERROR; + int expected = EXIT_CODE_BAD_PARAMETER_VALUE; if (!WIFEXITED(status) || WEXITSTATUS(status) != expected) { - logTest("Invalid case %d did not exit with code %d\n", which, expected); + logTest("Invalid case %d did not exit with code %d (was %d)\n", which, + expected, WIFEXITED(status) ? WEXITSTATUS(status) : -1); failures++; } } static void coincidentFertilizerCase(void) { - logTest("Running coincident fertilizer case\n"); + logTest("*** Running coincident fertilizer case\n"); Envi expected = {0}; for (int h = 0; h < 2; h++) { start(2, h ? "0 0 1 1" : NULL); @@ -276,7 +270,7 @@ static void sameEnvi(const Envi *a, const Envi *b) { #undef E } static void coincidentEventCase(int type, const char *arguments) { - logTest("Running coincident event case with type: %s params: %s\n", + logTest("*** Running coincident event case with type: %s params: %s\n", eventTypeToString(type), *arguments ? arguments : ""); Envi end; Fluxes ordinary; @@ -344,7 +338,7 @@ static void coincidentEventCase(int type, const char *arguments) { } static void restartCase(void) { - logTest("Running restart case\n"); + logTest("*** Running restart case\n"); Envi expected; for (int resumed = 0; resumed < 2; resumed++) { start(2, "0 0 1 1"); @@ -388,7 +382,7 @@ static void leafBudgetCases(void) { int caseNum = 0; for (int kind = 0; kind < 4; kind++) { // kind: 0 1 2 3 - // fracLeafFall 1 .25 .75 .75 + // fracLeafFall 1 .25 .40 .75 // leafOff a a b c // a: inserted at start // b: two events inserted at start @@ -407,12 +401,10 @@ static void leafBudgetCases(void) { // TODO: Add "1 1 0 0" and "0.5 0.5 0.5 0.5" cases for harvest logTest( "*** Running leaf budget case [%d] with kind: %d dark: %d " - "sign %d" - " account %d limited %d resorb %.1f harvest %s\n", + "sign %d account %d limited %d resorb %.1f harvest %s\n", caseNum++, kind, dark, sign, account, limited, resorb / 2.0, harvest ? "0 0 1 1" : "none"); start(2, harvest ? "0 0 1 1" : NULL); - logInfo("* Post start setup *\n"); envi.plantLeafC = 1; envi.plantWoodC = envi.fineRootC = envi.coarseRootC = 100; envi.plantCAccountingDelta = account; @@ -428,7 +420,10 @@ static void leafBudgetCases(void) { params.nVolatilizationFrac = params.nLeachingFrac = 0; params.nFixationFracMax = 0; params.leafNResorptionFrac = .5 * resorb; - params.fracLeafFall = kind == 0 ? 1 : kind == 1 ? .25 : .75; + params.fracLeafFall = kind == 0 ? 1 + : kind == 1 ? .25 + : kind == 2 ? .40 + : .75; if (dark) climate->par = 0; resetMeanTracker(meanNPP, 20 * sign); @@ -462,7 +457,7 @@ static void leafBudgetCases(void) { double expLeafOffLitter = kind == 3 ? 0 - : fmin(params.fracLeafFall, + : fmin(params.fracLeafFall * (1 + (kind == 2)), beforeLeafC + (fluxes.leafCreation - fluxes.leafLitter) * climate->length); @@ -474,8 +469,6 @@ static void leafBudgetCases(void) { near(envi.litterN - beforeLitterN, shed * (1 - params.leafNResorptionFrac) / params.leafCN, "leaf litter N transfer"); - logTest("beforeLitterN %f afterLitterN %f\n", beforeLitterN, - envi.litterN); } else { clearPlant(); climate = climate->nextClim; @@ -486,23 +479,10 @@ static void leafBudgetCases(void) { "no regrowth after complete harvest"); } finish(); - - logInfo(":\n"); - if (failures > 0) { - logTest("Complete harvest failures: %d\n", failures); - logTest("Press return to continue\n"); - getchar(); - } } // harvest loop } // resorb loop } // limited loop } // account loop - if (kind > 1) { - logTest("New sign value incoming\n"); - logTest("Complete harvest failures: %d\n", failures); - logTest("Press return to continue\n"); - getchar(); - } } // sign loop } // dark loop } // kind loop @@ -525,42 +505,42 @@ static void exportLossCase(void) { int main(void) { const char *harvest[] = {"0 0 1 1", "1 1 0 0", ".25 .75 .75 .25"}; const double above[] = {0, 1, .25}, below[] = {0, 1, .75}; - // for (int mode = 0; mode < 3; mode++) - // for (int account = -1; account <= 1; account++) - // for (int dark = 0; dark < 2; dark++) - // for (int limited = 0; limited <= (mode == 2); limited++) - // for (int route = 0; route < 3; route++) - // fullCase(mode, account, dark, limited, harvest[route], - // above[route], - // below[route]); - - fullCase(0, -1, 0, 0, harvest[0], above[0], below[0]); + for (int mode = 0; mode < 3; mode++) + for (int account = -1; account <= 1; account++) + for (int dark = 0; dark < 2; dark++) + for (int limited = 0; limited <= (mode == 2); limited++) + for (int route = 0; route < 3; route++) + fullCase(mode, account, dark, limited, harvest[route], above[route], + below[route]); + logTest("\n"); for (int account = -1; account <= 1; account++) { partialTimestepCase(account); partialCase(account, ".1 .2 .3 .4", .4); partialCase(account, "0 0 .999999999 .999999999", .999999999); } - logTest("Running invalid case checks; six errors expected\n"); - for (int which = 0; which < 6; which++) { - invalidCase(which); - } - + logTest("\n"); coincidentFertilizerCase(); coincidentEventCase(LEAFON, ""); coincidentEventCase(LEAFOFF, ""); coincidentEventCase(IRRIGATION, "2 0"); restartCase(); - logTest("\n\n\n"); + logTest("\n"); leafBudgetCases(); + exportLossCase(); - logTest("Press return to continue\n"); - getchar(); + // Leave these last, as the forking messes up the output otherwise + logTest("\n"); + logTest("*** Running invalid case checks; seven errors expected\n"); + for (int which = 0; which < 7; which++) { + invalidCase(which); + } + + logTest("\n"); + logTest("Complete harvest total failures: %d\n", failures); - exportLossCase(); - logTest("Complete harvest failures: %d\n", failures); return failures != 0; } From b691a19188007b71180ebf5beab469729d169b63 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Tue, 29 Sep 2026 12:37:37 -0400 Subject: [PATCH 07/14] All new test points passing --- tests/sipnet/test_modeling/testCompleteHarvest.c | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/sipnet/test_modeling/testCompleteHarvest.c b/tests/sipnet/test_modeling/testCompleteHarvest.c index d3d3713a..2fce3b54 100644 --- a/tests/sipnet/test_modeling/testCompleteHarvest.c +++ b/tests/sipnet/test_modeling/testCompleteHarvest.c @@ -380,6 +380,7 @@ static void restartCase(void) { static void leafBudgetCases(void) { int caseNum = 0; + char harvests[4][16] = {"none", "0 0 1 1", "0.5 0.5 0.5 0.5", "1 1 0 0"}; for (int kind = 0; kind < 4; kind++) { // kind: 0 1 2 3 // fracLeafFall 1 .25 .40 .75 @@ -397,14 +398,14 @@ static void leafBudgetCases(void) { // limited: starting minN = 0 if limited (1000 else) for (int resorb = 0; resorb < 3; resorb++) { /// resorb: leafNResorptionFrac = 0, .5, 1 (resorb/2) - for (int harvest = 0; harvest < 2; harvest++) { - // TODO: Add "1 1 0 0" and "0.5 0.5 0.5 0.5" cases for harvest + for (int harvest = 0; harvest < 4; harvest++) { + char *harvestStr = harvests[harvest]; logTest( "*** Running leaf budget case [%d] with kind: %d dark: %d " "sign %d account %d limited %d resorb %.1f harvest %s\n", caseNum++, kind, dark, sign, account, limited, resorb / 2.0, - harvest ? "0 0 1 1" : "none"); - start(2, harvest ? "0 0 1 1" : NULL); + harvestStr); + start(2, harvest ? harvestStr : NULL); envi.plantLeafC = 1; envi.plantWoodC = envi.fineRootC = envi.coarseRootC = 100; envi.plantCAccountingDelta = account; From 42acc96922f97d5d66f2411a8ea50febf8795404 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Tue, 29 Sep 2026 15:29:56 -0400 Subject: [PATCH 08/14] Put events guard around writeComputedEvents --- src/sipnet/sipnet.c | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 068e3259..000a63fd 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -819,10 +819,12 @@ void calcLeafOnOffFluxes(void) { } phenologyTrackers.didLeafGrowth = 1; - writeComputedEventOut(climate->year, climate->day, - eventTypeToString(LEAFON), 2, "leafOnCreation", - leafOn * len, "leafOnCreationFromWood", - leafOnFromWood * len); + if (ctx.events) { + writeComputedEventOut(climate->year, climate->day, + eventTypeToString(LEAFON), 2, "leafOnCreation", + leafOn * len, "leafOnCreationFromWood", + leafOnFromWood * len); + } } // check for end of growing season: From 9a5a9e25a93a00eff7c4a8d02e197a2a23e98a27 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Wed, 30 Sep 2026 10:54:11 -0400 Subject: [PATCH 09/14] Update for C/N param ordering --- tests/smoke/russell_2/events.out | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/smoke/russell_2/events.out b/tests/smoke/russell_2/events.out index 45b9422a..f387a652 100644 --- a/tests/smoke/russell_2/events.out +++ b/tests/smoke/russell_2/events.out @@ -1,11 +1,11 @@ year day type param_name=delta[,param_name=delta,...] 2016 47 leafon leafOnCreation=114.61,leafOnCreationFromWood=86.11 -2016 90 fert eventLitterC=5.00,eventSoilC=0.00,eventMinN=10.00,eventLitterN=15.00,eventInputC=5.00,eventInputN=25.00 +2016 90 fert eventLitterC=5.00,eventSoilC=0.00,eventInputC=5.00,eventMinN=10.00,eventLitterN=15.00,eventInputN=25.00 2016 92 irrig eventSoilWater=2.80,eventEvap=0.00 2016 96 irrig eventSoilWater=2.80,eventEvap=0.00 2016 100 irrig eventSoilWater=2.80,eventEvap=0.00 2016 104 irrig eventSoilWater=2.80,eventEvap=0.00 -2016 104 fert eventLitterC=5.00,eventSoilC=0.00,eventMinN=10.00,eventLitterN=15.00,eventInputC=5.00,eventInputN=25.00 +2016 104 fert eventLitterC=5.00,eventSoilC=0.00,eventInputC=5.00,eventMinN=10.00,eventLitterN=15.00,eventInputN=25.00 2016 108 irrig eventSoilWater=2.80,eventEvap=0.00 2016 112 irrig eventSoilWater=2.80,eventEvap=0.00 2016 116 irrig eventSoilWater=2.80,eventEvap=0.00 From a4eb1a964e96847a56aa40367a0eb0cc5685fed3 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Wed, 30 Sep 2026 10:54:21 -0400 Subject: [PATCH 10/14] Update event logging --- src/sipnet/events.c | 37 ++++++++++++++++++++++--------------- src/sipnet/sipnet.c | 2 +- 2 files changed, 23 insertions(+), 16 deletions(-) diff --git a/src/sipnet/events.c b/src/sipnet/events.c index 1e1a938b..80db2e21 100644 --- a/src/sipnet/events.c +++ b/src/sipnet/events.c @@ -390,9 +390,7 @@ void appendLog(EventNode *event, int numParams, ...) { for (int ind = 0; ind < numParams; ind++) { char *param = va_arg(args, char *); double val = va_arg(args, double); - char suffix = (ind == numParams - 1) ? '\n' : ','; - int success = - dsAppendFormatted(event->logLine, "%s=%-.2f%c", param, val, suffix); + int success = dsAppendFormatted(event->logLine, "%s=%-.2f,", param, val); if (!success) { logError("appending event log line failed\n"); exit(EXIT_CODE_INTERNAL_ERROR); @@ -407,8 +405,17 @@ void writeEventsOut(void) { const int climYear = climate->year; const int climDay = climate->day; while (gEvent != NULL && gEvent->year <= climYear && gEvent->day <= climDay) { + // Change last char to a newline if it is a comma + // Get the current length of the string + char *log = gEvent->logLine->buffer; + size_t len = strlen(log); + if (len > 0 && log[len - 1] == ',') { + log[len - 1] = '\n'; + } else { + dsAppend(gEvent->logLine, "\n"); + } fprintf(eventOutFile, "%4d %3d %-7s %s", gEvent->year, gEvent->day, - eventTypeToString(gEvent->type), gEvent->logLine->buffer); + eventTypeToString(gEvent->type), log); gEvent = gEvent->nextEvent; } } @@ -520,7 +527,7 @@ void processEventsForCarbon(EventNode *event) { fluxes.eventEvap += evapAmount / climLen; fluxes.eventSoilWater += soilAmount / climLen; - appendLog(gEvent, 2, "eventSoilWater", soilAmount, "eventEvap", + appendLog(event, 2, "eventSoilWater", soilAmount, "eventEvap", evapAmount); } break; @@ -545,7 +552,7 @@ void processEventsForCarbon(EventNode *event) { fluxes.eventInputC += inputC / climLen; // clang-format off - appendLog(gEvent, 5, + appendLog(event, 5, "eventLeafC", leafC, "eventWoodC", woodC, "eventFineRootC", fineRootC, @@ -631,7 +638,7 @@ void processEventsForCarbon(EventNode *event) { // clang-format off appendLog( - gEvent, 8, + event, 8, "eventSoilC", soilAdd, "eventLitterC", litterAdd, "eventLeafC", leafDelta, @@ -652,7 +659,7 @@ void processEventsForCarbon(EventNode *event) { // as there may be lingering effects from a prior tillage. eventTrackers.d_till_mod += tillParams->tillageEffect; - appendLog(gEvent, 1, "eventTrackers.d_till_mod", + appendLog(event, 1, "eventTrackers.d_till_mod", tillParams->tillageEffect); } break; @@ -669,7 +676,7 @@ void processEventsForCarbon(EventNode *event) { fluxes.eventInputC += orgC / climLen; // clang-format off - appendLog(gEvent, 3, + appendLog(event, 3, "eventLitterC", ctx.litterPool ? orgC : 0.0, "eventSoilC", ctx.litterPool ? 0.0 : orgC, "eventInputC", orgC); @@ -690,7 +697,7 @@ void processEventsForCarbon(EventNode *event) { // eventInputC // clang-format off - appendLog(gEvent, 2, + appendLog(event, 2, "eventLeafOnCreation", leafOnFlux * climLen, "eventLeafOnCreationFromWood", leafOnFluxFromWood * climLen); // clang-format on @@ -708,7 +715,7 @@ void processEventsForCarbon(EventNode *event) { double leafOff = envi.plantLeafC * params.fracLeafFall; fluxes.eventLeafOffLitterC += leafOff / climLen; - appendLog(gEvent, 1, "eventLeafOffLitter", leafOff); + appendLog(event, 1, "eventLeafOffLitter", leafOff); } break; case PLANTDEATH: // There should be no way to get here, but covering our bases... @@ -767,7 +774,7 @@ void processEventsForNitrogen(EventNode *event) { coarseRootC / params.woodCN; fluxes.eventInputN += inputN / climLen; - appendLog(gEvent, 1, "eventInputN", inputN); + appendLog(event, 1, "eventInputN", inputN); } break; case HARVEST: { // Harvest can both remove biomass and move biomass to the soil/litter @@ -803,7 +810,7 @@ void processEventsForNitrogen(EventNode *event) { // clang-format off appendLog( - gEvent, 3, + event, 3, "eventSoilOrgN", soilNAdd, "eventLitterN", litterNAdd, "eventOutputN", outputN); @@ -826,7 +833,7 @@ void processEventsForNitrogen(EventNode *event) { fluxes.eventInputN += (orgN + minN) / climLen; // clang-format off - appendLog(gEvent, 3, + appendLog(event, 3, "eventMinN", minN, "eventLitterN", orgN, "eventInputN", (orgN + minN)); @@ -866,7 +873,7 @@ void processEventsForNitrogen(EventNode *event) { double litterNAddFlux = fluxes.eventLeafOffLitterN - preLitter; // clang-format off - appendLog(gEvent, 2, + appendLog(event, 2, "eventLeafOffNResorption", leafNResorptionFlux * climLen, "eventLitterN", litterNAddFlux * climLen); // clang-format on diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 000a63fd..e24b4717 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -819,7 +819,7 @@ void calcLeafOnOffFluxes(void) { } phenologyTrackers.didLeafGrowth = 1; - if (ctx.events) { + if (ctx.events && leafOn > TINY) { writeComputedEventOut(climate->year, climate->day, eventTypeToString(LEAFON), 2, "leafOnCreation", leafOn * len, "leafOnCreationFromWood", From b2200fd7ce4dfb7d993465280d02a575041fe400 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Fri, 2 Oct 2026 13:34:01 -0400 Subject: [PATCH 11/14] Update for current fluxes/trackers --- src/sipnet/debug_log.c | 92 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 80 insertions(+), 12 deletions(-) diff --git a/src/sipnet/debug_log.c b/src/sipnet/debug_log.c index 1c62d77a..6f3119a8 100644 --- a/src/sipnet/debug_log.c +++ b/src/sipnet/debug_log.c @@ -4,11 +4,44 @@ #include "debug_log.h" +#include "events.h" #include "common/exitCodes.h" #include "common/logging.h" #include "common/context.h" #include "common/util.h" #include "state.h" + +#define NUM_LOGGED_ENVI_FIELDS 13 +#define NUM_LOGGED_FLUX_FIELDS 59 +#define NUM_LOGGED_TRACKER_FIELDS 33 +#define NUM_LOGGED_PHEN_TRACKER_FIELDS 3 +#define NUM_LOGGED_SURVIVAL_FIELDS 1 +#define NUM_LOGGED_EVENT_TRACKER_FIELDS 7 + +#define DEBUG_LAYOUT_ENVI_SIZE (8 * NUM_LOGGED_ENVI_FIELDS) +#define DEBUG_LAYOUT_FLUX_SIZE (8 * NUM_LOGGED_FLUX_FIELDS) +// The Trackers struct is not all doubles, but the int(s) get padded to 8 bytes +#define DEBUG_LAYOUT_TRACKER_SIZE (8 * NUM_LOGGED_TRACKER_FIELDS) +#define DEBUG_LAYOUT_PHEN_SIZE (4 * NUM_LOGGED_PHEN_TRACKER_FIELDS) +#define DEBUG_LAYOUT_SURVIVAL_SIZE (4 * NUM_LOGGED_SURVIVAL_FIELDS) +#define DEBUG_LAYOUT_EVENT_SIZE (8 * NUM_LOGGED_EVENT_TRACKER_FIELDS) + +_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"); +_Static_assert( + sizeof(PhenologyTrackers) == DEBUG_LAYOUT_PHEN_SIZE, + "Debug log schema drift: PhenologyTrackers changed; update debug_log.c"); +_Static_assert( + 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 { DEBUG_FIELD_INT = 0, DEBUG_FIELD_DOUBLE = 1 @@ -20,18 +53,13 @@ typedef struct DebugField { const void *value; } DebugField; -#define NUM_LOGGED_ENVI_FIELDS 13 -#define NUM_LOGGED_FLUX_FIELDS 57 -#define NUM_LOGGED_TRACKER_FIELDS 33 -#define NUM_LOGGED_PHEN_TRACKER_FIELDS 3 -#define NUM_LOGGED_SURVIVAL_FIELDS 1 - typedef struct DebugFieldArrays { DebugField enviDF[NUM_LOGGED_ENVI_FIELDS]; DebugField fluxDF[NUM_LOGGED_FLUX_FIELDS]; DebugField trackerDF[NUM_LOGGED_TRACKER_FIELDS]; DebugField phenoDF[NUM_LOGGED_PHEN_TRACKER_FIELDS]; DebugField survivalDF[NUM_LOGGED_SURVIVAL_FIELDS]; + DebugField eventDF[NUM_LOGGED_EVENT_TRACKER_FIELDS]; } DebugFieldArrays; static DebugFieldArrays *debugFields = NULL; @@ -61,7 +89,11 @@ void initDebugArrays() { debugFields->enviDF[ind++] = (DebugField){"soilOrgN", DEBUG_FIELD_DOUBLE, &envi.soilOrgN}, debugFields->enviDF[ind++] = (DebugField){"litterN", DEBUG_FIELD_DOUBLE, &envi.litterN}, debugFields->enviDF[ind++] = (DebugField){"plantStorageN", DEBUG_FIELD_DOUBLE, &envi.plantStorageN}, - debugFields->enviDF[ind ] = (DebugField){"plantCAccountingDelta", DEBUG_FIELD_DOUBLE,&envi.plantCAccountingDelta}; + debugFields->enviDF[ind++] = (DebugField){"plantCAccountingDelta", DEBUG_FIELD_DOUBLE,&envi.plantCAccountingDelta}; + if (ind != NUM_LOGGED_ENVI_FIELDS) { + logInternalError("Debug log array size mismatch: enviDF\n"); + exit(EXIT_CODE_INTERNAL_ERROR); + } ind = 0; debugFields->fluxDF[ind++] = (DebugField){"photosynthesis", DEBUG_FIELD_DOUBLE, &fluxes.photosynthesis}; @@ -90,6 +122,7 @@ void initDebugArrays() { 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}, debugFields->fluxDF[ind++] = (DebugField){"nOrgSoil", DEBUG_FIELD_DOUBLE, &fluxes.nOrgSoil}, @@ -117,10 +150,15 @@ void initDebugArrays() { debugFields->fluxDF[ind++] = (DebugField){"eventOutputN", DEBUG_FIELD_DOUBLE, &fluxes.eventOutputN}, debugFields->fluxDF[ind++] = (DebugField){"eventLeafOnCreation", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOnCreation}, debugFields->fluxDF[ind++] = (DebugField){"eventLeafOnCreationFromWood", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOnCreationFromWood}, - debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitter", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitterC}, + debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitterC", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitterC}, + debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitterN", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitterN}, debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffNResorption", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffNResorption}, debugFields->fluxDF[ind++] = (DebugField){"soilMethane", DEBUG_FIELD_DOUBLE, &fluxes.soilMethane}, - debugFields->fluxDF[ind ] = (DebugField){"litterMethane", DEBUG_FIELD_DOUBLE, &fluxes.litterMethane}; + debugFields->fluxDF[ind++] = (DebugField){"litterMethane", DEBUG_FIELD_DOUBLE, &fluxes.litterMethane}; + if (ind != NUM_LOGGED_FLUX_FIELDS) { + logInternalError("Debug log array size mismatch: fluxDF\n"); + exit(EXIT_CODE_INTERNAL_ERROR); + } ind = 0; debugFields->trackerDF[ind++] = (DebugField){"gpp", DEBUG_FIELD_DOUBLE, &trackers.gpp}, @@ -155,15 +193,40 @@ void initDebugArrays() { debugFields->trackerDF[ind++] = (DebugField){"nLeaching", DEBUG_FIELD_DOUBLE, &trackers.nLeaching}, debugFields->trackerDF[ind++] = (DebugField){"nFixation", DEBUG_FIELD_DOUBLE, &trackers.nFixation}, debugFields->trackerDF[ind++] = (DebugField){"nUptake", DEBUG_FIELD_DOUBLE, &trackers.nUptake}; - debugFields->trackerDF[ind ] = (DebugField){"meanNPP", DEBUG_FIELD_DOUBLE, &trackers.meanNPP}; + debugFields->trackerDF[ind++] = (DebugField){"meanNPP", DEBUG_FIELD_DOUBLE, &trackers.meanNPP}; + if (ind != NUM_LOGGED_TRACKER_FIELDS) { + logInternalError("Debug log array size mismatch: trackerDF\n"); + exit(EXIT_CODE_INTERNAL_ERROR); + } ind = 0; debugFields->phenoDF[ind++] = (DebugField){"didLeafGrowth", DEBUG_FIELD_INT, &phenologyTrackers.didLeafGrowth}, debugFields->phenoDF[ind++] = (DebugField){"didLeafFall", DEBUG_FIELD_INT, &phenologyTrackers.didLeafFall}, - debugFields->phenoDF[ind ] = (DebugField){"lastYear", DEBUG_FIELD_INT, &phenologyTrackers.lastYear}; + debugFields->phenoDF[ind++] = (DebugField){"lastYear", DEBUG_FIELD_INT, &phenologyTrackers.lastYear}; + if (ind != NUM_LOGGED_PHEN_TRACKER_FIELDS) { + logInternalError("Debug log array size mismatch: phenoDF\n"); + exit(EXIT_CODE_INTERNAL_ERROR); + } ind = 0; - debugFields->survivalDF[ind] = (DebugField){"isAlive", DEBUG_FIELD_INT, &plantSurvivalTracker.isAlive}; + debugFields->survivalDF[ind++] = (DebugField){"isAlive", DEBUG_FIELD_INT, &plantSurvivalTracker.isAlive}; + if (ind != NUM_LOGGED_SURVIVAL_FIELDS) { + logInternalError("Debug log array size mismatch: survivalDF\n"); + exit(EXIT_CODE_INTERNAL_ERROR); + } + + ind = 0; + debugFields->eventDF[ind++] = (DebugField){"d_till_mod", DEBUG_FIELD_DOUBLE, &eventTrackers.d_till_mod}; + debugFields->eventDF[ind++] = (DebugField){"ht.totalFracRemoved", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemoved}; + debugFields->eventDF[ind++] = (DebugField){"ht.totalFracTransferred", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferred}; + debugFields->eventDF[ind++] = (DebugField){"ht.totalFracRemovedAbove", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemovedAbove}; + debugFields->eventDF[ind++] = (DebugField){"ht.totalFracRemovedBelow", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemovedBelow}; + debugFields->eventDF[ind++] = (DebugField){"ht.totalFracTransferredAbove", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferredAbove}; + debugFields->eventDF[ind++] = (DebugField){"ht.totalFracTransferredBelow", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferredBelow}; + if (ind != NUM_LOGGED_EVENT_TRACKER_FIELDS) { + logInternalError("Debug log array size mismatch: eventDF\n"); + exit(EXIT_CODE_INTERNAL_ERROR); + } // clang-format on } @@ -279,6 +342,8 @@ void outputDebugHeaders(DebugLogFiles *debugLogFiles) { outputDebugFieldHeader(debugLogFiles->trackers, "s.", debugFields->survivalDF, NUM_LOGGED_SURVIVAL_FIELDS, 0); + outputDebugFieldHeader(debugLogFiles->trackers, "et.", debugFields->eventDF, + NUM_LOGGED_EVENT_TRACKER_FIELDS, 0); fprintf(debugLogFiles->trackers, "\n"); } } @@ -308,6 +373,9 @@ void outputDebugState(DebugLogFiles *debugLogFiles, int year, int day, outputDebugFieldValues(debugLogFiles->trackers, year, day, time, debugFields->survivalDF, NUM_LOGGED_SURVIVAL_FIELDS, 0); + outputDebugFieldValues(debugLogFiles->trackers, year, day, time, + debugFields->eventDF, + NUM_LOGGED_EVENT_TRACKER_FIELDS, 0); fprintf(debugLogFiles->trackers, "\n"); } } From a47cc0d7627f0be414e8b28f3bce32638e3d4a53 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Fri, 2 Oct 2026 13:34:49 -0400 Subject: [PATCH 12/14] Minor doc update --- src/sipnet/events.c | 5 +++-- src/sipnet/limitations.c | 3 ++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/sipnet/events.c b/src/sipnet/events.c index 80db2e21..6b7f9765 100644 --- a/src/sipnet/events.c +++ b/src/sipnet/events.c @@ -500,8 +500,9 @@ void processEventsForCarbon(EventNode *event) { // did not have a corresponding climate file record. if (event->year < climYear || event->day < climDay) { logError("Agronomic event found for year: %d day: %d that does not " - "have a corresponding record in the climate file\n", - event->year, event->day); + "have a corresponding record in the climate file " + "(current climate date is year %d day %d)\n", + event->year, event->day, climYear, climDay); exit(EXIT_CODE_INPUT_FILE_ERROR); } diff --git a/src/sipnet/limitations.c b/src/sipnet/limitations.c index a074c463..6a5d3a9e 100644 --- a/src/sipnet/limitations.c +++ b/src/sipnet/limitations.c @@ -185,12 +185,13 @@ static void checkNitrogenLimitation(void) { checkNegativeCreation(); // Reset and recalc event leaf-off N + // Nitrogen for "calculated" leaf-off events will be redone in + // calcNitrogenFluxes below fluxes.eventLeafOffNResorption = 0.0; fluxes.eventLeafOffLitterN = 0.0; calcLeafOffNEffects(fluxes.eventLeafOffLitterC * len, &fluxes.eventLeafOffNResorption, &fluxes.eventLeafOffLitterN); - // } // Reset N calculations calcNitrogenFluxes(); From 245e86c0a76d66d3d57cc4ae2dd9f97b6f077bad Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Fri, 2 Oct 2026 13:35:11 -0400 Subject: [PATCH 13/14] Update for harvest/leaf-off fixes --- .../events_output_header.out | 8 ++++---- .../events_output_no_header.out | 8 ++++---- .../sipnet/test_modeling/testCompleteHarvest.c | 13 +++++++++++++ tests/sipnet/test_modeling/testNitrogenCycle.c | 18 ++++++++++++++---- .../sipnet/test_modeling/testPlantMortality.c | 4 ++-- .../testDebugLogFiles.c | 13 ++++++------- tests/utils/helpers.c | 5 ++++- 7 files changed, 47 insertions(+), 22 deletions(-) diff --git a/tests/sipnet/test_events_infrastructure/events_output_header.out b/tests/sipnet/test_events_infrastructure/events_output_header.out index d00de444..0bda11a3 100644 --- a/tests/sipnet/test_events_infrastructure/events_output_header.out +++ b/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 eventLeafC=3.00,eventWoodC=4.00,eventFineRootC=5.00,eventCoarseRootC=6.00,eventInputC=18.00,eventInputN=0.00 +2023 65 plant eventLeafC=3.00,eventWoodC=4.00,eventFineRootC=5.00,eventCoarseRootC=6.00,eventInputC=18.00 2023 70 irrig eventSoilWater=5.00,eventEvap=0.00 -2023 200 harv eventSoilC=1.90,eventLitterC=3.56,eventLeafC=-5.93,eventWoodC=-4.75,eventAccountingC=-0.00,eventFineRootC=-3.73,eventCoarseRootC=-3.89,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=12.83,eventOutputN=0.00 -2024 65 plant eventLeafC=3.00,eventWoodC=5.00,eventFineRootC=7.00,eventCoarseRootC=9.00,eventInputC=24.00,eventInputN=0.00 +2023 200 harv eventSoilC=1.90,eventLitterC=3.56,eventLeafC=-5.93,eventWoodC=-4.75,eventAccountingC=-0.00,eventFineRootC=-3.73,eventCoarseRootC=-3.89,eventOutputC=12.83 +2024 65 plant eventLeafC=3.00,eventWoodC=5.00,eventFineRootC=7.00,eventCoarseRootC=9.00,eventInputC=24.00 2024 70 irrig eventSoilWater=2.50,eventEvap=2.50 -2024 200 harv eventSoilC=2.74,eventLitterC=1.51,eventLeafC=-1.39,eventWoodC=-1.63,eventAccountingC=-0.00,eventFineRootC=-2.52,eventCoarseRootC=-2.97,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=4.25,eventOutputN=0.00 +2024 200 harv eventSoilC=2.74,eventLitterC=1.51,eventLeafC=-1.39,eventWoodC=-1.63,eventAccountingC=-0.00,eventFineRootC=-2.52,eventCoarseRootC=-2.97,eventOutputC=4.25 diff --git a/tests/sipnet/test_events_infrastructure/events_output_no_header.out b/tests/sipnet/test_events_infrastructure/events_output_no_header.out index 3c727e29..553d7073 100644 --- a/tests/sipnet/test_events_infrastructure/events_output_no_header.out +++ b/tests/sipnet/test_events_infrastructure/events_output_no_header.out @@ -1,6 +1,6 @@ -2023 65 plant eventLeafC=10.00,eventWoodC=5.00,eventFineRootC=4.00,eventCoarseRootC=3.00,eventInputC=22.00,eventInputN=0.00 +2023 65 plant eventLeafC=10.00,eventWoodC=5.00,eventFineRootC=4.00,eventCoarseRootC=3.00,eventInputC=22.00 2023 70 irrig eventSoilWater=5.00,eventEvap=0.00 -2023 200 harv eventSoilC=12.40,eventLitterC=0.00,eventLeafC=-4.80,eventWoodC=-3.20,eventAccountingC=-0.00,eventFineRootC=-4.80,eventCoarseRootC=-4.80,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=5.20,eventOutputN=0.00 -2024 65 plant eventLeafC=10.00,eventWoodC=5.00,eventFineRootC=4.00,eventCoarseRootC=3.00,eventInputC=22.00,eventInputN=0.00 +2023 200 harv eventSoilC=12.40,eventLitterC=0.00,eventLeafC=-4.80,eventWoodC=-3.20,eventAccountingC=-0.00,eventFineRootC=-4.80,eventCoarseRootC=-4.80,eventOutputC=5.20 +2024 65 plant eventLeafC=10.00,eventWoodC=5.00,eventFineRootC=4.00,eventCoarseRootC=3.00,eventInputC=22.00 2024 70 irrig eventSoilWater=2.50,eventEvap=2.50 -2024 200 harv eventSoilC=12.14,eventLitterC=0.00,eventLeafC=-10.32,eventWoodC=-5.88,eventAccountingC=-0.00,eventFineRootC=-2.88,eventCoarseRootC=-2.48,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=9.42,eventOutputN=0.00 +2024 200 harv eventSoilC=12.14,eventLitterC=0.00,eventLeafC=-10.32,eventWoodC=-5.88,eventAccountingC=-0.00,eventFineRootC=-2.88,eventCoarseRootC=-2.48,eventOutputC=9.42 diff --git a/tests/sipnet/test_modeling/testCompleteHarvest.c b/tests/sipnet/test_modeling/testCompleteHarvest.c index 2fce3b54..d7c3db82 100644 --- a/tests/sipnet/test_modeling/testCompleteHarvest.c +++ b/tests/sipnet/test_modeling/testCompleteHarvest.c @@ -504,6 +504,9 @@ static void exportLossCase(void) { } int main(void) { + + logTest("Starting testCompleteHarvest\n"); + const char *harvest[] = {"0 0 1 1", "1 1 0 0", ".25 .75 .75 .25"}; const double above[] = {0, 1, .25}, below[] = {0, 1, .75}; for (int mode = 0; mode < 3; mode++) @@ -542,6 +545,16 @@ int main(void) { logTest("\n"); logTest("Complete harvest total failures: %d\n", failures); + logTest("\n"); + + int status = failures > 0; + + if (status) { + logTest("FAILED testCompleteHarvest with status %d\n", status); + exit(status); + } + + logTest("PASSED testCompleteHarvest\n"); return failures != 0; } diff --git a/tests/sipnet/test_modeling/testNitrogenCycle.c b/tests/sipnet/test_modeling/testNitrogenCycle.c index 95e419c0..1e061e51 100644 --- a/tests/sipnet/test_modeling/testNitrogenCycle.c +++ b/tests/sipnet/test_modeling/testNitrogenCycle.c @@ -2,6 +2,7 @@ #include "sipnet/events.c" #include "sipnet/nitrogen.c" #include "sipnet/limitations.c" +#include "utils/helpers.c" ///// // Setup and general test state management @@ -134,15 +135,17 @@ int testFertilization(void) { // init minN 2, nVol 0.1 initNVolatilizationState(initN, nVolFrac); - // fert event: 15 5 10 + // fert event: 15 orgN 5 orgC 10 minN double fertMinN = 10; initEvents("events_fert.in", "events.out", 0); setupEvents(); + EventNode *event = getCurrentEvent(); + processEventsForCarbon(event); + processEventsForNitrogen(event); calcNVolatilizationFlux(); - processEvents(); - updateNitrogenPools(); updatePoolsForEvents(); + updateNitrogenPools(); // Want to test: // envi: minN @@ -708,6 +711,8 @@ int testLeafTurnoverNResorption(void) { envi.minN = 1.0; envi.plantStorageN = 0.0; fluxes.leafOffNResorption = 2.0; + params.leafNResorptionFrac = 0.5; + updateNitrogenPools(); double expStorageN = 2.0 * climate->length; // 0.25 @@ -727,9 +732,14 @@ int testLeafTurnoverNResorption(void) { double availableN = minN2; // plus unclaimed plantStorageN, which is 0 here double maxUptake = demandFlux * climate->length; double reduction = availableN / maxUptake; + // enough so we don't limit leaf litter in checkNegativeCreation + envi.plantLeafC = 8.0; initNLimitationState(minN2, 0); fluxes.leafOffNResorption = resorpFlux; + // to make sure leafOffNResorption is not reduced by leafOffLitter + fluxes.leafOffLitter = + resorpFlux * params.leafCN / params.leafNResorptionFrac; doNFixUpLimitCalcs(); updateNitrogenPools(); @@ -738,7 +748,7 @@ int testLeafTurnoverNResorption(void) { "[turnover resorption] leafCreation"); status |= checkNLimitationFlux(fluxes.woodCreation, 500 * reduction, "[turnover resorption] woodCreation"); - status |= checkMinAndStorageN("turnover resorption", 0.0, + status |= checkMinAndStorageN("turnover resorption final", 0.0, resorpFlux * climate->length); return status; diff --git a/tests/sipnet/test_modeling/testPlantMortality.c b/tests/sipnet/test_modeling/testPlantMortality.c index 9e5d3b36..d7bace6d 100644 --- a/tests/sipnet/test_modeling/testPlantMortality.c +++ b/tests/sipnet/test_modeling/testPlantMortality.c @@ -44,8 +44,8 @@ void resetEnv(void) { envi.soilOrgN = 0.0; envi.litterN = 0.0; envi.plantStorageN = 0.0; - eventTrackers.harvestFracRemoved = 0.0; - eventTrackers.harvestFracTransferred = 0.0; + eventTrackers.harvestTrackers.totalFracRemoved = 0.0; + eventTrackers.harvestTrackers.totalFracTransferred = 0.0; } int checkPool(double calc, double exp, const char *label) { diff --git a/tests/sipnet/test_sipnet_infrastructure/testDebugLogFiles.c b/tests/sipnet/test_sipnet_infrastructure/testDebugLogFiles.c index 793712b2..d5d918a3 100644 --- a/tests/sipnet/test_sipnet_infrastructure/testDebugLogFiles.c +++ b/tests/sipnet/test_sipnet_infrastructure/testDebugLogFiles.c @@ -3,6 +3,7 @@ #include #include "common/logging.h" +#include "sipnet/events.h" #include "sipnet/state.h" #include "utils/tUtils.h" @@ -15,14 +16,12 @@ #define EXPECTED_ENVI_FIELDS ((int)(sizeof(Envi) / sizeof(double))) #define EXPECTED_FLUX_FIELDS ((int)(sizeof(Fluxes) / sizeof(double))) -// These are not all doubles -// Let's do a little more work for Trackers -#define NUM_TRACKER_INTS 1 -#define TRACKERS_SIZE (int)(sizeof(Trackers) - NUM_TRACKER_INTS * sizeof(int)) -#define EXPECTED_TRACKER_FIELDS \ - ((int)(TRACKERS_SIZE / sizeof(double)) + NUM_TRACKER_INTS) +// The Trackers struct is not all doubles, but the int(s) get padded to 8 bytes +#define EXPECTED_TRACKER_FIELDS ((int)(sizeof(Trackers) / sizeof(double))) #define EXPECTED_PHENOLOGY_FIELDS 3 #define EXPECTED_SURVIVAL_FIELDS 1 +#define EXPECTED_EVENT_TRACKER_FIELDS \ + (int)(sizeof(EventTrackers) / sizeof(double)) static int countTokens(const char *line) { int count = 0; @@ -114,7 +113,7 @@ int run(void) { status |= checkHeader(TRACKERS_FILE, 3 + EXPECTED_TRACKER_FIELDS + EXPECTED_PHENOLOGY_FIELDS + - EXPECTED_SURVIVAL_FIELDS, + EXPECTED_SURVIVAL_FIELDS + EXPECTED_EVENT_TRACKER_FIELDS, "t.gpp", "pt.lastYear"); int mainLines = countLines(SIPNET_OUT_FILE); diff --git a/tests/utils/helpers.c b/tests/utils/helpers.c index ff2e5b91..01fc1959 100644 --- a/tests/utils/helpers.c +++ b/tests/utils/helpers.c @@ -21,6 +21,9 @@ void prepTypesTest(void) { void procEvents(void) { resetFluxes(); - processEvents(); + EventNode *event = getCurrentEvent(); + processEventsForCarbon(event); + processEventsForNitrogen(event); + writeEventsOut(); updatePoolsForEvents(); } From fb63355e3d087321ca206e3671b115d950e57e22 Mon Sep 17 00:00:00 2001 From: Mike Longfritz Date: Fri, 2 Oct 2026 13:47:59 -0400 Subject: [PATCH 14/14] Update for #400 --- docs/CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 5fe7e618..c7e20d6b 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -26,6 +26,8 @@ sections to include in release notes: ### Added ### Fixed +- Terminating harvest events will now remove all biomass at end of time step, instead of allowing small amounts to remain in the plant pools due to growth in that step. (#400) +- Leaf-off mechanics refined to ensure that leaf-off events do not remove more leaf carbon than is available in the leaf pool. (#400) ### Changed - SIPNET will now error instead of warning when an environment pool goes negative. This is a change from previous behavior where SIPNET would log a warning and continue running. (#398)