Bugfix435: calculated relative humidity

NEWS.md

@@ -1,5 +1,15 @@
 # NEWS
 
+# SOILWAT2 v8.0.1
+* Fix the calculation of relative humidity (#435; @dschlaep).
+  Previously, relative humidity was incorrectly calculated if based on
+  vapor pressure or specific humidity.
+* Fix the count of days on which a missing weather value was replaced by a
+  non-missing value from the preceding day for the method `"LOCF"`
+  (last observation carried forward; #437; @dschlaep). Previously, any day
+  with a missing weather value was counted.
+
+
 # SOILWAT2 v8.0.0
 * Simulation output remains the same as the previous version.
   However, output of establishment/recruitment for two species is now

include/SW_Flow_lib_PET.h

@@ -96,15 +96,17 @@ double petfunc(
     LOG_INFO *LogInfo
 );
 
+double atmospheric_pressure(double elev);
+
+double psychrometric_constant(double pressure);
 
 double svp(double T, double *slope_svp_to_t);
 
 double svp2(double temp);
 
-double atmospheric_pressure(double elev);
-
-double psychrometric_constant(double pressure);
+double relativeHumidity1(double vp, double meanTemp);
 
+double relativeHumidity2(double huss, double meanTemp, double elevation);
 
 double actualVaporPressure1(double hurs, double meanTemp);
 
@@ -114,6 +116,7 @@ double actualVaporPressure2(
 
 double actualVaporPressure3(double dewpointTemp);
 
+double actualVaporPressure4(double huss, double elevation);
 
 #ifdef __cplusplus
 }

include/SW_Weather.h

@@ -150,6 +150,7 @@ void _read_weather_hist(
     unsigned int n_input_forcings,
     unsigned int *dailyInputIndices,
     Bool *dailyInputFlags,
+    double elevation,
     LOG_INFO *LogInfo
 );
 
@@ -168,6 +169,7 @@ void readAllWeather(
     RealD *cloudcov,
     RealD *windspeed,
     RealD *r_humidity,
+    double elevation,
     TimeInt cum_monthdays[],
     TimeInt days_in_month[],
     LOG_INFO *LogInfo

include/SW_datastructs.h

@@ -617,19 +617,58 @@ typedef struct {
 /*                   Weather structs                   */
 /* --------------------------------------------------- */
 
+/** Weather values of the current simulation day */
 typedef struct {
-    /* Weather values of the current simulation day */
-    RealD temp_avg, temp_max, temp_min, ppt, rain, cloudCover, windSpeed,
-        relHumidity, shortWaveRad, actualVaporPressure;
+    /** Daily near-surface average air temperature [C] */
+    double temp_avg;
+    /** Daily near-surface maximum air temperature [C] */
+    double temp_max;
+    /** Daily near-surface minimum air temperature [C] */
+    double temp_min;
+    /** Daily precipitation amount [cm] */
+    double ppt;
+    /** Daily precipitation amount that falls as rain [cm] */
+    double rain;
+    /** Daily cloud cover [%] */
+    double cloudCover;
+    /** Daily mean near-surface wind speed [m s-1] */
+    double windSpeed;
+    /** Daily mean near-surface relative humidity [%] */
+    double relHumidity;
+    /** Daily downward surface shortwave radiation:
+        global horizontal irradiation [MJ / m2] or flux density [W / m2],
+        see #SW_WEATHER.desc_rsds
+    */
+    double shortWaveRad;
+    /** Daily mean near-surface actual vapor pressure [kPa] */
+    double actualVaporPressure;
 } SW_WEATHER_NOW;
 
+/** Daily weather values for one calendar year */
 typedef struct {
-    /* Daily weather values for one year */
-    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
+    /** Daily maximum near-surface air temperature [C] */
+    double temp_max[MAX_DAYS];
+    /** Daily minimum near-surface air temperature [C] */
+    double temp_min[MAX_DAYS];
+    /** Daily average near-surface air temperature [C] */
+    double temp_avg[MAX_DAYS];
+    /** Daily precipitation amount [cm] */
+    double ppt[MAX_DAYS];
+    /** Daily cloud cover [%] */
+    double cloudcov_daily[MAX_DAYS];
+    /** Daily mean near-surface wind speed [m s-1] */
+    double windspeed_daily[MAX_DAYS];
+    /** Daily mean near-surface relative humidity [%] */
+    double r_humidity_daily[MAX_DAYS];
+    /** Daily downward surface shortwave radiation:
+        global horizontal irradiation [MJ / m2] or flux density [W / m2],
+        see #SW_WEATHER.desc_rsds
+    */
+    double shortWaveRad[MAX_DAYS];
+    /** Daily mean near-surface actual vapor pressure [kPa] */
+    double actualVaporPressure[MAX_DAYS];
+
+    // double temp_month_avg[MAX_MONTHS], temp_year_avg; // currently not used
 } SW_WEATHER_HIST;
 
 /* accumulators for output values hold only the */
@@ -925,18 +964,23 @@ typedef struct {
 /*                    Sky structs                      */
 /* --------------------------------------------------- */
 
+/** Across-year mean monthly climate variables */
 typedef struct {
-    RealD cloudcov[MAX_MONTHS],     /* monthly cloud cover (frac) */
-        windspeed[MAX_MONTHS],      /* windspeed (m/s) */
-        r_humidity[MAX_MONTHS],     /* relative humidity (%) */
-        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 snow_density_daily[MAX_DAYS + 1]; /* interpolated daily snow density
-                                               (kg/m3) */
-
+    /** Across-year mean monthly mean cloud cover [%] */
+    double cloudcov[MAX_MONTHS];
+    /** Across-year mean mean monthly mean near-surface wind speed [m s-1] */
+    double windspeed[MAX_MONTHS];
+    /** Across-year mean monthly mean near-surface relative humidity [%] */
+    double r_humidity[MAX_MONTHS];
+    /** Across-year mean monthly mean snow density [kg m-3] */
+    double snow_density[MAX_MONTHS];
+    /** Across-year mean monthly count of precipitation events */
+    /* currently used in interception functions */
+    double n_rain_per_day[MAX_MONTHS];
+
+    /** Across-year mean daily snow density [kg m-3],
+        interpolated from monthly values */
+    double snow_density_daily[MAX_DAYS + 1];
 } SW_SKY;
 
 /* =================================================== */

src/SW_Flow_lib_PET.c

@@ -1242,6 +1242,35 @@ double svp2(double temp) {
     return .6108 * exp((17.27 * temp) / (temp + 237.3));
 }
 
+/**
+@brief Calculate relative humidity from vapor pressure and temperature
+
+@param vp Daily mean vapor pressure [kPa]
+@param meanTemp Daily mean air temperature [C]
+
+@return Calculated relative humidity [0-100 %]
+*/
+double relativeHumidity1(double vp, double meanTemp) {
+    double tempSlope;
+    double svpVal = svp(meanTemp, &tempSlope);
+    return 100. * vp / svpVal;
+}
+
+/**
+@brief Calculate relative humidity from specific humidity and temperature
+
+@param huss Daily mean specific humidity [g kg-1]
+@param meanTemp Daily mean air temperature [C]
+@param elevation Site elevation [m above mean sea level]
+
+@return Calculated relative humidity [0-100 %]
+*/
+double relativeHumidity2(double huss, double meanTemp, double elevation) {
+    double vpVal = actualVaporPressure4(huss, elevation);
+
+    return relativeHumidity1(vpVal, meanTemp);
+}
+
 /**
 @brief Calculate actual vapor pressure based on relative humidity and mean
 temperature
@@ -1294,6 +1323,34 @@ double actualVaporPressure3(double dewpointTemp) {
     return svp2(dewpointTemp);
 }
 
+/**
+@brief Calculate actual vapor pressure from specific humidity and elevation
+
+@param huss Daily mean specific humidity [g kg-1]
+@param elevation Site elevation [m above mean sea level].
+
+@return Calculated actual vapor pressure [kPa]
+*/
+double actualVaporPressure4(double huss, double elevation) {
+    // Rd = specific gas constant for dry air (J kg−1 K−1)
+    // Rv = specific gas constant for water vapor (J kg−1 K−1)
+    // e = vapor pressure (Pa)
+    // p = air pressure (Pa)
+
+    double vpVal;
+
+    /* Rd / Rv = 287.052874 [J kg−1 K−1] / 461.50 [J kg−1 K−1] */
+    static const double RdToRv = 0.6220;
+    double p = atmospheric_pressure(elevation);
+
+    huss *= 1e-3; /* convert [g kg-1] to [kg kg-1] */
+
+    /* huss = q = Rd / Rv * e / (p - e * (1 - Rd / Rv)) */
+    vpVal = huss * p / (RdToRv + huss * (1. - RdToRv));
+
+    return vpVal;
+}
+
 /**
 @brief Daily potential evapotranspiration