diff --git a/ST_indivs.c b/ST_indivs.c index 95d6ecf4..4267dd07 100644 --- a/ST_indivs.c +++ b/ST_indivs.c @@ -225,7 +225,6 @@ void indiv_proportion_Kill(IndivType *ndv, int killType, RealF proportKilled) // insertIndivKill(ndv->id, killType); //kill indiv Proportionally or adjust their real size irrespective of being annual or perennial, both will have this effect - species_Update_Kills(ndv->myspecies, ndv->age); // saving killing year real size here that is going to use for calculating next year proportional recovery ndv->prv_yr_relsize = ndv->relsize; @@ -238,6 +237,8 @@ void indiv_proportion_Kill(IndivType *ndv, int killType, RealF proportKilled) if (ZERO(ndv->relsize) || LT(ndv->relsize, 0.0)) { ndv->relsize =0.0; + // increase mortality count only if relsize has become zero due to fire + species_Update_Kills(ndv->myspecies, ndv->age); } #undef xF_DELTA diff --git a/ST_mortality.c b/ST_mortality.c index 9a98dc77..5b22f092 100644 --- a/ST_mortality.c +++ b/ST_mortality.c @@ -146,6 +146,9 @@ void mort_Main( Bool *killed) { /* annuals are not subject to these sources of mortality and instead die in _kill_annuals */ if (g->max_age == 1) continue; + /* Calculate PR at the functional group level: resources required/resources available */ + g->pr = ZRO(g->res_avail) ? 0. : g->res_required / g->res_avail; + /* kill plants if low resources for consecutive years */ /* increment yrs_neg_pr if pr > 1, else zero it. */ /* one good year cancels all previous bad years. */ diff --git a/ST_params.c b/ST_params.c index 53f5d1ba..dc03fdb5 100644 --- a/ST_params.c +++ b/ST_params.c @@ -909,10 +909,12 @@ static void _rgroup_add1( char name[], RealF space, RealF density, RGroup[rg]->grp_num = rg; RGroup[rg]->max_stretch = (IntS) stretch; RGroup[rg]->max_spp_estab = (IntS) estab; - RGroup[rg]->max_density = density; - RGroup[rg]->max_per_sqm = density / Globals.plotsize; + // input of `density` is in units of [# / m2]; convert to units of [# / plot] + RGroup[rg]->max_density = density * Globals.plotsize; // density per plot + RGroup[rg]->max_per_sqm = density; // density per square-meter RGroup[rg]->use_mort = itob(mort); RGroup[rg]->slowrate = slow; + RGroup[rg]->space = space; RGroup[rg]->min_res_req = space; RGroup[rg]->est_annually = itob(estann); RGroup[rg]->startyr = styr; @@ -1107,7 +1109,8 @@ static void _species_init( void) { Species[sp]->received_prob = 0; Species[sp]->cohort_surv = cohort; Species[sp]->var = var; - Species[sp]->pseed = pseed / Globals.plotsize; + // input of `pseed` is in units of [# / m2]; convert to units of [# / plot] + Species[sp]->pseed = pseed * Globals.plotsize; /* Species[sp]->ann_mort_prob = (age > 0) ? -log(cohort)/age : 0.0; diff --git a/ST_resgroups.c b/ST_resgroups.c index 1749cd49..c16a508f 100644 --- a/ST_resgroups.c +++ b/ST_resgroups.c @@ -100,10 +100,6 @@ void rgroup_PartResources(void) { LogError(logfp, LOGWARN, "RGroup %s : res_avail is Zero and res_required > 0", g->name); } - /* Calculate PR at the functional group level: resources required/resources available */ - g->pr = ZRO(g->res_avail) ? 0. : g->res_required / g->res_avail; - //printf("g->pr = %f\n,Group = %s \n",RGroup[rg]->name, g->pr); - /* If relsize>0 and individuals are established, reset noplants from TRUE to FALSE */ if (GT(getRGroupRelsize(rg), 0.)) noplants = FALSE; @@ -257,7 +253,7 @@ static void _res_part_extra(RealF extra, RealF size[]) { continue; /* Check to avoid dividing by 0 */ - if (sum_size == 0.) + if (ZRO(sum_size)) req_prop = 0.; /* Calculate proportional biomass of each group out of the total biomass @@ -267,9 +263,10 @@ static void _res_part_extra(RealF extra, RealF size[]) { /* If the group can use extra resources, divide out extra based on * proportional biomass */ - if (g->use_extra_res) + if (g->use_extra_res) { g->res_extra = req_prop * extra; - + g->res_avail += g->res_extra; + } /* If the group can't use extra resources, set res_extra to 0 */ else g->res_extra = 0.; @@ -384,9 +381,8 @@ void rgroup_ResPartIndiv(void) { /* If individuals already have the resources they require do * not assign extra */ - if (ndv->res_avail == ndv->res_required) { + if (GE(ndv->res_avail, ndv->res_required)) { ndv->res_extra = 0.0; - //printf("ndv->res_extra = %f\n", ndv->res_extra); } /* Assign extra resource as the difference of what is required @@ -471,8 +467,7 @@ void rgroup_Grow(void) { continue; /* Modify growth rate by temperature calculated in Env_Generate() */ - if (s->tempclass != NoSeason) - tgmod = Env.temp_reduction[s->tempclass]; + tgmod = (s->tempclass == NoSeason) ? 1. : Env.temp_reduction[s->tempclass]; /* Now increase size of the individual plants of current species */ ForEachIndiv(ndv, s) { @@ -616,6 +611,8 @@ void rgroup_Establish(void) { /*------------------------------------------------------*/ IntS i, num_est; /* number of individuals of sp. that establish*/ + RealF + used_space = 0; /* sums up all of the space that is currently used */ GrpIndex rg; SppIndex sp; GroupType *g; @@ -634,6 +631,7 @@ void rgroup_Establish(void) { continue; g->regen_ok = TRUE; /* default */ + g->min_res_req = g->space; /* reset min_res_req, if it was modified last year */ if (Globals.currYear < RGroup[rg]->startyr) { g->regen_ok = FALSE; @@ -652,7 +650,6 @@ void rgroup_Establish(void) { if (Species[sp]->max_age == 1) { //printf("Globals.currYear = %hu, call to _add_annuals sp=%d Species[sp]->lastyear_relsize : %.5f \n", Globals.currYear, sp, Species[sp]->lastyear_relsize); num_est = _add_annuals(rg, sp, Species[sp]->lastyear_relsize); - // printf("g->seedbank annuals=%d \n",g->seedbank); } /* Establishment for species that belong to perennial functional groups*/ @@ -670,8 +667,30 @@ void rgroup_Establish(void) { } /* end ForEachGroupSpp() */ } + + // sum up min_res_req in case we need to redistribute min_res_req + if (g->est_count > 0) { + used_space += g->min_res_req; + } else { + // this group needs no resources because it is not established. + g->min_res_req = 0; + } } /* end ForEachGroup() */ + + // If there is unused (or too much used) space we need to redistribute + if (!EQ(used_space, 1.0)) { + + ForEachGroup(rg) { + g = RGroup[rg]; + /* Redistribute the unused (or too much used) space to all groups that + are established */ + if (g->est_count > 0) { + g->min_res_req = g->min_res_req / used_space; + } + } + } } + /***********************************************************/ void rgroup_IncrAges(void) { @@ -717,17 +736,18 @@ void rgroup_IncrAges(void) /* Sums relsize for all individuals in all species in RGroup rg. param rg = RGroup index. Return: RGroup relsize. */ + RealF getRGroupRelsize(GrpIndex rg){ Int n; SppIndex sp; double sum = 0.0; ForEachEstSpp( sp, rg, n){ - sum += getSpeciesRelsize(sp); + sum += getSpeciesRelsize(sp); } if(RGroup[rg]->est_count > 0){ - return (RealF) sum / (RealF) RGroup[rg]->est_count; + return (RealF) sum; } else { return 0; } diff --git a/ST_structs.h b/ST_structs.h index ad38c18f..c1027e08 100644 --- a/ST_structs.h +++ b/ST_structs.h @@ -90,7 +90,7 @@ struct species_st { estabs, /* number of individuals established in iter */ *seedprod, /* annuals: array of previous years' seed production (size = viable_yrs)*/ seedbank, - pseed; + pseed; /* average number of seeds produced by annual species per 1g of biomass, per 1m^2 and per year (internally re-calculated as seeds per 1 g biomass per plot and per year) */ RealF lastyear_relsize, /* relsize from the previous year, used for annual establishment */ extragrowth, /* amt of superfluous growth from extra resources */ received_prob, //the chance that this species received seeds this year... only applicable if using seed dispersal and gridded option @@ -180,9 +180,10 @@ struct resourcegroup_st { grazingfrq, /* grazing effect on group at this frequency: <1=prob, >1=# years */ grazingfreq_startyr;/* start year for grazing frequency*/ SppIndex *species; /*list of spp belonging to this grp*/ - RealF min_res_req, /* input from table */ + RealF space, /* input from table */ + min_res_req, /* input space from table, rescaled if one or more rgroups is not established */ max_density, /* number of mature plants per plot allowed */ - max_per_sqm, /* convert density and plotsize to max plants/m^2 */ + max_per_sqm, /* density of mature plants in units of plants / m^2 */ max_bmass, /* sum of mature biomass for all species in group */ killfreq, /* kill group at this frequency: <1=prob, >1=# years */ ignition, /* cheatgrass biomass (g/m2) that triggers potential ignition of a wildfire */ diff --git a/sxw_resource.c b/sxw_resource.c index 4770544e..94a35b8f 100644 --- a/sxw_resource.c +++ b/sxw_resource.c @@ -79,7 +79,7 @@ extern extern transp_t transp_window; -extern +extern pcg32_random_t resource_rng; // _transp_contribution_by_group() has different behavior for production years and setup years. @@ -128,13 +128,13 @@ void _sxw_update_resource(void) { * The first step is to get the current biomass for each STEPPE functional group. * Second, we re-calculate the relative active roots in each layer in each month * using the updated functional group biomass. Third, we divide transpiration to - * each STEPPE functional group based on the matching of active roots in each - * soil layer and month with transpiration in each layer and month. Finally, we + * each STEPPE functional group based on the matching of active roots in each + * soil layer and month with transpiration in each layer and month. Finally, we * scale resources available (cm) to resources in terms of grams of biomass */ - + RealF *sizes; GrpIndex g; - + sizes = (RealF *)Mem_Calloc(Globals.max_rgroups, sizeof(RealF), "_sxw_update_resource"); ForEachGroup(g) @@ -145,14 +145,14 @@ void _sxw_update_resource(void) { continue; sizes[g] = RGroup_GetBiomass(g); } - + /* Update the active relative roots based on current biomass values */ _sxw_update_root_tables(sizes); - + /* Assign transpiration (resource availability) to each STEPPE functional group */ _transp_contribution_by_group(_resource_cur); - - /* Scale transpiration resources by a constant, bvt, to convert resources + + /* Scale transpiration resources by a constant, bvt, to convert resources * (cm) to biomass that can be supported by those resources (g/cm) */ ForEachGroup(g) { @@ -161,7 +161,7 @@ void _sxw_update_resource(void) { //printf("for groupName= %s, resource_cur post multiplication: %f\n\n",Rgroup[g]->name, _resource_cur[g]); } /* _print_debuginfo(); */ - + Mem_Free(sizes); } @@ -174,7 +174,8 @@ void _sxw_update_root_tables( RealF sizes[] ) { GrpIndex g; LyrIndex l, nLyrs; TimeInt p; - RealD x; + RealD x, sum_all_veg_types; + int t; /* Set some things to zero where 4 refers to Tree, Shrub, Grass, Forb */ Mem_Set(_roots_active_sum, 0, NVEGTYPES * SXW.NPds * SXW.NTrLyrs * sizeof(RealD)); @@ -194,18 +195,26 @@ void _sxw_update_root_tables( RealF sizes[] ) { } } - /* Rescale _roots_active_sum to represent the relative "activity" of a + /* Rescale _roots_active_sum to represent the relative "activity" of a * STEPPE group's roots in a given layer in a given month */ + /* Details: for each soil layer `l` and each month (trperiod) `p`, the sum + of `_roots_active_rel[Iglp(g, l, p)]` across rgroups `g` must be 1; + see its use by function _transp_contribution_by_group() */ ForEachGroup(g) { nLyrs = getNTranspLayers(RGroup[g]->veg_prod_type); for (l = 0; l < nLyrs; l++) { ForEachTrPeriod(p) { + sum_all_veg_types = 0; + ForEachVegType(t) { + sum_all_veg_types += _roots_active_sum[Itlp(t, l, p)]; + } + _roots_active_rel[Iglp(g, l, p)] = - ZRO(_roots_active_sum[Itlp(RGroup[g]->veg_prod_type, l, p)]) ? + ZRO(sum_all_veg_types) ? 0. : - _roots_active[Iglp(g, l, p)] / _roots_active_sum[Itlp(RGroup[g]->veg_prod_type, l, p)]; + _roots_active[Iglp(g, l, p)] / sum_all_veg_types; } } } @@ -214,7 +223,7 @@ void _sxw_update_root_tables( RealF sizes[] ) { static void _transp_contribution_by_group(RealF use_by_group[]) { /*======================================================*/ - /* use_by_group is the amount of transpiration (cm) assigned to each STEPPE + /* use_by_group is the amount of transpiration (cm) assigned to each STEPPE * functional group. Must call _update_root_tables() before this. * Compute each group's amount of transpiration from SOILWAT2 * based on its biomass, root distribution, and phenological @@ -224,13 +233,11 @@ static void _transp_contribution_by_group(RealF use_by_group[]) { GrpIndex g; TimeInt p; + int nLyrs, t; LyrIndex l; - int t; - RealD *transp; + RealD *transp, proportion_total_resources, average_proportion; RealF sumUsedByGroup = 0., sumTranspTotal = 0., TranspRemaining = 0.; - RealF transp_ratio; - added_transp = 0; - RealF transp_ratio_sd; + RealF transp_ratio, added_transp = 0, transp_ratio_sd; // year 0 is a set up year. No need to calculate transpiration. // if there are multiple iterations the last year will run twice; @@ -252,13 +259,16 @@ static void _transp_contribution_by_group(RealF use_by_group[]) { ForEachGroup(g) //Steppe functional group { use_by_group[g] = 0.; - transp = SXW.transpVeg[RGroup[g]->veg_prod_type]; + transp = SXW.transpTotal; //Loops through each month and calculates amount of transpiration for each STEPPE functional group //according to whether that group has active living roots in each soil layer for each month - + /* Details: for each soil layer `l` and each month (trperiod) `p`, the + sum of `_roots_active_rel[Iglp(g, l, p)]` across rgroups `g` must + be 1; only if they sum to 1, can they can be used as proper weights + for `transp[Ilp(l, p)]` to be split up among rgroups */ ForEachTrPeriod(p) { - int nLyrs = getNTranspLayers(RGroup[g]->veg_prod_type); + nLyrs = getNTranspLayers(RGroup[g]->veg_prod_type); for (l = 0; l < nLyrs; l++) { use_by_group[g] += (RealF) (_roots_active_rel[Iglp(g, l, p)] * transp[Ilp(l, p)]); } @@ -269,16 +279,32 @@ static void _transp_contribution_by_group(RealF use_by_group[]) { //printf(" sumUsedByGroup in transp=%f \n",sumUsedByGroup); } + // Rescale transpiration based on space (min_res_req) + ForEachGroup(g){ + if(use_by_group[g] != 0){ //Prevents NaN errors + //Proportion of total traspiration that this group would recieve without accounting for space + proportion_total_resources = use_by_group[g] / sumUsedByGroup; + //Average of the two proportions. This equally weights proportion_total_resources and min_res_req. + average_proportion = (proportion_total_resources + RGroup[g]->min_res_req) / 2; + //Recalculate this group's resources. + use_by_group[g] = sumUsedByGroup * average_proportion; + } + } + + //sumUsedByGroup needs to be synced because of potention floating point errors + sumUsedByGroup = 0; + ForEachGroup(g){ + sumUsedByGroup += use_by_group[g]; + } + //Very small amounts of transpiration remain and not perfectly partitioned to functional groups. //This check makes sure any remaining transpiration is divided proportionately among groups. - ForEachTrPeriod(p) { for (t = 0; t < SXW.NSoLyrs; t++) sumTranspTotal += SXW.transpTotal[Ilp(t, p)]; } TranspRemaining = sumTranspTotal - sumUsedByGroup; - //printf(" sumTranspTotal=%f, sumUsedByGroup=%f TranspRemaining=%f"\n", sumTranspTotal, sumUsedByGroup, TranspRemaining); /* ------------- Begin testing to see if additional transpiration is necessary ------------- */ @@ -336,7 +362,7 @@ static void _transp_contribution_by_group(RealF use_by_group[]) { RealF min = transp_window.ratio_average - transp_ratio_sd; RealF max = transp_window.ratio_average + transp_ratio_sd; - // This transpiration will be added + // This transpiration will be added added_transp = (1 - transp_ratio / RandUniFloatRange(min, max, &resource_rng)) * transp_window.average; if(added_transp < 0){ LogError(logfp, LOGNOTE, "sxw_resource: Added transpiration less than 0.\n"); @@ -362,10 +388,7 @@ static void _transp_contribution_by_group(RealF use_by_group[]) { // replace the sum of squares with what we just calculated transp_window.SoS_array[transp_window.oldest_index] = ssqr; - //printf("Year %d: ratio = %f, mean = %f, sos = %f\n",Globals.currYear, - //transp_ratio+added_transp_ratio,transp_window.ratio_average, transp_window.sum_of_sqrs); - - /* Adds the additional transpiration to the remaining transpiration + /* Adds the additional transpiration to the remaining transpiration * so it can be distributed proportionally to the functional groups. */ TranspRemaining += added_transp; } @@ -384,7 +407,7 @@ static void _transp_contribution_by_group(RealF use_by_group[]) { //printf("for groupName= %s, after sum use_by_group[g]= %f \n",RGroup[g]->name,use_by_group[g]); } } - + //remember the last year. When setting up for a new iteration the same year will appear twice, and we want to skip it the second time - lastYear = Globals.currYear; + lastYear = Globals.currYear; } diff --git a/sxw_soilwat.c b/sxw_soilwat.c index 7583dc59..1981d979 100644 --- a/sxw_soilwat.c +++ b/sxw_soilwat.c @@ -93,8 +93,6 @@ void _sxw_sw_setup (RealF sizes[]) { v->veg[k].lai_conv_daily[doy] = 0.; } } - - SW_VPD_init(); } void _sxw_sw_run(void) { diff --git a/testing.sagebrush.master/Stepwat_Inputs/Input/rgroup.in b/testing.sagebrush.master/Stepwat_Inputs/Input/rgroup.in index 9528aaac..866ff65f 100644 --- a/testing.sagebrush.master/Stepwat_Inputs/Input/rgroup.in +++ b/testing.sagebrush.master/Stepwat_Inputs/Input/rgroup.in @@ -14,9 +14,9 @@ # name = a name to give the group (12 chars) # space = proportion of the resource space used by the group. # this must sum to 1.0 -# density = relative number of mature individuals per plot, reassigned to max_density in the model code -# maxest = maximum species that can establish in a given year, -# this value cannot be greater than the number of species assigned +# density = number of mature individuals per square-meter, converted to `max_density` (in units of plants per plot) in the model code +# maxest = maximum species that can establish in a given year, +# this value cannot be greater than the number of species assigned # to the resource group in the species.in file # slow = slow growth rate, defines when to count stretched years # stretch = max yrs resources can be stretched before low-resource mortality occurs. diff --git a/testing.sagebrush.master/Stepwat_Inputs/Input/species.in b/testing.sagebrush.master/Stepwat_Inputs/Input/species.in index e18907b5..68e94f7e 100644 --- a/testing.sagebrush.master/Stepwat_Inputs/Input/species.in +++ b/testing.sagebrush.master/Stepwat_Inputs/Input/species.in @@ -50,7 +50,7 @@ # along with the mean (pestab), is used to calculate the two shape parameters: alpha and beta. # Default value: 0.0001. NOTE: var must be < pestab * (1-pestab) or alpha and/or beta will be negative, # we will no longer be meeting the assumptions of a beta distribution, and the code will fail. -# pseed = the average number of seeds produced by annual species per 1g of biomass, per 1m^2 and per year. +# pseed = the average number of seeds produced by annual species per 1g of biomass, per 1m^2 and per year (internally re-calculated as seeds per 1 g biomass per plot and per year) ###PARAMETER SOURCES #TM - eind values are density values taken from Adler datasets from Idaho. Species vuoc (1090), bogr, spcr, brar (2080), are from Karl et al 1999 @@ -78,21 +78,21 @@ # ALSO tested max age for ARTR - found that model is not very sensitive to max age but literature shows 100-150 max age # name rg age irate ratep slow disturb pestab eind minbio maxbio clonal vegindv tclass cosurv onoff var pseed - artr 1 100 0.10000 0.900 4 3 0.3000 2 2.000 400.00 n 0 2 0.0100 1 0.0 0 + artr 1 100 0.35000 0.900 4 3 0.1500 1 2.000 1000.00 n 0 2 0.0100 1 0.0 0 trdu 2 1 0.94700 0.900 2 3 0.0300 5 0.035 4.700 n 0 2 0.0000 0 0.0 0 cryp 2 1 0.94700 0.900 2 3 0.2000 10 0.035 4.700 n 0 2 0.0000 1 0.0001 100 amre 3 1 0.94700 0.900 2 3 0.0300 5 0.035 4.429 n 0 2 0.0000 0 0.0 0 chen 3 1 0.94700 0.900 2 2 0.2000 8 0.035 4.429 n 0 1 0.0000 1 0.0001 100 acmi 4 10 0.42600 0.900 2 3 0.1500 6 0.035 5.000 n 0 2 0.0100 0 0.0 0 - phho 4 10 0.42600 0.900 2 3 0.1500 15 0.035 5.000 n 0 2 0.0100 1 0.0 0 + phho 4 10 0.42600 0.900 2 3 0.1500 15 0.035 22.000 n 0 2 0.0100 1 0.0 0 arfr 5 5 0.42600 0.900 2 3 0.1000 5 0.035 7.480 n 0 2 0.0100 1 0.0 0 lipu 5 35 0.42600 0.900 2 3 0.1000 4 0.035 7.480 n 0 2 0.0100 0 0.0 0 brte 6 1 0.94700 0.900 2 3 0.1500 20 0.020 5.000 n 0 2 0.0000 1 0.0001 1000 vuoc 6 1 0.94700 0.900 2 1 0.0300 16 0.020 5.000 n 0 2 0.0000 0 0.0 0 pssp 7 35 0.47400 0.900 2 3 0.1500 20 0.500 15.096 y 3 2 0.0100 1 0.0 0 koma 7 5 0.47400 0.900 2 3 0.1500 2 0.500 15.096 y 3 2 0.0100 0 0.0 0 - bogr 8 35 0.47400 0.900 2 2 0.0500 2 0.605 15.096 y 3 1 0.0100 0 0.0 0 - spcr 8 29 0.47400 0.900 2 3 0.0500 3 0.605 6.000 y 0 1 0.0100 1 0.0 0 + bogr 8 35 0.47400 0.900 2 2 0.0500 20 0.605 15.096 y 3 1 0.0100 1 0.0 0 + spcr 8 29 0.47400 0.900 2 3 0.0500 3 0.605 6.000 y 0 1 0.0100 0 0.0 0 gusa 9 35 0.73700 0.900 2 3 0.1500 2 0.986 12.126 n 0 2 0.0100 0 0.0 0 chvi 9 17 0.47400 0.900 2 2 0.1500 2 0.986 70.726 n 0 2 0.0100 1 0.0 0 oppo 10 30 0.28900 0.900 4 4 0.0500 1 2.250 15.00 y 3 1 0.0100 1 0.0 0