Merge pull request #321 from DrylandEcology/feature_veg_estimation

Feature veg estimation, climate predictors, and read all weather

This merge combines
- feature_read_weather (close #311): read and process daily weather before main simulation loop
- feature_climate_predictors (close #317): calculate long-term climate summaries
- feature_veg_estimation (close #318): estimate fractional land cover
  representing a potential natural vegetation based on climate relationships

.github/workflows/check_doc.yml

@@ -13,7 +13,7 @@ jobs:
 
     steps:
     - name: Checkout repository and submodules
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
       with:
         submodules: recursive
 

.github/workflows/main_nix.yml

@@ -20,7 +20,7 @@ jobs:
 
     steps:
     - name: Checkout repository and submodules
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
       with:
         submodules: recursive
 
@@ -50,14 +50,14 @@ jobs:
 
     steps:
     - name: Checkout repository and submodules
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
       with:
         submodules: recursive
 
     - name: Generate coverage report
       run: make clean cov cov_run
 
     - name: Upload coverage to Codecov
-      uses: codecov/codecov-action@v2
+      uses: codecov/codecov-action@v3
       with:
         fail_ci_if_error: false # just report, don't fail checks

.github/workflows/main_win.yml

@@ -26,12 +26,12 @@ jobs:
       shell: bash
 
     - name: Checkout repository and submodules
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
       with:
         submodules: recursive
 
     - name: Install cygwin (windows)
-      uses: cygwin/cygwin-install-action@v1
+      uses: cygwin/cygwin-install-action@v3
       with:
         packages: gcc-core gcc-g++ make
 

NEWS.md

@@ -6,9 +6,61 @@
     * Automatic conversion between matric and bulk density as needed
       using the new input `type_soilDensityInput`.
 
+* Daily weather inputs that force a simulation are now processed
+  all at once; previously, values were processed for one year at a time during
+  the main simulation loop (issue #311; @dschlaep, @N1ckP3rsl3y).
+    * Daily weather inputs are now obtained by `readAllWeather()`
+      via `SW_WTH_read()` during `SW_CTL_read_inputs_from_disk()`, i.e.,
+      the same time as other inputs are read from files.
+    * Then, weather values are "finalized", i.e., missing values are imputed
+      (e.g., by the weather generator) and scaled with monthly parameters,
+      by `finalizeAllWeather()` via `SW_WTH_finalize_all_weather()`;
+      this must occur before the simulation is "initialized"
+      by `SW_CTL_init_run()`.
+
+* SOILWAT2 gains the ability to calculate long-term climate summaries
+  (issue #317; @N1ckP3rsl3y, @dschlaep).
+    * New `calcSiteClimate()` calculates monthly and annual time
+      series of climate variables from daily weather.
+    * New `averageClimateAcrossYears()` calculates long-term climate summaries.
+    * Both functions are based on `rSOILWAT2::calc_SiteClimate()`
+      which was previously coded in R.
+    * This version fixes issues from the previous R version:
+       * Mean annual temperature is now the mean across years of
+         means across days within year of mean daily temperature.
+       * Years at locations in the southern hemisphere are now adjusted to start
+         on July 1 of the previous calendar year.
+       * Variables `Month7th_PPT_mm` and `MinTemp_of2ndMonth_C` are now adjusted
+         for location by hemisphere.
+
+* SOILWAT2 gains the ability to estimate fractional land cover
+  representing a potential natural vegetation based on climate relationships
+  (using new input `veg_method`) instead of reading land cover values
+  from input files (issue #318; @N1ckP3rsl3y, @dschlaep).
+    * New `estimatePotNatVegComposition()` estimates
+      fractional land cover representing a potential natural vegetation
+      based on climate relationships.
+      This function is based on `rSOILWAT2::estimate_PotNatVeg_composition()`
+      which was previously coded in R.
+    * New `estimateVegetationFromClimate()`
+      (which is called by `SW_VPD_init_run()`) uses `veg_method` to determine
+      at run time if a simulation utilizes `averageClimateAcrossYears()` and
+      `estimatePotNatVegComposition()` to set land cover values
+      instead of using the cover values provided in the input file.
+    * This version fixes issues from the previous R version:
+       * The `C4` grass correction based on Teeri & Stowe 1976 is now applied
+         as documented (`rSOILWAT2` issue #218).
+       * The sum of all grass components, if fixed, is now incorporated into
+         the total sum of all fixed components (`rSOILWAT2` issue #219).
+
+
 ## Changes to inputs
 * SOILWAT2 gains `type_soilDensityInput` as new user input (`siteparam.in`)
-  with default value 0 (matric soil density) that reproduces previous behavior.
+  with default value 0 (i.e., matric soil density)
+  that reproduces previous behavior.
+* SOILWAT2 gains `veg_method` as new user input (`"veg.in"`)
+  with default value 0 (i.e., land cover are obtained from input files)
+  that reproduces previous behavior.
 
 
 # SOILWAT2 v6.6.0

SW_Control.c

@@ -58,7 +58,6 @@
 static void _begin_year(void) {
 	// SW_F_new_year() not needed
 	SW_MDL_new_year(); // call first to set up time-related arrays for this year
-	SW_WTH_new_year();
 	// SW_MKV_new_year() not needed
 	SW_SKY_new_year(); // Update daily climate variables from monthly values
 	//SW_SIT_new_year() not needed
@@ -77,7 +76,6 @@ static void _begin_day(void) {
 
 static void _end_day(void) {
 	_collect_values();
-	SW_WTH_end_day();
 	SW_SWC_end_day();
 }
 
@@ -250,23 +248,28 @@ void SW_CTL_read_inputs_from_disk(void) {
   if (debug) swprintf(" > 'model'");
   #endif
 
-  SW_WTH_read();
+  SW_WTH_setup();
   #ifdef SWDEBUG
-  if (debug) swprintf(" > 'weather'");
+  if (debug) swprintf(" > 'weather setup'");
   #endif
 
   SW_SKY_read();
   #ifdef SWDEBUG
   if (debug) swprintf(" > 'climate'");
   #endif
 
-  if (SW_Weather.use_weathergenerator) {
+  if (SW_Weather.generateWeatherMethod == 2) {
     SW_MKV_setup();
     #ifdef SWDEBUG
     if (debug) swprintf(" > 'weather generator'");
     #endif
   }
 
+  SW_WTH_read();
+  #ifdef SWDEBUG
+  if (debug) swprintf(" > 'weather read'");
+  #endif
+
   SW_VPD_read();
   #ifdef SWDEBUG
   if (debug) swprintf(" > 'veg'");

SW_Flow.c

@@ -373,7 +373,7 @@ void SW_Water_Flow(void) {
 			 calculate soil temperature at end of each day
 		*/
 		SW_ST_setup_run(
-			w->now.temp_avg[Today],
+			w->now.temp_avg,
 			lyrSWCBulk,
 			lyrSWCBulk_Saturated,
 			lyrbDensity,
@@ -408,15 +408,15 @@ void SW_Water_Flow(void) {
 		x,
 		SW_Sky.cloudcov_daily[doy],
 		SW_Sky.r_humidity_daily[doy],
-		w->now.temp_avg[Today],
+		w->now.temp_avg,
 		&sw->H_oh,
 		&sw->H_ot,
 		&sw->H_gh
 	);
 
 	sw->pet = SW_Site.pet_scale * petfunc(
 		sw->H_gt,
-		w->now.temp_avg[Today],
+		w->now.temp_avg,
 		SW_Site.altitude,
 		x,
 		SW_Sky.r_humidity_daily[doy],
@@ -444,7 +444,7 @@ void SW_Water_Flow(void) {
 	}
 
 	/* Rainfall interception */
-	h2o_for_soil = w->now.rain[Today]; /* ppt is partioned into ppt = snow + rain */
+	h2o_for_soil = w->now.rain; /* ppt is partioned into ppt = snow + rain */
 
   ForEachVegType(k)
   {
@@ -837,13 +837,13 @@ void SW_Water_Flow(void) {
 	// soil_temperature function computes the soil temp for each layer and stores it in lyravgLyrTemp
 	// doesn't affect SWC at all (yet), but needs it for the calculation, so therefore the temperature is the last calculation done
 	if (SW_Site.use_soil_temp) {
-		soil_temperature(w->now.temp_avg[Today], sw->pet, sw->aet, x, lyrSWCBulk,
+		soil_temperature(w->now.temp_avg, sw->pet, sw->aet, x, lyrSWCBulk,
 			lyrSWCBulk_Saturated, lyrbDensity, lyrWidths, lyroldavgLyrTemp, lyravgLyrTemp, surfaceAvg,
 			SW_Site.n_layers, SW_Site.bmLimiter,
 			SW_Site.t1Param1, SW_Site.t1Param2, SW_Site.t1Param3, SW_Site.csParam1,
 			SW_Site.csParam2, SW_Site.shParam, sw->snowdepth, SW_Site.Tsoil_constant,
 			SW_Site.stDeltaX, SW_Site.stMaxDepth, SW_Site.stNRGR, sw->snowpack[Today],
-			&SW_Soilwat.soiltempError, w->now.temp_max[Today], w->now.temp_min[Today],
+			&SW_Soilwat.soiltempError, w->now.temp_max, w->now.temp_min,
             sw->H_gt, sw->maxLyrTemperature, sw->minLyrTemperature, &w->surfaceMax, &w->surfaceMin);
 	}
 

SW_Main.c

@@ -83,6 +83,9 @@ int main(int argc, char **argv) {
 	// read user inputs
 	SW_CTL_read_inputs_from_disk();
 
+	// finalize daily weather
+	SW_WTH_finalize_all_weather();
+
 	// initialize simulation run (based on user inputs)
 	SW_CTL_init_run();
 

SW_Markov.c

@@ -297,6 +297,15 @@ void SW_MKV_today(TimeInt doy0, RealD *tmax, RealD *tmin, RealD *rain) {
 	short debug = 0;
 	#endif
 
+	#ifdef SWDEBUG
+	if (debug) {
+		swprintf(
+			"mkv(before): yr=%u/doy0=%u: ppt=%.3f, tmax=%.3f, tmin=%.3f\n",
+			SW_Model.year, doy0, *rain, *tmax, *tmin
+		);
+	}
+	#endif
+
 	/* Calculate Precipitation:
 		prop = probability that it precipitates today depending on whether it
 			was wet (precipitated) yesterday `wetprob` or
@@ -335,8 +344,10 @@ void SW_MKV_today(TimeInt doy0, RealD *tmax, RealD *tmin, RealD *rain) {
 
 	#ifdef SWDEBUG
 	if (debug) {
-		swprintf("mkv: yr=%d/doy0=%d/week=%d: ppt=%.3f, tmax=%.3f, tmin=%.3f\n",
-			SW_Model.year, doy0, week, *rain, *tmax, *tmin);
+		swprintf(
+			"mkv(after): yr=%u/doy0=%u/week=%u: ppt=%.3f, tmax=%.3f, tmin=%.3f\n",
+			SW_Model.year, doy0, week, *rain, *tmax, *tmin
+		);
 	}
 	#endif
 
@@ -525,14 +536,22 @@ Bool SW_MKV_read_cov(void) {
 void SW_MKV_setup(void) {
   SW_MKV_construct();
 
-  if (!SW_MKV_read_prob()) {
-    LogError(logfp, LOGFATAL, "Markov weather requested but could not open %s",
-      SW_F_name(eMarkovProb));
+  if (!SW_MKV_read_prob() && SW_Weather.generateWeatherMethod == 2) {
+    LogError(
+      logfp,
+      LOGFATAL,
+      "Weather generator requested but could not open %s",
+      SW_F_name(eMarkovProb)
+    );
   }
 
-  if (!SW_MKV_read_cov()) {
-    LogError(logfp, LOGFATAL, "Markov weather requested but could not open %s",
-      SW_F_name(eMarkovCov));
+  if (!SW_MKV_read_cov() && SW_Weather.generateWeatherMethod == 2) {
+    LogError(
+      logfp,
+      LOGFATAL,
+      "Weather generator requested but could not open %s",
+      SW_F_name(eMarkovCov)
+    );
   }
 }
 

SW_Model.c

@@ -43,7 +43,6 @@
 
 #include "SW_Defines.h"
 #include "SW_Files.h"
-#include "SW_Weather.h"
 #include "SW_Site.h" // externs SW_Site
 #include "SW_SoilWater.h"  /* for setup_new_year() */
 #include "SW_Times.h"

SW_Output.c

@@ -354,17 +354,17 @@ static void sumof_wth(SW_WEATHER *v, SW_WEATHER_OUTPUTS *s, OutKey k)
 	{
 
 	case eSW_Temp:
-		s->temp_max += v->now.temp_max[Today];
-		s->temp_min += v->now.temp_min[Today];
-		s->temp_avg += v->now.temp_avg[Today];
+		s->temp_max += v->now.temp_max;
+		s->temp_min += v->now.temp_min;
+		s->temp_avg += v->now.temp_avg;
 		//added surfaceAvg for sum
         s->surfaceAvg += v->surfaceAvg;
         s->surfaceMax += v->surfaceMax;
         s->surfaceMin += v->surfaceMin;
 		break;
 	case eSW_Precip:
-		s->ppt += v->now.ppt[Today];
-		s->rain += v->now.rain[Today];
+		s->ppt += v->now.ppt;
+		s->rain += v->now.rain;
 		s->snow += v->snow;
 		s->snowmelt += v->snowmelt;
 		s->snowloss += v->snowloss;

SW_Output.h

@@ -59,7 +59,6 @@
 #include "Times.h"
 #include "SW_Defines.h"
 #include "SW_SoilWater.h"
-#include "SW_Weather.h"
 #include "SW_VegProd.h"
 
 #ifdef __cplusplus

SW_SoilWater.c

@@ -194,7 +194,7 @@ void SW_WaterBalance_Checks(void)
   runoff = w->snowRunoff + w->surfaceRunoff;
   runon = w->surfaceRunon;
   snowmelt = w->snowmelt;
-  rain = w->now.rain[Today];
+  rain = w->now.rain;
 
   arriving_water = rain + snowmelt + runon;
 

SW_VegEstab.c

@@ -231,11 +231,11 @@ void SW_VES_checkestab(void) {
 static void _checkit(TimeInt doy, unsigned int sppnum) {
 
 	SW_VEGESTAB_INFO *v = SW_VegEstab.parms[sppnum];
-	SW_WEATHER_2DAYS *wn = &SW_Weather.now;
+	SW_WEATHER_NOW *wn = &SW_Weather.now;
 	SW_SOILWAT *sw = &SW_Soilwat;
 
 	IntU i;
-	RealF avgtemp = wn->temp_avg[Today], /* avg of today's min/max temp */
+	RealF avgtemp = wn->temp_avg, /* avg of today's min/max temp */
 	avgswc; /* avg_swc today */
 
 	if (doy == SW_Model.firstdoy) {