Add a constructor for the SW_CARBON structure, remove legacy defensiv…

…e code

Close #29: CO2 multipliers incorrect when scenario doesn't exist
Close #34: Properly initialize the CO2 variables

SW_Carbon.c

@@ -44,6 +44,28 @@ SW_VEGPROD SW_VegProd;  // Declared here, externed elsewhere
 /*             Public Function Definitions             */
 /* --------------------------------------------------- */
 
+/**
+ * Initializes the multiplier values in the SW_CARBON structure
+ * @note The spin-up year has been known to have the multipliers equal
+ *       to 0 without this constructor
+ */
+void SW_CBN_construct(void)
+{
+	memset(&SW_Carbon, 0, sizeof(SW_Carbon));
+
+	SW_CARBON *c = &SW_Carbon;
+	int i;
+
+	for (i = 0; i < MAX_CO2_YEAR; i++)
+	{
+		c->co2_multipliers[0][i] = 1.0;
+		c->co2_multipliers[1][i] = 1.0;
+	}
+
+	c->co2_biomass_mult = 1.0;
+	c->co2_wue_mult = 1.0;
+}
+
 /* A description on how these 'onGet' and 'onSet' functions work...
  * Summary: onGet instantiates the 'swCarbon' class and returns the object, and is only used once
  *          onSet extracts the value of the given object, and is used on both calls to SOILWAT2
@@ -79,15 +101,15 @@ void onSet_swCarbon(SEXP object) {
   strcpy(c->scenario, CHAR(STRING_ELT(GET_SLOT(object, install("Scenario")), 0)));  // e.g. c->scenario = "RCP85"
 }
 #endif
-
-
+
 /**
  * Calculates the multipliers for biomass and the transpiration
  * equation. If a multiplier is disabled, its value is set to 1. All the years
  * in carbon.in are calculated. RSOILWAT will pass in the settings directly.
  * SOILWAT will read siteparam.in for settings.
  */
 void calculate_CO2_multipliers(void) {
+  int i;
   FILE *f;
   char scenario[64];
   double ppm;
@@ -98,32 +120,27 @@ void calculate_CO2_multipliers(void) {
   MyFileName     = SW_F_name(eCarbon);
   f              = OpenFile(MyFileName, "r");
 
-  // Read carbon.in
+  /* Calculate multipliers by reading carbon.in */
   while (GetALine(f, inbuf)) {
-	// Scan for the year first, because if the year is 0 it marks a change in the scenario
+	// Read the year standalone because if it's 0 it marks a change in the scenario,
+	// in which case we'll need to read in a string instead of an int
 	sscanf(inbuf, "%d", &year);
-	  
-	// We found a scenario, do we want this one?
+
+	// Find scenario
 	if (year == 0)
 	{
       sscanf(inbuf, "%d %63s", &year, scenario);
-	  continue;
+	  continue;  // Skip to the PPM values
 	}
-
-	// NO, keep searching
 	if (strcmp(scenario, c->scenario) != 0)
-	  continue;
+	  continue;  // Keep searching for the right scenario
 
-	// YES, calculate multipliers
-	sscanf(inbuf, "%d %lf", &year, &ppm);
-    c->co2_multipliers[0][year] = 1.0;
-    c->co2_multipliers[1][year] = 1.0;
+    sscanf(inbuf, "%d %lf", &year, &ppm);
 	if (c->use_bio_mult) c->co2_multipliers[0][year] = v->co2_biomass_1  * pow(ppm, v->co2_biomass_2);
 	if (c->use_sto_mult) c->co2_multipliers[1][year] = v->co2_stomatal_1 * pow(ppm, v->co2_stomatal_2);
   }
 }
 
-
 /**
  * Applies CO2 effects to supplied biomass data. Two parameters are needed so that
  * we do not have a compound effect on the biomass.
@@ -134,7 +151,6 @@ void apply_CO2(double new_biomass[], double biomass[]) {
   int i;
   SW_CARBON  *c  = &SW_Carbon;
 
-  if (c->co2_biomass_mult == 0) c->co2_biomass_mult = 1.0;  // The start-up year has a multiplier of 0 since it has not calculated the multipliers yet
   for (i = 0; i < 12; i++) {
     new_biomass[i] = (biomass[i] * c->co2_biomass_mult);
   }

SW_Carbon.h

@@ -8,21 +8,22 @@
  */
 #ifndef CARBON
   #define CARBON
-
-  typedef struct {
+  #define MAX_CO2_YEAR 3000  // Include all years that our PPM sources might provide data for
+
+  typedef struct {	
     int
-    addtl_yr,                  // Added to SW_Model.year to get the future year we're simulating
-    use_sto_mult,              // Determine which multipliers we will be calculating...
+    addtl_yr,                          // Added to SW_Model.year to get the future year we're simulating
+    use_sto_mult,                      // Determine which multipliers we will be calculating...
     use_bio_mult;
 
     double
-    carbon[2],                 // Hold misc. data
-    co2_biomass_mult,          // The biomass multiplier (yearly)
-    co2_wue_mult,              // The stomatal multiplier (yearly)
-    co2_multipliers[2][3000];  // Holds the above multipliers for every year, accessed directly (e.g. biomass multiplier for 1982 is co2_multipliers[1][1982])
-
+    carbon[2],                         // Hold misc. data
+    co2_multipliers[2][MAX_CO2_YEAR],  // Holds the above multipliers for every year, accessed directly (e.g. biomass multiplier for 1982 is co2_multipliers[1][1982])
+	co2_biomass_mult,                  // The biomass multiplier (yearly)
+    co2_wue_mult;                      // The stomatal multiplier (yearly)
+	
 	char
-	scenario[64];              // The scenario we are extracting PPM from
+	scenario[64];                      // The scenario we are extracting PPM from
   } SW_CARBON;
 
   /* Function Declarations */
@@ -31,7 +32,7 @@
     void onSet_swCarbon(SEXP object);
 
   #endif
-
+  void SW_CBN_construct(void);
   void apply_CO2(double* new_biomass, double *biomass);
   void calculate_CO2_multipliers(void);
 

SW_Control.c

@@ -95,6 +95,7 @@ void SW_CTL_init_model(const char *firstfile) {
 	SW_OUT_construct();
 	SW_SWC_construct();
 	SW_FLW_construct();
+	SW_CBN_construct();
 
 	_read_inputs();
 
@@ -200,8 +201,8 @@ static void _read_inputs(void) {
 		//Rprintf("swOutput\n");
 		onSet_SW_SWC(GET_SLOT(InputData,install("swc")));
 		//Rprintf("swSWC\n");
-    onSet_swCarbon(GET_SLOT(InputData, install("carbon")));
-    //Rprintf("swCarbon\n");
+		onSet_swCarbon(GET_SLOT(InputData, install("carbon")));
+		//Rprintf("swCarbon\n");
 	}
 #endif
   // This will work in both RSOILWAT and SOILWAT

SW_Flow_lib.c

@@ -853,8 +853,6 @@ void pot_transp(double *bstrate, double swpavg, double biolive, double biodead,
 			shadeaf = 1.0;
 		}
 
-		if (c->co2_wue_mult == 0) c->co2_wue_mult = 1.0;  // The start-up year has a multiplier of 0 since it has not calculated the multipliers yet
-
 		*bstrate = watrate(swpavg, petday, swp_shift, swp_shape,
 							   swp_inflec, swp_range) * shadeaf * petday * fbst * c->co2_wue_mult;
 	}

SW_Output.c

@@ -1209,7 +1209,7 @@ void SW_OUT_write_today(void)
 	int i;
 
 	// Adjust the model year to match simulation years for the output
-  SW_Model.year += SW_Carbon.addtl_yr;
+    SW_Model.year += SW_Carbon.addtl_yr;
 
 	ForEachOutKey(k)
 	{