Merge pull request #338 from DrylandEcology/read_more_daily_input

Expanded daily input variables

Daily weather inputs, in addition to the previous variables maximum air temperature, minimum air temperature, and precipitation amount, can now process the following variables (close #341):
    * Cloud cover (can be replaced by shortwave radiation)
    * Wind speed (can be replaced by wind components)
    * Wind speed eastward component (optional)
    * Wind speed northward component (optional)
    * Relative humidity (can be replaced by max/min humidity, specific humidity, dew point temperature, or vapor pressure)
    * Maximum relative humidity (optional)
    * Minimum relative humidity (optional)
    * Specific humidity (optional)
    * Dew point temperature (optional)
    * Actual vapor pressure (optional)
    * Downward surface shortwave radiation (optional)

NEWS.md

@@ -32,6 +32,26 @@
       by `finalizeAllWeather()` via `SW_WTH_finalize_all_weather()`;
       this must occur before the simulation is "initialized"
       by `SW_CTL_init_run()`.
+    * Imputation of missing daily weather values by `generateMissingWeather()`
+      can be done either by last-observation-carried-forward `"LOCF"`
+      (which handles all daily weather variables) or by the Markov weather
+      generator (which handles only temperature and precipitation).
+
+* Daily weather inputs, in addition to the previous variables
+  maximum air temperature, minimum air temperature, and precipitation amount,
+  can now process the following variables (issue #341; @dschlaep, @N1ckP3rsl3y):
+    * Cloud cover (can be replaced by shortwave radiation)
+    * Wind speed (can be replaced by wind components)
+    * Wind speed eastward component (optional)
+    * Wind speed northward component (optional)
+    * Relative humidity (can be replaced by max/min humidity, specific humidity
+      dew point temperature, or vapor pressure)
+    * Maximum relative humidity (optional)
+    * Minimum relative humidity (optional)
+    * Specific humidity (optional)
+    * Dew point temperature (optional)
+    * Actual vapor pressure (optional)
+    * Downward surface shortwave radiation (optional)
 
 * SOILWAT2 gains the ability to calculate long-term climate summaries
   (issue #317; @N1ckP3rsl3y, @dschlaep).
@@ -71,6 +91,12 @@
 
 
 ## Changes to inputs
+* New inputs via `"weathsetup.in"` determine whether monthly or daily inputs
+  for cloud cover, relative humidity, and wind speed are utilized;
+  describe which daily weather variables are contained in the weather input
+  files `weath.YYYY`; and describe units of (optiona) input shortwave radiation.
+* New (optional) variables (columns) in weather input files `weath.YYYY` that
+  are described via `"weathsetup.in"`.
 * 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

include/SW_Defines.h

@@ -59,7 +59,6 @@ extern "C" {
 
 #define SW_MISSING     999.     /**< Value to use as MISSING */
 
-
 // Euler's constant
 #ifdef M_E
   #define swE M_E
@@ -117,6 +116,26 @@ extern "C" {
 #define SW_FORBS 2
 #define SW_GRASS 3
 
+/*
+   Indices to daily input flags/indices (dailyInputFlags & dailyInputIndices in SW_WEATHER)
+   The order of these indices must match the order of weather input flags within `weathsetup.in`
+*/
+#define MAX_INPUT_COLUMNS 14    /**< Maximum number of columns that can be input in a weath.YYYY file*/
+#define TEMP_MAX      0
+#define TEMP_MIN      1
+#define PPT           2
+#define CLOUD_COV     3
+#define WIND_SPEED    4
+#define WIND_EAST     5
+#define WIND_NORTH    6
+#define REL_HUMID     7
+#define REL_HUMID_MAX 8
+#define REL_HUMID_MIN 9
+#define SPEC_HUMID    10
+#define TEMP_DEWPOINT 11
+#define ACTUAL_VP     12
+#define SHORT_WR      13
+
 
 /* output period specifiers */
 #define SW_DAY   "DY"

include/SW_Flow_lib_PET.h

@@ -35,11 +35,15 @@ double overcast_attenuation_KastenCzeplak1980(double cloud_cover);
 double overcast_attenuation_Angstrom1924(double sunshine_fraction);
 double clearsky_directbeam(double P, double e_a, double int_sin_beta);
 double clearnessindex_diffuse(double K_b);
+double actual_horizontal_transmissivityindex(double tau);
 
-double solar_radiation(unsigned int doy,
+double solar_radiation(
+  unsigned int doy,
   double lat, double elev, double slope, double aspect,
-  double albedo, double cloud_cover, double rel_humidity, double air_temp_mean,
-  double *H_oh, double *H_ot, double *H_gh);
+  double albedo, double *cloud_cover, double e_a,
+  double rsds, unsigned int desc_rsds,
+  double *H_oh, double *H_ot, double *H_gh
+);
 
 
 double blackbody_radiation(double T);
@@ -49,9 +53,14 @@ double petfunc(double H_g, double avgtemp, double elev,
   double reflec, double humid, double windsp, double cloudcov);
 
 double svp(double T, double *slope_svp_to_t);
+double svp2(double temp);
 double atmospheric_pressure(double elev);
 double psychrometric_constant(double pressure);
 
+double actualVaporPressure1(double hurs, double meanTemp);
+double actualVaporPressure2(double maxHurs, double minHurs, double maxTemp, double minTemp);
+double actualVaporPressure3(double dewpointTemp);
+
 #ifdef __cplusplus
 }
 #endif

include/SW_Sky.h

@@ -32,10 +32,7 @@ typedef struct {
           snow_density [MAX_MONTHS], /* snow density (kg/m3) */
           n_rain_per_day[MAX_MONTHS]; /* number of precipitation events per month (currently used in interception functions) */
 
-    RealD cloudcov_daily     [MAX_DAYS+1], /* interpolated daily cloud cover (frac) */
-          windspeed_daily    [MAX_DAYS+1], /* interpolated daily windspeed (m/s) */
-          r_humidity_daily   [MAX_DAYS+1], /* interpolated daily relative humidity (%) */
-          snow_density_daily	[MAX_DAYS+1];	/* interpolated daily snow density (kg/m3) */
+    RealD snow_density_daily	[MAX_DAYS+1];	/* interpolated daily snow density (kg/m3) */
 
 } SW_SKY;
 
@@ -51,7 +48,6 @@ extern SW_SKY SW_Sky;
 /*             Global Function Declarations            */
 /* --------------------------------------------------- */
 void SW_SKY_read(void);
-void SW_SKY_init_run(void);
 void SW_SKY_new_year(void);
 
 #ifdef __cplusplus

include/SW_Weather.h

@@ -41,12 +41,15 @@ extern "C" {
 
 typedef struct {
 	/* Weather values of the current simulation day */
-	RealD temp_avg, temp_max, temp_min, ppt, rain;
+	RealD temp_avg, temp_max, temp_min, ppt, rain, cloudCover, windSpeed, relHumidity,
+    shortWaveRad, actualVaporPressure;
 } SW_WEATHER_NOW;
 
 typedef struct {
 	/* Daily weather values for one year */
-	RealD temp_max[MAX_DAYS], temp_min[MAX_DAYS], temp_avg[MAX_DAYS], ppt[MAX_DAYS];
+	RealD temp_max[MAX_DAYS], temp_min[MAX_DAYS], temp_avg[MAX_DAYS], ppt[MAX_DAYS],
+    cloudcov_daily[MAX_DAYS], windspeed_daily[MAX_DAYS], r_humidity_daily[MAX_DAYS],
+    shortWaveRad[MAX_DAYS], actualVaporPressure[MAX_DAYS];
 	// RealD temp_month_avg[MAX_MONTHS], temp_year_avg; // currently not used
 } SW_WEATHER_HIST;
 
@@ -153,11 +156,20 @@ typedef struct {
     scale_temp_min[MAX_MONTHS],
     scale_skyCover[MAX_MONTHS],
     scale_wind[MAX_MONTHS],
-    scale_rH[MAX_MONTHS];
+    scale_rH[MAX_MONTHS],
+    scale_actVapPress[MAX_MONTHS],
+    scale_shortWaveRad[MAX_MONTHS];
 	char name_prefix[MAX_FILENAMESIZE - 5]; // subtract 4-digit 'year' file type extension
 	RealD snowRunoff, surfaceRunoff, surfaceRunon, soil_inf, surfaceAvg;
 	RealD snow, snowmelt, snowloss, surfaceMax, surfaceMin;
 
+  Bool use_cloudCoverMonthly, use_windSpeedMonthly, use_humidityMonthly;
+  Bool dailyInputFlags[MAX_INPUT_COLUMNS];
+
+  unsigned int dailyInputIndices[MAX_INPUT_COLUMNS],
+ 			   n_input_forcings, // Number of input columns found in weath.YYYY
+    desc_rsds; /**< Description of units and definition of daily inputs of observed shortwave radiation, see `solar_radiation()` */
+
 	/* This section is required for computing the output quantities.  */
 	SW_WEATHER_OUTPUTS
 		*p_accu[SW_OUTNPERIODS], // output accumulator: summed values for each time period
@@ -185,6 +197,17 @@ extern SW_WEATHER SW_Weather;
 /*             Global Function Declarations            */
 /* --------------------------------------------------- */
 void SW_WTH_setup(void);
+void set_dailyInputIndices(
+  Bool dailyInputFlags[MAX_INPUT_COLUMNS],
+  unsigned int dailyInputIndices[MAX_INPUT_COLUMNS],
+  unsigned int *n_input_forcings
+);
+void check_and_update_dailyInputFlags(
+  Bool use_cloudCoverMonthly,
+  Bool use_humidityMonthly,
+  Bool use_windSpeedMonthly,
+  Bool *dailyInputFlags
+);
 void SW_WTH_read(void);
 void averageClimateAcrossYears(SW_CLIMATE_YEARLY *climateOutput, int numYears,
                                SW_CLIMATE_CLIM *climateAverages);
@@ -204,14 +227,26 @@ void deallocateClimateStructs(SW_CLIMATE_YEARLY *climateOutput,
 void _read_weather_hist(
   TimeInt year,
   SW_WEATHER_HIST *yearWeather,
-  char weather_prefix[]
+  char weather_prefix[],
+  unsigned int n_input_forcings,
+  unsigned int *dailyInputIndices,
+  Bool *dailyInputFlags
 );
 void readAllWeather(
   SW_WEATHER_HIST **allHist,
   int startYear,
   unsigned int n_years,
   Bool use_weathergenerator_only,
-  char weather_prefix[]
+  char weather_prefix[],
+  Bool use_cloudCoverMonthly,
+  Bool use_humidityMonthly,
+  Bool use_windSpeedMonthly,
+  unsigned int n_input_forcings,
+  unsigned int *dailyInputIndices,
+  Bool *dailyInputFlags,
+  RealD *cloudcov,
+  RealD *windspeed,
+  RealD *r_humidity
 );
 void finalizeAllWeather(SW_WEATHER *w);
 
@@ -221,7 +256,12 @@ void scaleAllWeather(
   unsigned int n_years,
   double *scale_temp_max,
   double *scale_temp_min,
-  double *scale_precip
+  double *scale_precip,
+  double *scale_skyCover,
+  double *scale_wind,
+  double *scale_rH,
+  double *scale_actVapPress,
+  double *scale_shortWaveRad
 );
 void generateMissingWeather(
   SW_WEATHER_HIST **allHist,
@@ -230,6 +270,7 @@ void generateMissingWeather(
   unsigned int method,
   unsigned int optLOCF_nMax
 );
+void checkAllWeather(SW_WEATHER *weather);
 void allocateAllWeather(SW_WEATHER *w);
 void deallocateAllWeather(SW_WEATHER *w);
 void _clear_hist_weather(SW_WEATHER_HIST *yearWeather);

include/Times.h

@@ -100,7 +100,8 @@ TimeInt yearto4digit(TimeInt yr);
 
 Bool isleapyear(const TimeInt year);
 
-void interpolate_monthlyValues(double monthlyValues[], double dailyValues[]);
+void interpolate_monthlyValues(double monthlyValues[], Bool interpAsBase1,
+                               double dailyValues[]);
 
 
 #ifdef __cplusplus

src/SW_Control.c

@@ -156,7 +156,7 @@ void SW_CTL_init_run(void) {
 	SW_WTH_init_run();
 	// SW_MKV_init_run() not needed
 	SW_PET_init_run();
-	SW_SKY_init_run();
+	// SW_SKY_init_run() not needed
 	SW_SIT_init_run();
 	SW_VES_init_run(); // must run after `SW_SIT_init_run()`
 	SW_VPD_init_run();

src/SW_Flow.c

@@ -406,9 +406,10 @@ void SW_Water_Flow(void) {
 		SW_Site.slope,
 		SW_Site.aspect,
 		x,
-		SW_Sky.cloudcov_daily[doy],
-		SW_Sky.r_humidity_daily[doy],
-		w->now.temp_avg,
+        &w->now.cloudCover,
+        w->now.actualVaporPressure,
+        w->now.shortWaveRad,
+        w->desc_rsds,
 		&sw->H_oh,
 		&sw->H_ot,
 		&sw->H_gh
@@ -419,9 +420,9 @@ void SW_Water_Flow(void) {
 		w->now.temp_avg,
 		SW_Site.altitude,
 		x,
-		SW_Sky.r_humidity_daily[doy],
-		SW_Sky.windspeed_daily[doy],
-		SW_Sky.cloudcov_daily[doy]
+        w->now.relHumidity,
+        w->now.windSpeed,
+        w->now.cloudCover
 	);