Merge pull request #550 from DrylandEcology/forage_output

Forage output

ST_functions.h

@@ -64,7 +64,7 @@ SppIndex Species_Name2Index (const char *name);
 void Species_Kill (const SppIndex sp, int killType);
 void Species_Proportion_Kill (const SppIndex sp, int killType, RealF proportionKilled );
 void Species_Proportion_Recovery (const SppIndex sp, int killType, RealF proportionRecovery,RealF proportionKilled);
-void Species_Proportion_Grazing(const SppIndex sp,RealF proportionGrazing );
+RealF Species_Proportion_Grazing(const SppIndex sp,RealF proportionGrazing );
 void Species_Annual_Kill(const SppIndex sp, int killType);
 IntS Species_NumEstablish( SppIndex sp);
 

ST_grid.c

@@ -524,6 +524,15 @@ static void _allocate_accumulators(void){
 			   gridCells[i][j].<accumulator> because the ST_stats accumulators are local. */
 			load_cell(i,j);
 
+			if (BmassFlags.graz) {
+				gridCells[i][j]._Grazed = (StatType*)Mem_Calloc(SuperGlobals.max_rgroups, sizeof(StatType),
+					"allocate_accumulators(Grazed)");
+				ForEachGroup(rg) {
+					gridCells[i][j]._Grazed[rg].s = (struct accumulators_st*)Mem_Calloc(SuperGlobals.runModelYears,
+						sizeof(struct accumulators_st), "_allocate_accumulators(Grazed)");
+				}
+			}
+
   			if (BmassFlags.dist) {
     			gridCells[i][j]._Dist = (StatType*) Mem_Calloc(1, sizeof(StatType), "_allocate_accumulators(Dist)");
 		    	gridCells[i][j]._Dist->s = (struct accumulators_st *)
@@ -798,7 +807,7 @@ void load_cell(int row, int col){
 	                       gridCells[row][col]._Grp, gridCells[row][col]._Gsize, gridCells[row][col]._Gpr,
 						   gridCells[row][col]._Gmort, gridCells[row][col]._Gestab, gridCells[row][col]._Spp,
 						   gridCells[row][col]._Indv, gridCells[row][col]._Smort, gridCells[row][col]._Sestab,
-						   gridCells[row][col]._Sreceived, gridCells[row][col]._Gwf, gridCells[row][col].stats_init);
+						   gridCells[row][col]._Sreceived, gridCells[row][col]._Grazed, gridCells[row][col]._Gwf, gridCells[row][col].stats_init);
 
 	/* Copy this cell's SXW variables into the local variables in sxw.c */
 	copy_sxw_variables(gridCells[row][col].mySXW, gridCells[row][col].mySXWResources, gridCells[row][col].myTranspWindow);
@@ -851,7 +860,7 @@ void unload_cell(){
 	Plot = NULL;
 	Globals = NULL;
 	// Nullify the accumulators
-	stat_Copy_Accumulators(NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,FALSE);
+	stat_Copy_Accumulators(NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,FALSE);
 	// Nullify sxw
 	copy_sxw_variables(NULL,NULL,NULL);
 }
@@ -1338,8 +1347,8 @@ void _Output_AllCellAvgBmass(const char * filename){
 	/* One accumulator for every accumulator in ST_stats.c */
 	float ppt, pptstd, pptsos, temp, tempstd, tempsos, dist, wildfire, grp[SuperGlobals.max_rgroups], grpstd[SuperGlobals.max_rgroups], 
 		  grpsos[SuperGlobals.max_rgroups], gsize[SuperGlobals.max_rgroups], gpr[SuperGlobals.max_rgroups], 
-		  gprsos[SuperGlobals.max_rgroups], gprstd[SuperGlobals.max_rgroups], prescribedfire[SuperGlobals.max_rgroups],
-		  spp[SuperGlobals.max_spp_per_grp * SuperGlobals.max_rgroups],
+		  gprsos[SuperGlobals.max_rgroups], gprstd[SuperGlobals.max_rgroups], graze[SuperGlobals.max_rgroups],
+		  prescribedfire[SuperGlobals.max_rgroups], spp[SuperGlobals.max_spp_per_grp * SuperGlobals.max_rgroups],
 		  indv[SuperGlobals.max_spp_per_grp * SuperGlobals.max_rgroups];
 
 	char buf[2048], tbuf[2048];	// Two buffers: one for accumulating and one for formatting.
@@ -1380,6 +1389,7 @@ void _Output_AllCellAvgBmass(const char * filename){
 			gprsos[rg] = 0;
 			gprstd[rg] = 0; 
 			prescribedfire[rg] = 0;
+			graze[rg] = 0;
 		}
 		ForEachSpecies(sp){
 			spp[sp] = 0;
@@ -1391,6 +1401,7 @@ void _Output_AllCellAvgBmass(const char * filename){
 			for(j = 0; j < grid_Cols; ++j){ // For each column
 				/* ------------- Accumulate requested output ----------------- */
 				nobs++;
+
 				if(BmassFlags.ppt) {
 
 					float old_ppt_ave = ppt;
@@ -1406,6 +1417,11 @@ void _Output_AllCellAvgBmass(const char * filename){
 				}
 				if(BmassFlags.dist) dist += gridCells[i][j]._Dist->s[year].nobs;
 				if(BmassFlags.grpb) {
+					if (BmassFlags.graz) {
+						ForEachGroup(rg) {
+							graze[rg] = get_running_mean(nobs, graze[rg], gridCells[i][j]._Grazed[rg].s[year].ave);
+						}
+					}
 					if(BmassFlags.wildfire){
 						wildfire += gridCells[i][j]._Gwf->wildfire[year];
 					}
@@ -1495,6 +1511,10 @@ void _Output_AllCellAvgBmass(const char * filename){
 					sprintf(tbuf, "%f%c", prescribedfire[rg], sep);
 					strcat(buf, tbuf);
 				}
+				if (BmassFlags.graz) {
+					sprintf(tbuf, "%f%c", graze[rg], sep);
+					strcat(buf, tbuf);
+				}
 			}
 		}
 		if(BmassFlags.sppb){

ST_grid.h

@@ -93,7 +93,7 @@ typedef struct grid_cell_st
 	/* ---------------- accumulators -------------------- */
 	StatType *_Dist, *_Ppt, *_Temp,
   		*_Grp, *_Gsize, *_Gpr, *_Gmort, *_Gestab,
-  		*_Spp, *_Indv, *_Smort, *_Sestab, *_Sreceived;
+  		*_Spp, *_Indv, *_Smort, *_Sestab, *_Sreceived, *_Grazed;
 	FireStatsType *_Gwf;
 	Bool stats_init;
 	/* -------------- end accumulators ------------------ */

ST_indivs.c

@@ -44,7 +44,7 @@ Bool indiv_Kill_Partial( MortalityType code,
 void indiv_Kill_Complete( IndivType *ndv, int killType);
 void indiv_proportion_Kill( IndivType *ndv, int killType,RealF proportionKilled);
 void indiv_proportion_Recovery( IndivType *ndv, int killType,RealF proportionRecovery,RealF proportionKilled);
-void indiv_proportion_Grazing( IndivType *ndv, RealF proportionGrazing);
+RealF indiv_proportion_Grazing( IndivType *ndv, RealF proportionGrazing);
 
 /*********** Locally Used Function Declarations ************/
 /***********************************************************/
@@ -290,15 +290,15 @@ void indiv_proportion_Kill(IndivType *ndv, int killType, RealF proportKilled)
  * Implement grazing for each individual. Also keep up with survivorship data.
  * 
  * \param ndv A pointer to the individual.
- * \param proportionGrazing Value between 0 and 1. The proportion of biomass to remove.
+ * \param proportionGrazing Value between 0 and 1. The proportion of biomass to remove / reduction in relsize.
  * 
- * \sideeffect ndv->relsize is adjusted.
+ * \sideeffect ndv->relsize is adjusted and the reduction in individual relsizes due to grazing is returned.
  * 
  * \sa Species_Proportion_Grazing()
  * 
  * \ingroup INDIVIDUAL
  */
-void indiv_proportion_Grazing( IndivType *ndv, RealF proportionGrazing)
+RealF indiv_proportion_Grazing( IndivType *ndv, RealF proportionGrazing)
 {
 #define xF_DELTA (20*F_DELTA)
 #define xD_DELTA (20*D_DELTA)
@@ -319,6 +319,9 @@ void indiv_proportion_Grazing( IndivType *ndv, RealF proportionGrazing)
 #undef xF_DELTA
 #undef xD_DELTA
 #undef ZERO
+
+//return the reduction in individual relative size removed by grazing this year
+return grazing_reduce * -1;
 }
 
 /**

ST_initialization.c

@@ -83,7 +83,7 @@ void runInitialization(void){
     /* Dummy accumulators to ensure we do not collect statistics */
 	StatType *dummy_Dist, *dummy_Ppt, *dummy_Temp,
   		*dummy_Grp, *dummy_Gsize, *dummy_Gpr, *dummy_Gmort, *dummy_Gestab,
-  		*dummy_Spp, *dummy_Indv, *dummy_Smort, *dummy_Sestab, *dummy_Sreceived;
+  		*dummy_Spp, *dummy_Indv, *dummy_Smort, *dummy_Sestab, *dummy_Sreceived, *dummy_Grazed;
 	FireStatsType *dummy_Gwf;
 
 	/* For iterating over gridCells */
@@ -126,7 +126,7 @@ void runInitialization(void){
                     to accumulate stats while in spinup. We need to load in dummy accumulators to ensure
                     we ignore everything that happens in spinup. */
                 stat_Copy_Accumulators(dummy_Dist, dummy_Ppt, dummy_Temp, dummy_Grp, dummy_Gsize, dummy_Gpr, dummy_Gmort, dummy_Gestab,
-                                        dummy_Spp, dummy_Indv, dummy_Smort, dummy_Sestab, dummy_Sreceived, dummy_Gwf, TRUE);
+                                        dummy_Spp, dummy_Indv, dummy_Smort, dummy_Sestab, dummy_Sreceived, dummy_Grazed, dummy_Gwf, TRUE);
 
                 Globals->currYear = year;
 

ST_mortality.c

@@ -406,6 +406,7 @@ void grazing_EndOfYear( void){
 	{
 		IntU grazingyr =0;
 		g = RGroup[rg];
+        RGroup[rg]->res_grazed = 0;
 
 		if (Globals->currYear < RGroup[rg]->startyr)
 		{
@@ -442,7 +443,7 @@ void grazing_EndOfYear( void){
 				}
 
 				/* Remove plant biomass to implement grazing using the proportion_grazing specified in inputs */
-				Species_Proportion_Grazing(RGroup[rg]->est_spp[i],RGroup[rg]->proportion_grazing );
+				RGroup[rg]->res_grazed += Species_Proportion_Grazing(RGroup[rg]->est_spp[i],RGroup[rg]->proportion_grazing );
 			}
 		}
 	}

ST_params.c

@@ -438,7 +438,7 @@ static void _bmassflags_init( void) {
 
    FILE *fin;
    Int x, i,
-       nitems=16; /* number of items expected in first input line */
+       nitems=17; /* number of items expected in first input line */
 
    /*   code      controls: */
    char u[5],  /* summary? if 'n' don't init and don't print */
@@ -456,7 +456,8 @@ static void _bmassflags_init( void) {
         w[5],  /* wildfire count */
         m[5],  /* prescribed fire count */
         s[5],  /* biomass for each species */
-        n[5];  /* number of individuals for each species */
+        n[5],  /* number of individuals for each species */
+        b[5];  /* if 'y', output grazed biomass */
    char z;
 
    MyFileName = Parm_name(F_BMassFlag);
@@ -466,8 +467,8 @@ static void _bmassflags_init( void) {
      LogError(logfp, LOGFATAL, "%s: No data found!\n", MyFileName);
    }
 
-   x = sscanf( inbuf, "%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s",
-                      u, a, h, f, y, d, p, c, t, g, q, r, w, m, s, n );
+   x = sscanf( inbuf, "%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s",
+                      u, a, h, f, y, d, p, c, t, g, q, r, w, m, b, s, n );
 
    /* don't bother initializing the rest if first flag is 'n' */
    BmassFlags.summary  = (Bool)(*u=='y'||*u=='Y');
@@ -527,16 +528,17 @@ static void _bmassflags_init( void) {
             BmassFlags.prescribedfire     = (Bool)(*m=='y'||*m=='Y');
             break;
        case 13:
-            BmassFlags.sppb   = (Bool)(*s=='y'||*s=='Y');
+            BmassFlags.graz = (Bool)(*b == 'y' || *b == 'Y');
             break;
        case 14:
-            BmassFlags.indv   = (Bool)(*n=='y'||*n=='Y');
+            BmassFlags.sppb = (Bool)(*s == 'y' || *s == 'Y');
             break;
        case 15:
+            BmassFlags.indv = (Bool)(*n == 'y' || *n == 'Y');
             break;
      }
-
    }
+
    CloseFile(&fin);
 
    /* remove old output and/or create the output directories if needed */

ST_resgroups.c

@@ -1172,6 +1172,7 @@ void copy_rgroup(const GroupType* src, GroupType* dest){
     dest->regen_ok = src->regen_ok;
     dest->res_avail = src->res_avail;
     dest->res_extra = src->res_extra;
+    dest->res_grazed = src->res_grazed;
     dest->res_required = src->res_required;
     dest->rgroupFractionOfVegTypeBiomass = src->rgroupFractionOfVegTypeBiomass;
     dest->slowrate = src->slowrate;

ST_species.c

@@ -49,7 +49,7 @@ void indiv_Kill_Complete(IndivType *ndv, int killType);
 void indiv_proportion_Kill(IndivType *ndv, int killType, RealF proportionKilled);
 void indiv_proportion_Recovery(IndivType *ndv, int killType,
 		RealF proportionRecovery, RealF proportionKilled);
-void indiv_proportion_Grazing(IndivType *ndv, RealF proportionGrazing);
+RealF indiv_proportion_Grazing(IndivType *ndv, RealF proportionGrazing);
 
 void _delete(IndivType *ndv);
 void save_annual_species_relsize(void);
@@ -504,33 +504,35 @@ void Species_Proportion_Kill(const SppIndex sp, int killType,
 }
 
 /**
- * \brief Performs grazing on a given species at a given proportion.
+ * \brief Performs grazing on a given species at a given grazing intensity.
  * 
- * This function will perform grazing on every individual in the given species.
- * It will also graze the extra (superfluous) growth stored at the species level.
+ * This function will implement grazing on every plant individual.
+ * It will also implement grazing on the extra (superfluous) growth that is stored at the species level.
  * 
  * Initial programming by AT on 1/11/15.
  * Extra growth grazing added by Chandler Haukap 14/8/18.
+ * Saving of forage utilization for output added by Alexis Belt 5/23.
  * 
  * \param sp is the index in \ref Species of the species.
  * \param proportionGrazing is the proportion of relative size to remove. Value
  *                          between 0 and 1.
  * 
  * \sideeffect All individuals in the given species' relsizes will be reduced
- *             proportionally.\n
- *             \ref Species[sp]->relsize will be reduced by the same proportion.
+ *             proportionally and saved.\n
+ *             \ref Species[sp]->relsize will be reduced by the same proportion and saved.
  * 
  * \sa indiv_proportion_Grazing() which is the function this function calls to 
  *                                perform grazing on each individual.
  * 
  * \ingroup SPECIES
  */
-void Species_Proportion_Grazing(const SppIndex sp, RealF proportionGrazing)
+RealF Species_Proportion_Grazing(const SppIndex sp, RealF proportionGrazing)
 {
 	//CH- extra growth is only stored at the species level. This will graze extra
 	//    growth for the whole species.
-	//    loss represents the proportion of extragrowth that is eaten by livestock
+	//    loss represents the proportion of extragrowth that is removed by livestock grazing
 	RealF loss = Species[sp]->extragrowth * proportionGrazing;
+	RealF indivGrazed = 0.0;
 
 	//CH- Remove the loss from Species extra growth.
 	Species[sp]->extragrowth -= loss;	// remove the loss from extragrowth
@@ -540,9 +542,11 @@ void Species_Proportion_Grazing(const SppIndex sp, RealF proportionGrazing)
 	while (p) //while p points to an individual
 	{
 		t = p->Next; //must store Next since p might be deleted at any time.
-		indiv_proportion_Grazing(p, proportionGrazing);
+		indivGrazed += indiv_proportion_Grazing(p, proportionGrazing); //store the total reduction in individual plant size resulting from grazing
 		p = t; //move to the next plant.
 	}
+	//res_grazed = species biomass removed by livestock grazing this year
+	return Species[sp]->mature_biomass * (loss + indivGrazed);
 }
 
 /**