Merge pull request #136 from DrylandEcology/feature_VegArrays

Feature VegArrays
- vegetation type related variables are now arrays
- pre-defined macro loops for easier coding:
    * ForEachVegType
    * ForEachVegTypeBottomUp

SW_Carbon.c

@@ -209,6 +209,7 @@ void SW_CBN_read(void)
  * be simulated.
  */
 void calculate_CO2_multipliers(void) {
+  int k;
   TimeInt year,
     simendyr = SW_Model.endyr + SW_Model.addtl_yr;
   double ppm;
@@ -235,28 +236,25 @@ void calculate_CO2_multipliers(void) {
     }
 
     // Calculate multipliers per PFT
-    if (c->use_bio_mult)
-    {
-      v->grass.co2_multipliers[BIO_INDEX][year] = v->grass.co2_bio_coeff1 * pow(ppm, v->grass.co2_bio_coeff2);
-      v->shrub.co2_multipliers[BIO_INDEX][year] = v->shrub.co2_bio_coeff1 * pow(ppm, v->shrub.co2_bio_coeff2);
-      v->tree.co2_multipliers[BIO_INDEX][year] = v->tree.co2_bio_coeff1 * pow(ppm, v->tree.co2_bio_coeff2);
-      v->forb.co2_multipliers[BIO_INDEX][year] = v->forb.co2_bio_coeff1 * pow(ppm, v->forb.co2_bio_coeff2);
+    if (c->use_bio_mult) {
+      ForEachVegType(k) {
+        v->veg[k].co2_multipliers[BIO_INDEX][year] = v->veg[k].co2_bio_coeff1 * pow(ppm, v->veg[k].co2_bio_coeff2);
+      }
     }
 
     #ifdef SWDEBUG
     if (debug) {
       swprintf("Shrub: use%d: bio_mult[%d] = %1.3f / coeff1 = %1.3f / coeff2 = %1.3f / ppm = %3.2f\n",
-        c->use_bio_mult, year, v->shrub.co2_multipliers[BIO_INDEX][year],
-        v->shrub.co2_bio_coeff1, v->shrub.co2_bio_coeff2, ppm);
+        c->use_bio_mult, year, v->veg[SW_SHRUB].co2_multipliers[BIO_INDEX][year],
+        v->veg[SW_SHRUB].co2_bio_coeff1, v->veg[SW_SHRUB].co2_bio_coeff2, ppm);
     }
     #endif
 
-    if (c->use_wue_mult)
-    {
-      v->grass.co2_multipliers[WUE_INDEX][year] = v->grass.co2_wue_coeff1 * pow(ppm, v->grass.co2_wue_coeff2);
-      v->shrub.co2_multipliers[WUE_INDEX][year] = v->shrub.co2_wue_coeff1 * pow(ppm, v->shrub.co2_wue_coeff2);
-      v->tree.co2_multipliers[WUE_INDEX][year] = v->tree.co2_wue_coeff1 * pow(ppm, v->tree.co2_wue_coeff2);
-      v->forb.co2_multipliers[WUE_INDEX][year] = v->forb.co2_wue_coeff1 * pow(ppm, v->forb.co2_wue_coeff2);
+    if (c->use_wue_mult) {
+      ForEachVegType(k) {
+        v->veg[k].co2_multipliers[WUE_INDEX][year] = v->veg[k].co2_wue_coeff1 *
+          pow(ppm, v->veg[k].co2_wue_coeff2);
+      }
     }
   }
 }

SW_Defines.h

@@ -33,7 +33,6 @@
 #define MAX_LAYERS  25
 #define MAX_TRANSP_REGIONS 4
 #define MAX_ST_RGR 100
-#define NVEGTYPES 4
 
 #define MAX_NYEAR 2500  /**< An integer representing the max calendar year that is supported. The number just needs to be reasonable, it is an artifical limit. */
 
@@ -64,17 +63,27 @@
 #define SW_MIN 0
 #define SW_MAX 1
 
+/* indices to vegetation types */
+#define NVEGTYPES 4
+#define SW_TREES 0
+#define SW_SHRUB 1
+#define SW_FORBS 2
+#define SW_GRASS 3
+
+
 
 /*------------ DON'T CHANGE ANYTHING BELOW THIS LINE ------------*/
 /* Macros to simplify and add consistency to common tasks */
 /* Note the loop var must be declared as LyrIndex */
 #define ForEachSoilLayer(i)     for((i)=0; (i) < SW_Site.n_layers;      (i)++)
 #define ForEachEvapLayer(i)     for((i)=0; (i) < SW_Site.n_evap_lyrs;   (i)++)
-#define ForEachTreeTranspLayer(i)   for((i)=0; (i) < SW_Site.n_transp_lyrs_tree; (i)++)
-#define ForEachShrubTranspLayer(i)   for((i)=0; (i) < SW_Site.n_transp_lyrs_shrub; (i)++)
-#define ForEachGrassTranspLayer(i)   for((i)=0; (i) < SW_Site.n_transp_lyrs_grass; (i)++)
-#define ForEachForbTranspLayer(i)   for((i)=0; (i) < SW_Site.n_transp_lyrs_forb; (i)++)
+#define ForEachTreeTranspLayer(i)   for((i)=0; (i) < SW_Site.n_transp_lyrs[SW_TREES]; (i)++)
+#define ForEachShrubTranspLayer(i)   for((i)=0; (i) < SW_Site.n_transp_lyrs[SW_SHRUB]; (i)++)
+#define ForEachGrassTranspLayer(i)   for((i)=0; (i) < SW_Site.n_transp_lyrs[SW_GRASS]; (i)++)
+#define ForEachForbTranspLayer(i)   for((i)=0; (i) < SW_Site.n_transp_lyrs[SW_FORBS]; (i)++)
 #define ForEachTranspRegion(r)  for((r)=0; (r) < SW_Site.n_transp_rgn;  (r)++)
+#define ForEachVegType(k)  for ((k) = 0; (k) < NVEGTYPES; (k)++)
+#define ForEachVegTypeBottomUp(k)  for ((k) = NVEGTYPES - 1; (k) >= 0; (k)--)
 /* define m as Months */
 #define ForEachMonth(m)         for((m)=Jan; (m) <= Dec;  (m)++)
 

SW_Flow_lib.c

@@ -120,24 +120,24 @@ static ST_RGR_VALUES stValues; // keeps track of the soil_temperature values
 /* --------------------------------------------------- */
 
 /**
-	\fn void grass_intercepted_water(double *pptleft, double *wintgrass, double ppt, double vegcov, double scale, double a, double b, double c, double d) {
-	\brief Calculate the water intercepted by grasses.
+	\fn void veg_intercepted_water(double *pptleft, double *wintgrass, double ppt, double x, double scale, double a, double b, double c, double d) {
+	\brief Calculate the water intercepted by vegetation.
 
 	Equations based on Corbet and Crouse (1968). \cite Corbett1968
 
-  \param pptleft.
-	\param wintgrass.
-	\param ppt. daily precipitation
-	\param vegcov. vegetation cover of grass for the day (based on monthly biomass
-  values, see the routine "initprod").
-	\param scale. Scale paramater. Used to represent snow depth.
-	\param a. a parameter for intercept of grass interception equation.
-	\param b. b parameter for intercept of grass interception equation.
-	\param c. c parameter for slope of grass interception equation.
-	\param d. d parameter for slope of grass interception equation.
-
-	\return pptleft. Amount of precipitation left after interception.
-	\return wintgrass. Amount of precipitation interception by grass.
+  \param pptleft
+	\param wintgrass
+	\param ppt daily precipitation
+	\param x vegetation cover or LAI for the day (based on monthly biomass
+	       values, see the routine "initprod").
+	\param scale Scale paramater. Used to represent snow depth.
+	\param a a parameter for intercept of grass interception equation.
+	\param b b parameter for intercept of grass interception equation.
+	\param c c parameter for slope of grass interception equation.
+	\param d d parameter for slope of grass interception equation.
+
+	\return pptleft Amount of precipitation left after interception.
+	\return wintgrass Amount of precipitation interception by grass.
 */
 
 
@@ -158,112 +158,25 @@ static ST_RGR_VALUES stValues; // keeps track of the soil_temperature values
 
  21-Oct-03 (cwb) added MAX_WINTLIT line
  **********************************************************************/
-void grass_intercepted_water(double *pptleft, double *wintgrass, double ppt, double vegcov, double scale, double a, double b, double c, double d) {
-	double intcpt, slope;
-
-	if (GT(vegcov, 0.) && GT(ppt, 0.)) {
-		intcpt = b * vegcov + a;
-		slope = d * vegcov + c;
-
-		*wintgrass = (intcpt + slope * ppt) * scale;
-
-		*wintgrass = fmin(*wintgrass, ppt);
-		*wintgrass = fmin(*wintgrass, MAX_WINTSTCR);
-		*pptleft = fmax( ppt - *wintgrass, 0.0);
-	} else { /*  no precip., so obviously nothing is intercepted by standing crop. */
-		*pptleft = ppt;
-		*wintgrass = 0.0;
-	}
-}
-
-void shrub_intercepted_water(double *pptleft, double *wintshrub, double ppt, double vegcov, double scale, double a, double b, double c, double d) {
-	/**********************************************************************
-	 PURPOSE: Calculate the water intercepted by shrubs
-	 **********************************************************************/
-	double intcpt, slope;
-
-	if (GT(vegcov, 0.) && GT(ppt, 0.)) {
-		intcpt = b * vegcov + a;
-		slope = d * vegcov + c;
-
-		*wintshrub = (intcpt + slope * ppt) * scale;
-
-		*wintshrub = fmin(*wintshrub, ppt);
-		*wintshrub = fmin(*wintshrub, MAX_WINTSTCR);
-		*pptleft = fmax( ppt - *wintshrub, 0.0);
-	} else { /*  no precip., so obviously nothing is intercepted by standing crop. */
-		*pptleft = ppt;
-		*wintshrub = 0.0;
-	}
-}
-
-void tree_intercepted_water(double *pptleft, double *wintfor, double ppt, double LAI, double scale, double a, double b, double c, double d) {
-	/**********************************************************************
-	 PURPOSE: Calculate water intercepted by forests
-
-	 HISTORY:
-	 11/16/2010	(drs)
-
-	 INPUTS:
-	 ppt     - precip. for the day in cm
-	 LAI	- forest LAI in cm/cm
-	 scale - scale interception with fraction of tree vegetation component or with snowdepth-scaler
-
-	 OUTPUT:
-	 pptleft -  precip. left after interception by forest in cm.
-	 wintfor - amount of water intercepted by forest in cm.
-	 **********************************************************************/
+void veg_intercepted_water(double *pptleft, double *wintveg, double ppt, double x,
+  double scale, double a, double b, double c, double d)
+{
 	double intcpt, slope;
 
-	if (GT(LAI, 0.) && GT(ppt, 0.)) {
-		intcpt = b * LAI + a;
-		slope = d * LAI + c;
-
-		*wintfor = (intcpt + slope * ppt) * scale;
-
-		*wintfor = fmin(*wintfor, ppt);
-		*wintfor = fmin(*wintfor, MAX_WINTFOR);
-		*pptleft = fmax( ppt - *wintfor, 0.0);
-	} else { /*  no precip., so obviously nothing is intercepted by forest. */
-		*pptleft = ppt;
-		*wintfor = 0.0;
-	}
-}
-
-void forb_intercepted_water(double *pptleft, double *wintforb, double ppt, double vegcov, double scale, double a, double b, double c, double d) {
-	/**********************************************************************
-	 PURPOSE: Calculate water intercepted by forbs
-
-
-	 HISTORY:
-	 07/09/2013	(clk)
-
-	 INPUTS:
-
-	 ppt     - precip. for the day in cm
-	 vegcov	- vegetation cover for the day (based on monthly biomass
-	 values, see the routine "initprod")
-	 scale - scale interception with fraction of forb vegetation component or with snowdepth-scaler
-
-	 OUTPUT:
-
-	 pptleft -  precip. left after interception by forb in cm.
-	 wintforb - amount of water intercepted by forb in cm.
-	 **********************************************************************/
-	double intcpt, slope;
+	if (GT(x, 0.) && GT(ppt, 0.)) {
+		intcpt = b * x + a;
+		slope = d * x + c;
 
-	if (GT(vegcov, 0.) && GT(ppt, 0.)) {
-		intcpt = b * vegcov + a;
-		slope = d * vegcov + c;
+		*wintveg = (intcpt + slope * ppt) * scale;
 
-		*wintforb = (intcpt + slope * ppt) * scale;
+		*wintveg = fmin(*wintveg, ppt);
+		*wintveg = fmin(*wintveg, MAX_WINTSTCR(x));
+		*pptleft = fmax(ppt - *wintveg, 0.0);
 
-		*wintforb = fmin(*wintforb, ppt);
-		*wintforb = fmin(*wintforb, MAX_WINTSTCR);
-		*pptleft = fmax( ppt - *wintforb, 0.0);
-	} else { /*  no precip., so obviously nothing is intercepted by forest. */
+	} else {
+	  /*  no precip. or no biomass, so obviously nothing is intercepted by standing crop. */
 		*pptleft = ppt;
-		*wintforb = 0.0;
+		*wintveg = 0.0;
 	}
 }
 
@@ -636,7 +549,7 @@ void transp_weighted_avg(double *swp_avg, unsigned int n_tr_rgns, unsigned int n
 	* @param blivelai the live biomass leaf area index
 	* @param lai_param
 	*/
-void grass_EsT_partitioning(double *fbse, double *fbst, double blivelai, double lai_param) {
+void EsT_partitioning(double *fbse, double *fbst, double blivelai, double lai_param) {
 	/**********************************************************************
 	 PURPOSE: Calculate fraction of water loss from bare soil
 	 evaporation and transpiration
@@ -645,6 +558,10 @@ void grass_EsT_partitioning(double *fbse, double *fbst, double blivelai, double
 	 4/30/92  (SLC)
 	 24-Oct-03 (cwb) changed exp(-blivelai*bsepar1) + bsepar2;
 	 to exp(-blivelai);
+	 08/22/2011	(drs)	For trees: according to a regression based on a review by
+	  Daikoku, K., S. Hattori, A. Deguchi, Y. Aoki, M. Miyashita, K. Matsumoto, J. Akiyama,
+	  S. Iida, T. Toba, Y. Fujita, and T. Ohta. 2008. Influence of evaporation from the
+	  forest floor on evapotranspiration from the dry canopy. Hydrological Processes 22:4083-4096.
 
 	 INPUTS:
 	 blivelai - live biomass leaf area index
@@ -666,50 +583,6 @@ void grass_EsT_partitioning(double *fbse, double *fbst, double blivelai, double
 	*fbst = 1. - (*fbse);
 }
 
-void shrub_EsT_partitioning(double *fbse, double *fbst, double blivelai, double lai_param) {
-	/**********************************************************************
-	 PURPOSE: Calculate fraction of water loss from bare soil
-	 evaporation and transpiration
-	 **********************************************************************/
-
-	double bsemax = 0.995;
-
-	*fbse = exp(-lai_param * blivelai);
-
-	*fbse = fmin(*fbse, bsemax);
-	*fbst = 1. - (*fbse);
-}
-
-void tree_EsT_partitioning(double *fbse, double *fbst, double blivelai, double lai_param) {
-	/**********************************************************************
-	 PURPOSE: Calculate fraction of water loss from bare soil
-	 evaporation and transpiration
-
-	 08/22/2011	(drs)	According to a regression based on a review by Daikoku, K., S. Hattori, A. Deguchi, Y. Aoki, M. Miyashita, K. Matsumoto, J. Akiyama, S. Iida, T. Toba, Y. Fujita, and T. Ohta. 2008. Influence of evaporation from the forest floor on evapotranspiration from the dry canopy. Hydrological Processes 22:4083-4096.
-	 **********************************************************************/
-
-	double bsemax = 0.995;
-
-	*fbse = exp(-lai_param * blivelai);
-
-	*fbse = fmin(*fbse, bsemax);
-	*fbst = 1. - (*fbse);
-}
-
-void forb_EsT_partitioning(double *fbse, double *fbst, double blivelai, double lai_param) {
-	/**********************************************************************
-	 PURPOSE: Calculate fraction of water loss from bare soil
-	 evaporation and transpiration
-
-	 **********************************************************************/
-
-	double bsemax = 0.995;
-
-	*fbse = exp(-lai_param * blivelai);
-
-	*fbse = fmin(*fbse, bsemax);
-	*fbst = 1. - (*fbse);
-}
 
 void pot_soil_evap(double *bserate, unsigned int nelyrs, double ecoeff[], double totagb, double fbse, double petday, double shift, double shape, double inflec, double range,
 		double width[], double swc[], double Es_param_limit) {

SW_Flow_lib.h

@@ -40,7 +40,7 @@
 /* Standing crop can only intercept so much precip
  * This is the limiter used inside stdcrop_intercepted()
  */
-#define MAX_WINTSTCR (vegcov * .1)
+#define MAX_WINTSTCR(x) ((x) * .1)
 #define MAX_WINTFOR (ppt)
 
 /* Litter can only intercept so much precip.
@@ -80,13 +80,9 @@ typedef struct {
 /* =================================================== */
 /*                Function Definitions                 */
 /* --------------------------------------------------- */
-void grass_intercepted_water(double *pptleft, double *wintgrass, double ppt, double vegcov, double scale, double a, double b, double c, double d);
 
-void shrub_intercepted_water(double *pptleft, double *wintshrub, double ppt, double vegcov, double scale, double a, double b, double c, double d);
-
-void tree_intercepted_water(double *pptleft, double *wintfor, double ppt, double LAI, double scale, double a, double b, double c, double d);
-
-void forb_intercepted_water(double *pptleft, double *wintforb, double ppt, double vegcov, double scale, double a, double b, double c, double d);
+void veg_intercepted_water(double *pptleft, double *wintveg, double ppt, double x,
+  double scale, double a, double b, double c, double d);
 
 void litter_intercepted_water(double *pptleft, double *wintlit, double blitter, double scale, double a, double b, double c, double d);
 
@@ -100,10 +96,7 @@ double svapor(double temp);
 
 void transp_weighted_avg(double *swp_avg, unsigned int n_tr_rgns, unsigned int n_layers, unsigned int tr_regions[], double tr_coeff[], double swc[]);
 
-void grass_EsT_partitioning(double *fbse, double *fbst, double blivelai, double lai_param);
-void shrub_EsT_partitioning(double *fbse, double *fbst, double blivelai, double lai_param);
-void tree_EsT_partitioning(double *fbse, double *fbst, double blivelai, double lai_param);
-void forb_EsT_partitioning(double *fbse, double *fbst, double blivelai, double lai_param);
+void EsT_partitioning(double *fbse, double *fbst, double blivelai, double lai_param);
 
 void pot_soil_evap(double *bserate, unsigned int nelyrs, double ecoeff[], double totagb, double fbse, double petday, double shift, double shape, double inflec, double range,
 		double width[], double swc[], double Es_param_limit);