Merge pull request #322 from DrylandEcology/feature_swrc

Feature swrc

* close #315 

Multiple soil water release curves (`SWRC`) are now implemented and can be
selected with new input `swrc_name`. Implemented `SWRCs` currently include
`"Campbell1974"`, `"vanGenuchten1980"`, and `"FXW"`. New input `has_swrcp`
determines if parameters for a `SWRC` are estimated at run-time via an
implemented pedotransfer function (`PTF`) based on new input `ptf_name` or
if they are provided as inputs via new input file `"swrc_params.in"`.
`rSOILWAT2` implements additional pedotransfer functions. See documentation
entry of `"swrc_ptf"` for additional details and for guidance on how to
implement additional `SWRCs` and `PTFs` (issue #315; @dschlaep).

.LSAN_suppr.txt

@@ -5,3 +5,6 @@
 
 # https://github.com/DrylandEcology/SOILWAT2/issues/205
 leak:SW_VPD_construct
+
+# on Apple arm64: https://github.com/google/sanitizers/issues/1501
+leak:realizeClassWithoutSwift

NEWS.md

@@ -1,6 +1,21 @@
 # NEWS
 
 # SOILWAT2 v7.0.0-9000
+* This version produces nearly identical simulation output
+  as the previous release under default values for the new inputs.
+  Small deviations arise due to a fix in the handling of soil moisture values
+  when between field capacity and saturation.
+
+* Multiple soil water release curves (`SWRC`) are now implemented and can be
+  selected with new input `swrc_name`. Implemented `SWRCs` currently include
+  `"Campbell1974"`, `"vanGenuchten1980"`, and `"FXW"`. New input `has_swrcp`
+  determines if parameters for a `SWRC` are estimated at run-time via an
+  implemented pedotransfer function (`PTF`) based on new input `ptf_name` or
+  if they are provided as inputs via new input file `"swrc_params.in"`.
+  `rSOILWAT2` implements additional pedotransfer functions. See documentation
+  entry of `"swrc_ptf"` for additional details and for guidance on how to
+  implement additional `SWRCs` and `PTFs` (issue #315; @dschlaep).
+
 * Soil density inputs can now represent either matric or bulk density
   (issue #280; @dschlaep).
     * Automatic conversion between matric and bulk density as needed
@@ -56,6 +71,15 @@
 
 
 ## Changes to inputs
+* New inputs via `"siteparam.in"` select a soil water release curve `swrc_name`
+  and determine parameter source `has_swrcp`, i.e.,
+  either estimated via selected pedotransfer function `ptf_name` or
+  values obtained from new input file `"swrc_params.in"`.
+  Default values `"Campbell1974"`, `"Cosby1984AndOthers"`, and 0
+  (i.e., use `PTF` to estimate paramaters) reproduce previous behavior.
+* New input file `"swrc_params.in"` to provide parameters of the selected
+  soil water release curve (if not estimated via a pedotransfer function)
+  for each soil layer.
 * SOILWAT2 gains `type_soilDensityInput` as new user input (`siteparam.in`)
   with default value 0 (i.e., matric soil density)
   that reproduces previous behavior.

SW_Control.c

@@ -275,9 +275,19 @@ void SW_CTL_read_inputs_from_disk(void) {
   if (debug) swprintf(" > 'veg'");
   #endif
 
-  SW_SIT_read(); // inputs also soil layer data
+  SW_SIT_read();
   #ifdef SWDEBUG
-  if (debug) swprintf(" > 'site' + 'soils'");
+  if (debug) swprintf(" > 'site'");
+  #endif
+
+  SW_LYR_read();
+  #ifdef RSWDEBUG
+  if (debug) swprintf(" > 'soils'");
+  #endif
+
+  SW_SWRC_read();
+  #ifdef SWDEBUG
+  if (debug) swprintf(" > 'swrc parameters'");
   #endif
 
   SW_VES_read();

SW_Defines.h

@@ -51,13 +51,21 @@ extern "C" {
 #define SLOW_DRAIN_DEPTH 15. /* numerator over depth in slow drain equation */
 
 /* some basic constants */
-#define MAX_LAYERS  25
-#define MAX_TRANSP_REGIONS 4
-#define MAX_ST_RGR 100
+#define MAX_LAYERS  25 /**< Maximum number of soil layers */
+#define MAX_TRANSP_REGIONS 4 /**< Maximum number of transpiration regions */
+#define MAX_ST_RGR 100 /**< Maximum number of soil temperature nodes */
 
 #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. */
 
-#define SW_MISSING     999.     /* value to use as MISSING */
+#define SW_MISSING     999.     /**< Value to use as MISSING */
+
+
+// Euler's constant
+#ifdef M_E
+  #define swE M_E
+#else
+  #define swE         2.71828182845904523536028747135266249
+#endif
 
 
 /* M_PI and M_PI_2 from <math.h> if implementation conforms to POSIX extension
@@ -80,7 +88,6 @@ extern "C" {
 #define deg_to_rad    0.0174532925199433 /**< Convert arc-degrees to radians, i.e., x * deg_to_rad with deg_to_rad = pi / 180 */
 #define rad_to_deg    57.29577951308232 /**< Convert radians to arc-degrees, i.e., x * rad_to_deg with rad_to_deg = 180 / pi */
 
-#define BARCONV     1024.
 #define SEC_PER_DAY	86400. // the # of seconds in a day... (24 hrs * 60 mins/hr * 60 sec/min = 86400 seconds)
 
 
@@ -104,7 +111,7 @@ extern "C" {
 #define SW_MAX 1
 
 /* indices to vegetation types */
-#define NVEGTYPES 4
+#define NVEGTYPES 4 /**< Number of vegetation types implemented */
 #define SW_TREES 0
 #define SW_SHRUB 1
 #define SW_FORBS 2

SW_Files.c

@@ -148,53 +148,53 @@ void SW_F_read(const char *s) {
     #endif
 
 		switch (lineno) {
-		case 5:
+		case 6:
 			strcpy(weather_prefix, inbuf);
 			break;
-		case 13:
+		case 14:
 			strcpy(output_prefix, inbuf);
 			break;
-		case 15:
+		case 16:
 			InFiles[eOutputDaily] = Str_Dup(inbuf);
 			++fileno;
 			SW_CSV_F_INIT(InFiles[eOutputDaily]);
 			break;
-		case 16:
+		case 17:
 			InFiles[eOutputWeekly] = Str_Dup(inbuf);
 			++fileno;
 			SW_CSV_F_INIT(InFiles[eOutputWeekly]);
 			//printf("filename: %s \n",InFiles[eOutputWeekly]);
 			break;
-		case 17:
+		case 18:
 			InFiles[eOutputMonthly] = Str_Dup(inbuf);
 			++fileno;
 			SW_CSV_F_INIT(InFiles[eOutputMonthly]);
 			//printf("filename: %s \n",InFiles[eOutputMonthly]);
 			break;
-		case 18:
+		case 19:
 			InFiles[eOutputYearly] = Str_Dup(inbuf);
 			++fileno;
 			SW_CSV_F_INIT(InFiles[eOutputYearly]);
 			break;
-		case 19:
+		case 20:
 			InFiles[eOutputDaily_soil] = Str_Dup(inbuf);
 			++fileno;
 			SW_CSV_F_INIT(InFiles[eOutputDaily_soil]);
 			//printf("filename: %s \n",InFiles[eOutputDaily]);
 			break;
-		case 20:
+		case 21:
 			InFiles[eOutputWeekly_soil] = Str_Dup(inbuf);
 			++fileno;
 			SW_CSV_F_INIT(InFiles[eOutputWeekly_soil]);
 			//printf("filename: %s \n",InFiles[eOutputWeekly]);
 			break;
-		case 21:
+		case 22:
 			InFiles[eOutputMonthly_soil] = Str_Dup(inbuf);
 			++fileno;
 			SW_CSV_F_INIT(InFiles[eOutputMonthly_soil]);
 			//printf("filename: %s \n",InFiles[eOutputMonthly]);
 			break;
-		case 22:
+		case 23:
 			InFiles[eOutputYearly_soil] = Str_Dup(inbuf);
 			++fileno;
 			SW_CSV_F_INIT(InFiles[eOutputYearly_soil]);

SW_Files.h

@@ -20,16 +20,30 @@
 extern "C" {
 #endif
 
-#define SW_NFILES 22
+#define SW_NFILES 23
 
 /* The number of enum elements between eNoFile and
  * eEndFile (not inclusive) must match SW_NFILES.
  * also, these elements must match the order of
- * input from files.in.
+ * input from `files.in`.
  */
 typedef enum {
-	eNoFile = -1, eFirst = 0, eModel, eLog, eSite, eLayers, eWeather,
-	eMarkovProb, eMarkovCov, eSky, eVegProd, eVegEstab, eCarbon, eSoilwat,
+	eNoFile = -1,
+	/* List of all input files */
+	eFirst = 0,
+	/* Description of a model run */
+	eModel, eLog,
+	/* Description of simulated site */
+	eSite, eLayers, eSWRCp,
+	/* Weather and climate forcing */
+	eWeather, eMarkovProb, eMarkovCov, eSky,
+	/* Description of vegetation */
+	eVegProd, eVegEstab,
+	/* Description of CO2 effects */
+	eCarbon,
+	/* (optional) soil moisture measurements */
+	eSoilwat,
+	/* Simulation outputs */
 	eOutput, eOutputDaily, eOutputWeekly, eOutputMonthly, eOutputYearly,
 	eOutputDaily_soil, eOutputWeekly_soil, eOutputMonthly_soil, eOutputYearly_soil,
 	eEndFile

SW_Flow_lib.c

@@ -372,7 +372,7 @@ void transp_weighted_avg(double *swp_avg, unsigned int n_tr_rgns, unsigned int n
 
 		for (i = 0; i < n_layers; i++) {
 			if (tr_regions[i] == r) {
-				swp += tr_coeff[i] * SW_SWCbulk2SWPmatric(SW_Site.lyr[i]->fractionVolBulk_gravel, swc[i], i);
+				swp += tr_coeff[i] * SW_SWRC_SWCtoSWP(swc[i], SW_Site.lyr[i]);
 				sumco += tr_coeff[i];
 			}
 		}
@@ -484,10 +484,10 @@ void pot_soil_evap(double *bserate, unsigned int nelyrs, double ecoeff[], double
 	  }
 		x = width[i] * ecoeff[i];
 		sumwidth += x;
-		avswp += x * SW_SWCbulk2SWPmatric(SW_Site.lyr[i]->fractionVolBulk_gravel, swc[i], i);
+		avswp += x * SW_SWRC_SWCtoSWP(swc[i], SW_Site.lyr[i]);
 	}
 
-  // Note: avswp = 0 if swc = 0 because that is the return value of SW_SWCbulk2SWPmatric
+  // Note: avswp = 0 if swc = 0 because that is the return value of SW_SWRC_SWCtoSWP
 	avswp /= (ZRO(sumwidth)) ? 1 : sumwidth;
 
 	/*  8/27/92 (SLC) if totagb > Es_param_limit, assume soil surface is
@@ -542,7 +542,7 @@ void pot_soil_evap_bs(double *bserate, unsigned int nelyrs, double ecoeff[], dou
 	for (i = 0; i < nelyrs; i++) {
 		x = width[i] * ecoeff[i];
 		sumwidth += x;
-		avswp += x * SW_SWCbulk2SWPmatric(SW_Site.lyr[i]->fractionVolBulk_gravel, swc[i], i);
+		avswp += x * SW_SWRC_SWCtoSWP(swc[i], SW_Site.lyr[i]);
 	}
 
 	avswp /= sumwidth;
@@ -742,12 +742,18 @@ void remove_from_soil(double swc[], double qty[], double *aet, unsigned int nlyr
 	 **********************************************************************/
 
 	unsigned int i;
-	double swpfrac[MAX_LAYERS], sumswp = 0.0, swc_avail, q;
+	double swpfrac[MAX_LAYERS], sumswp = 0.0, swc_avail, q, tmpswp;
 
     SW_SOILWAT *st = &SW_Soilwat;
 
 	for (i = 0; i < nlyrs; i++) {
-		swpfrac[i] = coeff[i] / SW_SWCbulk2SWPmatric(SW_Site.lyr[i]->fractionVolBulk_gravel, swc[i], i);
+		tmpswp = SW_SWRC_SWCtoSWP(swc[i], SW_Site.lyr[i]);
+		if (GT(tmpswp, 0.)) {
+			swpfrac[i] = coeff[i] / tmpswp;
+		} else {
+			// saturated conditions -> divide by SWP [-bar] at field capacity
+			swpfrac[i] = coeff[i] / 0.333;
+		}
 		sumswp += swpfrac[i];
 	}
 
@@ -969,7 +975,7 @@ void hydraulic_redistribution(
 			swc[i] - fmin(lyr[i]->swcBulk_wiltpt, lyr[i]->swcBulk_atSWPcrit[vegk])
 		);
 
-		swp[i] = SW_SWCbulk2SWPmatric(lyr[i]->fractionVolBulk_gravel, swc[i], i);
+		swp[i] = SW_SWRC_SWCtoSWP(swc[i], SW_Site.lyr[i]);
 
 		/* Ryel et al. 2002: eq. 7 relative soil-root conductance */
 		relCondroot[i] = fmin(

SW_Output.h

@@ -234,6 +234,7 @@ extern IntUS ncol_OUT[SW_OUTNKEYS];
 
 extern char const *key2str[];
 extern char const *pd2longstr[];
+extern char const *styp2str[];
 
 
 

SW_Output_get_functions.c

@@ -979,8 +979,8 @@ void get_swpMatric_text(OutPeriod pd)
 	ForEachSoilLayer(i)
 	{
 		/* swpMatric at this point is identical to swcBulk */
-		val = SW_SWCbulk2SWPmatric(SW_Site.lyr[i]->fractionVolBulk_gravel,
-			vo->swpMatric[i], i);
+		val = SW_SWRC_SWCtoSWP(vo->swpMatric[i], SW_Site.lyr[i]);
+
 
 		snprintf(str, OUTSTRLEN, "%c%.*f", _Sep, OUT_DIGITS, val);
 		strcat(sw_outstr, str);
@@ -1006,8 +1006,7 @@ void get_swpMatric_mem(OutPeriod pd)
 	ForEachSoilLayer(i)
 	{
 		/* swpMatric at this point is identical to swcBulk */
-		p[iOUT(i, pd)] = SW_SWCbulk2SWPmatric(
-			SW_Site.lyr[i]->fractionVolBulk_gravel, vo->swpMatric[i], i);
+		p[iOUT(i, pd)] = SW_SWRC_SWCtoSWP(vo->swpMatric[i], SW_Site.lyr[i]);
 	}
 }
 
@@ -1031,8 +1030,7 @@ void get_swpMatric_agg(OutPeriod pd)
 	ForEachSoilLayer(i)
 	{
 		/* swpMatric at this point is identical to swcBulk */
-		val = SW_SWCbulk2SWPmatric(
-			SW_Site.lyr[i]->fractionVolBulk_gravel, vo->swpMatric[i], i);
+		val = SW_SWRC_SWCtoSWP(vo->swpMatric[i], SW_Site.lyr[i]);
 		do_running_agg(p, psd, iOUT(i, pd), Globals->currIter, val);
 	}