diff --git a/ST_indivs.c b/ST_indivs.c index d6ca9a42..811c67ed 100644 --- a/ST_indivs.c +++ b/ST_indivs.c @@ -308,12 +308,19 @@ void indiv_proportion_Recovery( IndivType *ndv, int killType,RealF proportionRec // using individual killing year old real size and reduction for making base for calculating proportional recovery RealF prev_reduction = ndv->prv_yr_relsize * proportionKilled; - RealF increase = prev_reduction * proportionRecovery; + RealF increase = prev_reduction * proportionRecovery; ndv->relsize = ndv->relsize + increase; Species_Update_Newsize(ndv->myspecies, increase); if (ZERO(ndv->relsize) || LT(ndv->relsize, 0.0)) { + // this should never happend because `increase` should always be positive + LogError(logfp, LOGWARN, "'indiv_proportion_Recovery': an individual of "\ + "%s reached relsize of 0 and is removed (increase = %.3f): for " \ + "killType = %d, proportionKilled = %.2f, proportionRecovery = %.2f", + Species[ndv->myspecies]->name, increase, + killType, proportionKilled, proportionRecovery); + _delete(ndv); } #undef xF_DELTA @@ -343,8 +350,26 @@ void indiv_Kill_Complete( IndivType *ndv, int killType) { // if(!UseGrid) // insertIndivKill(ndv->id,killType); species_Update_Kills(ndv->myspecies, ndv->age); + +if (Species[ndv->myspecies]->res_grp == 6 && Species[ndv->myspecies]->sp_num == 11) { +printf("'indiv_Kill_Complete' after 'species_Update_Kills': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, ndv->relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, ndv->relsize, Species[11]->est_count); +} + Species_Update_Newsize(ndv->myspecies, -ndv->relsize); - _delete(ndv); + +if (Species[ndv->myspecies]->res_grp == 6 && Species[ndv->myspecies]->sp_num == 11) { +printf("'indiv_Kill_Complete' after 'Species_Update_Newsize': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); +} + + _delete(ndv); // `_delete` updates the Species[ndv->myspecies]->est_count, i.e., removes one individual } @@ -352,7 +377,8 @@ void indiv_Kill_Complete( IndivType *ndv, int killType) { void _delete (IndivType *ndv) { /*======================================================*/ /* PURPOSE */ -/* Local routine to remove the data object of an individual. +/* Local routine to remove the data object of an individual and update + * the number of individuals `Species[ndv->myspecies]->est_count` * Called from indiv_Kill_Complete(). */ /* HISTORY */ @@ -387,9 +413,13 @@ void _delete (IndivType *ndv) { ndv->Next->Prev = ndv->Prev; } - if( --s->est_count == 0) + // update Species[ndv->myspecies]->est_count, i.e., + // remove one individual from tally + if( --s->est_count == 0) { + // if there are no individual left of this species, then remove species + // from resource group and update `est_count` of the resource group rgroup_DropSpecies(sp); - + } if ((s->est_count > 0 && s->IndvHead == NULL) || (s->est_count == 0 && s->IndvHead != NULL)) diff --git a/ST_main.c b/ST_main.c index 09f83a58..8a3e2423 100644 --- a/ST_main.c +++ b/ST_main.c @@ -195,45 +195,121 @@ int main(int argc, char **argv) { /* ------ Begin running the model ------ */ for (year = 1; year <= Globals.runModelYears; year++) { - //printf("------------------------Repetition/year = %d / %d\n", iter, year); + printf("------------------------Repetition/year = %d / %d\n", iter, year); Globals.currYear = year; rgroup_Establish(); + printf("'main' after 'rgroup_Establish': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after 'rgroup_Establish'"); Env_Generate(); rgroup_PartResources(); + printf("'main' after 'rgroup_PartResources': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after 'rgroup_PartResources'"); #ifdef STEPWAT if (!isnull(SXW.debugfile) ) SXW_PrintDebug(0); #endif rgroup_Grow(); + printf("'main' after 'rgroup_Grow': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after 'rgroup_Grow'"); mort_Main(&killedany); + printf("'main' after 'mort_Main': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after 'mort_Main'"); rgroup_IncrAges(); + printf("'main' after 'rgroup_IncrAges': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after 'rgroup_IncrAges'"); // Added functions for Grazing and mort_end_year as proportional killing effect before exporting biomass end of the year grazing_EndOfYear(); - save_annual_species_relsize(); + printf("'main' after 'grazing_EndOfYear': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after 'grazing_EndOfYear'"); + + + save_annual_species_relsize(); + printf("'main' after 'save_annual_species_relsize': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after 'save_annual_species_relsize'"); + mort_EndOfYear(); + printf("'main' after 'mort_EndOfYear': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after 'mort_EndOfYear'"); stat_Collect(year); + printf("'main' after 'stat_Collect': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after 'stat_Collect'"); if (BmassFlags.yearly) output_Bmass_Yearly(year); // Moved kill annual and kill extra growth after we export biomass, and recovery of biomass after fire before the next year _kill_annuals(); - proportion_Recovery(); + printf("'main' after '_kill_annuals': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after '_kill_annuals'"); + + proportion_Recovery(); + printf("'main' after 'proportion_Recovery': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + check_sizes("'main' after 'proportion_Recovery'"); + _kill_extra_growth(); // Check that relsizes match up at end of year after extra growth is removed // may want to wrap this in #ifdef DEBUG once problem is fixed check_sizes("'main' at end of year"); + printf("'main' after '_kill_extra_growth': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); } /* end model run for this year*/ if (MortFlags.summary) { @@ -550,7 +626,7 @@ void check_sizes(const char *chkpt) { if (LT(diff, fabs(spsize - Species[sp]->relsize)) ) { LogError(stdout, LOGWARN, "%s (%d:%d): SP: \"%s\" size error: " - "SP=%.9f, ndv=%.9f\n", + "SP=%.3f, ndv=%.3f", chkpt, Globals.currIter, Globals.currYear, Species[sp]->name, Species[sp]->relsize, spsize); } @@ -558,7 +634,7 @@ void check_sizes(const char *chkpt) { if ( LT(diff, fabs(rgsize -RGroup[rg]->relsize)) ) { LogError(stdout, LOGWARN, "%s (%d:%d): RG \"%s\" size error: " - "RG=%.9f, ndv=%.9f\n", + "RG=%.3f, ndv=%.3f", chkpt, Globals.currIter, Globals.currYear, RGroup[rg]->name, RGroup[rg]->relsize, rgsize); } diff --git a/ST_mortality.c b/ST_mortality.c index aadf05e2..625ecaec 100644 --- a/ST_mortality.c +++ b/ST_mortality.c @@ -148,6 +148,13 @@ void mort_Main( Bool *killed) { if ( GT(g->pr, 1.0) ) { if (++g->yrs_neg_pr >= g->max_stretch) _no_resources( rg); + +if (rg == 6) printf("'mort_Main' after '_no_resources': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + } else { g->yrs_neg_pr = 0; } @@ -157,7 +164,21 @@ void mort_Main( Bool *killed) { /* Take care of mortality types 1 and 2*/ if ( g->use_mort ) { _age_independent( sp ); + +if (rg == 6 && sp == 11) printf("'mort_Main' after '_age_independent': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + _slow_growth( sp ); + +if (rg == 6 && sp == 11) printf("'mort_Main' after '_slow_growth': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + } /* now deal with succulents problems*/ if (g->succulent @@ -176,6 +197,13 @@ void mort_Main( Bool *killed) { case Burrow: _burrow( sp); break; + +if (rg == 6 && sp == 11) printf("'mort_Main' after 'Plot.disturbance': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); + } } /* end for each species*/ @@ -598,14 +626,52 @@ static void _age_independent( const SppIndex sp) { kills[++k] = ndv; } +if (Species[sp]->res_grp == 6 && sp == 11) { + int kk = 0; + printf("'_age_independent': individuals of %s/%s with relsize = %.2f\n", + RGroup[6]->name, Species[11]->name, Species[11]->relsize); + + ForEachIndiv (ndv, Species[sp]) { + printf("\t%d, relsize = %.2f\n", + kk, ndv->relsize); + kk++; + } + + printf("'_age_independent': individuals of %s/%s slated to kill\n", + RGroup[6]->name, Species[11]->name); + + for( n=0; n <= k; n++ ) { + printf("\t%d, relsize = %.2f\n", n, kills[n]->relsize); + } +} + + for( n=0; n <= k; n++ ) { + printf("'_age_independent' calling 'indiv_Kill_Complete': n=%d <= k=%d\n", n, k); indiv_Kill_Complete(kills[n], 9); } +if (Species[sp]->res_grp == 6 && sp == 11) { + printf("'_age_independent' after 'indiv_Kill_Complete': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); +} if (k >= 0) _SomeKillage = TRUE; Mem_Free(kills); + + +if (Species[sp]->res_grp == 6 && sp == 11) { + printf("'_age_independent' after 'Mem_Free': \n" \ + "\t%s, relsize = %.2f, est_count = %d\n" \ + "\t%s, relsize = %.2f, est_count = %d\n", + RGroup[6]->name, RGroup[6]->relsize, RGroup[6]->est_count, + Species[11]->name, Species[11]->relsize, Species[11]->est_count); +} + } /***********************************************************/ diff --git a/ST_resgroups.c b/ST_resgroups.c index 8c552240..df6dff19 100644 --- a/ST_resgroups.c +++ b/ST_resgroups.c @@ -132,13 +132,13 @@ void rgroup_PartResources(void) { static RealF _add_annuals(const GrpIndex rg, const SppIndex sp, const RealF lastyear_relsize) { /*======================================================*/ - /* Establishment for annual species, which includes a function that modifies - * the seedbank through seed addition and deletion, a function that determines - * the amount of viable seed in the seedbank if regen_ok is true, a random draw - * from the beta distribution that is used to determine the number of seeds - * that will emerge as seedlings this year, and removal of those germinated + /* Establishment for annual species, which includes a function that modifies + * the seedbank through seed addition and deletion, a function that determines + * the amount of viable seed in the seedbank if regen_ok is true, a random draw + * from the beta distribution that is used to determine the number of seeds + * that will emerge as seedlings this year, and removal of those germinated * seeds from the seedbank.*/ - + IntU i, num_est; //number of individuals that will establish this year RealF viable_seeds; //number of viable seeds in the seedbank float var; //random number draw from beta distribution for calculation of num_est @@ -148,30 +148,30 @@ static RealF _add_annuals(const GrpIndex rg, const SppIndex sp, const RealF last g = RGroup[rg]; s = Species[sp]; - /*Increment viable years for seeds and implement the decay process (seed mortality). + /*Increment viable years for seeds and implement the decay process (seed mortality). * Then add seeds to the seedbank*/ _add_annual_seedprod(sp, lastyear_relsize); /*Get viable seeds from seed bank*/ viable_seeds = (g->regen_ok) ? _get_annual_maxestab(sp) : 0; - /* Create a beta random number draw based on alpha and beta for each species + /* Create a beta random number draw based on alpha and beta for each species * (calculated based on mean (s->seedling_estab_prob and variance (s->var)) */ var = RandBeta(Species[sp]->alpha, Species[sp]->beta); - - /*Determine number of seedlings to add. If the number of seeds calculated + + /*Determine number of seedlings to add. If the number of seeds calculated * from the random draw is larger than max_seed_estab, max_seed_estab is used instead*/ - num_est = min(viable_seeds * var, s->max_seed_estab); + num_est = min(viable_seeds * var, s->max_seed_estab); //printf("Species name=%s , num_est =%u \n",s->name, num_est); - - /*Multiple the proportion of seeds in each viable year array by the total - number of seeds that germinated as seedlings and subtract those seeds from + + /*Multiple the proportion of seeds in each viable year array by the total + number of seeds that germinated as seedlings and subtract those seeds from the relevant seedprod array.*/ for (i = 0; i <= s->viable_yrs; i++) { - //printf("Species name=%s , old calculated value s->seedprod[%hu]= %d \n", s->name, i, s->seedprod[i]); - s->seedprod[i] = s->seedprod[i] - round(num_est * s->seedprod[i] / viable_seeds); + //printf("Species name=%s , old calculated value s->seedprod[%hu]= %d \n", s->name, i, s->seedprod[i]); + s->seedprod[i] = s->seedprod[i] - round(num_est * s->seedprod[i] / viable_seeds); //printf("Species name=%s , so new calculated value s->seedprod[%hu]= %d \n", s->name, i, s->seedprod[i]); - } + } return num_est; } @@ -198,7 +198,7 @@ static void _add_annual_seedprod(SppIndex sp, RealF lastyear_relsize) { SpeciesType *s = Species[sp]; IntU i; - /*Age of seeds is increased by one year, seed mortality occurs, and seeds produced in the + /*Age of seeds is increased by one year, seed mortality occurs, and seeds produced in the previous year are added to the seedbank at relative age 0 */ for (i = s->viable_yrs - 1; i > 0; i--) { //printf("Species name=%s , seedprod before decay s->seedprod[%hu]= %d \n", s->name, i, s->seedprod[i]); @@ -208,10 +208,10 @@ static void _add_annual_seedprod(SppIndex sp, RealF lastyear_relsize) { // printf("Species name=%s ,old Array array 0 index i=%u, value =%hu \n", s->name, i, s->seedprod[i]); - /* If the current year is year 1 of the simulation, then the number of seeds - * added is a random number draw between 1 and the maximum number of seedlings - * that can establish. Otherwise, this year's seed production is a function - * of the number of seeds produced per unit biomass multiplied by species biomass + /* If the current year is year 1 of the simulation, then the number of seeds + * added is a random number draw between 1 and the maximum number of seedlings + * that can establish. Otherwise, this year's seed production is a function + * of the number of seeds produced per unit biomass multiplied by species biomass * (maximum species biomass * last year's species relative size). */ if (Globals.currYear == 1) { s->seedprod[0] = RandUniRange(1, s->max_seed_estab); @@ -576,33 +576,44 @@ static void _extra_growth(GrpIndex rg) { } /*END ForEachIndiv */ - //printf("s->relsize before = %f\n, Species = %s \n", Species[sp]->name, Species[sp]->relsize); - //printf("s->extragrowth out of loop = %f\n, Species = %s \n", Species[sp]->name,s->extragrowth); + /*printf("%s(%d): relsize before = %.3f, extragrowth = %.3f " \ + "relsize + extra = %.3f || %s(%d): relsize before = %.3f\n", + Species[sp]->name, sp, Species[sp]->relsize, s->extragrowth, + Species[sp]->relsize + s->extragrowth, + RGroup[rg]->name, rg, RGroup[rg]->relsize); + */ Species_Update_Newsize(sp, Species[sp]->extragrowth); - //printf("s->relsize after = %f\n, Species = %s \n", Species[sp]->name, Species[sp]->relsize); + + /*printf("%s: relsize after = %.3f || %s: relsize after = %.3f\n", + Species[sp]->name, Species[sp]->relsize, + RGroup[rg]->name, RGroup[rg]->relsize); + */ + } /* ENDFOR j (for each species)*/ } + + /***********************************************************/ void rgroup_Establish(void) { /*======================================================*/ /* PURPOSE */ - /* Determines which and how many species can establish in a given year. For + /* Determines which and how many species can establish in a given year. For * each species in each perennial functional group, check that a uniform - * random number between 0 and 1 is less than the species' establishment - * probability. a) If so, return a random number of individuals up to the - * maximum allowed to establish for the species. This is the number of individuals - * in this species that will establish this year. b) If not, continue with - * the next species. Establishment for species of annual functional groups + * random number between 0 and 1 is less than the species' establishment + * probability. a) If so, return a random number of individuals up to the + * maximum allowed to establish for the species. This is the number of individuals + * in this species that will establish this year. b) If not, continue with + * the next species. Establishment for species of annual functional groups * occurs differently. See notes at the top of _add_annuals */ /* HISTORY */ /* Chris Bennett @ LTER-CSU 6/15/2000 */ - /* The probability of establishment emulates the occurrence of microsite - * conditions that allow for establishment. - * 7-Nov-03 (cwb) Adding the new algorithm to handle annuals. It's more - * complicated than before (which didn't really work) so annuals are now added - * in the PartResources()function. Only perennials are added here. Also, there's + /* The probability of establishment emulates the occurrence of microsite + * conditions that allow for establishment. + * 7-Nov-03 (cwb) Adding the new algorithm to handle annuals. It's more + * complicated than before (which didn't really work) so annuals are now added + * in the PartResources()function. Only perennials are added here. Also, there's * now a parameter to define the start year of establishment for perennials. * KAP: Annual establishment now occurs here instead of in PartResources*/ /*------------------------------------------------------*/