Merge pull request #322 from UW-Hydro/hotfix/4.2.c

4.2.c hotfix release

docs/Documentation/VICDisagg.md

@@ -15,13 +15,14 @@ Create a [global parameter file](GlobalParam.md). Because VIC will not actually
 *   `STARTYEAR`,` STARTMONTH`, `STARTDAY`, `STARTHOUR`: Set these to the start date for your disaggregated forcings
 *   `ENDYEAR`, `ENDMONTH`, `ENDDAY`: Set these to the end date for your disaggregated forcings
 *   All the variables in the forcing section should be the same as before, i.e. these describe the **input** (daily) forcings that you are **reading**: `FORCING1`, `FORCING2` (if applicable), `FORCE_FORMAT`, `FORCE_ENDIAN`, `N_TYPES`, `FORCE_TYPE` (there must be one of these for each input variable, e.g. PREC, TMAX, TMIN, WIND), `FORCE_DT`, `FORCEYEAR`, `FORCEMONTH`, `FORCEDAY`, `FORCEHOUR`, `GRID_DECIMAL`, `WIND_H`, `MEASURE_H`, and `ALMA_INPUT`
+*   `VEGLIB`,`VEGLIB_VEGCOVER`,`VEGPARAM`,`ROOT_ZONES`,`VEGPARAM_LAI`,`VEGPARAM_VEGCOVER`,`VEGPARAM_ALB`: Set to the same values as you would for a full VIC simulation
 *   `RESULT_DIR`: Set to the name of the directory where the disaggregated forcings should be written
 *   `OUT_STEP`: Set to 0
 *   `ALMA_OUTPUT`: For standard VIC forcings, set to FALSE; for ALMA-compliant forcings (often required by other models) set to TRUE
 *   `BINARY_OUTPUT`: Set this to FALSE to produce ASCII forcings, TRUE to produce BINARY forcings
 *   `SKIPYEAR`: We recommend setting this to 0
 *   `OUTPUT_FORCE`: This must be set to TRUE
-*   `N_OUTFILES`, `OUTFILE`, `OUTVAR`: These can be omitted; by default, VIC will produce 1 output file per grid cell, names "full_data__lat___lon_," where lat, lon = latitude and longitude of te grid cell's center. These default output files will contain the following variables:
+*   `N_OUTFILES`, `OUTFILE`, `OUTVAR`: These can be omitted; by default, VIC will produce 1 output file per grid cell, names "full_data_lat_lon_," where lat, lon = latitude and longitude of te grid cell's center. These default output files will contain the following variables:
 
 | Name          | Units (ALMA_OUTPUT FALSE)     | Units (ALMA_OUTPUT TRUE)  |
 |-----------    |---------------------------    |-------------------------- |

src/ChangeLog

@@ -24,6 +24,17 @@ Usage:
 
 
 
+ -------------------------------------------------------------------------------
+ ***** Description of changes between VIC 4.2.b and VIC 4.2.c *****
+ -------------------------------------------------------------------------------
+
+ Bug Fixes:
+ ----------
+
+Fixed water balance errors in lake model (Github Issues: #308, #316)
+Fixed OUTPUT_FORCE behavior related to veg_hist features (Github Issues: #305, #307)
+Updated bare soil evaporation behavior related to veg_hist update from 4.2 (Github Issues: #306, #315)
+
 
 -------------------------------------------------------------------------------
 ***** Description of changes between VIC 4.2.a and VIC 4.2.b *****
@@ -781,7 +792,7 @@ Incorrect handling of user-supplied tskc (cloud fraction) for LW_CLOUD==LW_CLOUD
 	Previous versions of VIC (before 4.1.2) used a full-sky longwave
 	formulation taken from two formulas in the Bras hydrology text.  For
 	the new Deardorff full-sky longwave formulation, the dependence on
-	cloud fraction is different from the old Bras formulation.   In 
+	cloud fraction is different from the old Bras formulation.   In
 	4.1.2 (and 4.1.2.a-b), the new Deardorff formulation did not account
 	for the possibility of user-supplied cloud fraction; if the user
 	supplied cloud fraction as an input forcing, the resulting longwave
@@ -6235,4 +6246,3 @@ soils in the upper Mississippi River basin, J. Geophys. Res., 104(D16),
 Thornton, P.E. and S.W. Running, An improved algorithm for estimating
 incident daily solar radiation from measurements of temperature, humidity,
 and precipitation, Ag. For. Met., 93, 211-228, 1999.
-

src/arno_evap.c

@@ -97,7 +97,7 @@ double arno_evap(layer_data_struct *layer,
 
   /* Calculate the potential bare soil evaporation (mm/time step) */
 
-  Epot = penman(air_temp, elevation, rad, vpd, ra, 0.0, 0.0) * delta_t / SEC_PER_DAY;
+  Epot = penman(air_temp, elevation, rad, vpd, ra, 0.0, RARC_SOIL) * delta_t / SEC_PER_DAY;
 
   /**********************************************************************/
   /*  Compute temporary infiltration rate based on given soil_moist.    */

src/calc_veg_params.c

@@ -15,7 +15,7 @@ double calc_veg_displacement(double height) {
 
   double value;
 
-  value = 0.67 * height;
+  value = RATIO_DH_HEIGHT_VEG * height;
 
   return (value);
 
@@ -31,7 +31,7 @@ double calc_veg_height(double displacement) {
 
   double value;
 
-  value = displacement / 0.67;
+  value = displacement / RATIO_DH_HEIGHT_VEG;
 
   return (value);
 
@@ -48,7 +48,7 @@ double calc_veg_roughness(double height) {
 
   double value;
 
-  value = 0.123 * height;
+  value = RATIO_RL_HEIGHT_VEG * height;
 
   return (value);
 

src/check_files.c

@@ -23,9 +23,9 @@ void check_files(filep_struct     *filep,
   extern FILE          *open_file(char string[], char type[]);
 
   filep->soilparam   = open_file(fnames->soil, "r");
+  filep->veglib      = open_file(fnames->veglib, "r");
+  filep->vegparam    = open_file(fnames->veg, "r");
   if (!options.OUTPUT_FORCE) {
-    filep->veglib      = open_file(fnames->veglib, "r");
-    filep->vegparam    = open_file(fnames->veg, "r");
     if(options.SNOW_BAND>1)
       filep->snowband    = open_file(fnames->snowband, "r");
     if ( options.LAKES )

src/full_energy.c

@@ -330,8 +330,8 @@ int  full_energy(int                  gridcell,
 
         /* Initialize wind speeds */
         tmp_wind[0] = atmos->wind[NR];
-        tmp_wind[1] = -999.;
-        tmp_wind[2] = -999.;
+        tmp_wind[1] = ERROR;
+        tmp_wind[2] = ERROR;
 
         /* Set surface descriptive variables */
         if (p < N_PET_TYPES_NON_NAT) {
@@ -388,7 +388,7 @@ int  full_energy(int                  gridcell,
           if (dmy->hour == 0) {
             calc_Nscale_factors(veg_lib[veg_class].NscaleFlag,
                                 veg_con[iveg].CanopLayerBnd,
-                                veg_lib[veg_class].LAI[dmy[rec].month-1],
+                                veg_con[iveg].LAI[dmy[rec].month-1],
                                 soil_con->lat,
                                 soil_con->lng,
                                 soil_con->time_zone_lng,
@@ -543,9 +543,11 @@ int  full_energy(int                  gridcell,
     oldsnow = lake_var->snow.swq;
     snowprec = gauge_correction[SNOW] * (atmos->prec[NR] - rainonly);
     rainprec = gauge_correction[SNOW] * rainonly;
-    atmos->out_prec += (snowprec + rainprec) * lake_con->Cl[0] * lakefrac;
-    atmos->out_rain += rainprec * lake_con->Cl[0] * lakefrac;
-    atmos->out_snow += snowprec * lake_con->Cl[0] * lakefrac;
+    Cv = veg_con[iveg].Cv;
+    Cv *= lakefrac;
+    atmos->out_prec += (snowprec + rainprec) * Cv;
+    atmos->out_rain += rainprec * Cv;
+    atmos->out_snow += snowprec * Cv;
 
     ErrorFlag = solve_lake(snowprec, rainprec, atmos->air_temp[NR],
                            atmos->wind[NR], atmos->vp[NR] / 1000.,

src/func_surf_energy_bal.c

@@ -295,11 +295,15 @@ double func_surf_energy_bal(double Ts, va_list ap)
   double D1_plus;
   double *transp;
   double Ra_bare[3];
+  double Ra_veg[3];
   double tmp_wind[3];
   double tmp_height;
   double tmp_displacement[3];
   double tmp_roughness[3];
   double tmp_ref_height[3];
+  double ga_veg;
+  double ga_bare;
+  double ga_average;
 
   /************************************
     Read variables from variable list 
@@ -687,15 +691,71 @@ double func_surf_energy_bal(double Ts, va_list ap)
 			    TMean, Tair, wind[UnderStory], roughness[UnderStory]);
   else
     Ra_used[0] = HUGE_RESIST;
-
+
+  /*************************************************
+    Compute aerodynamic resistance for the case of exposed soil between
+    plants (or gaps in canopy).  Assume plants and exposed soil are well-
+    mixed, i.e., exposed soil is neither as disconnected from the
+    atmosphere as soil under veg or canopy, nor as exposed as a large
+    area of exposed soil.  Rather, it is subject to some wind attenuation
+    from the surrounding plants.  Thus, compute as the area-weighted
+    average of "pure" understory and "pure" exposed conditions.
+
+    NOTE: can't average the resistances; must convert to conductances,
+    average the conductances, and then convert the average conductance
+    to a resistance.
+
+    NOTE 2: since a resistance of 0 corresponds to an infinite conductance,
+    if either Ra_veg (resistance under veg) or Ra_bare (resistance over
+    exposed soil) are 0, then Ra_used must necessarily be 0 as well.
+  *************************************************/
+  if (veg_var->vegcover < 1) {
+    Ra_veg[0] = Ra_used[0];
+    Ra_veg[1] = Ra_used[1];
+    Ra_veg[2] = Ra_used[2];
+    /** If Ra_veg is non-zero, use it to compute area-weighted average **/
+    if (Ra_veg[0] > 0) {
+      /** aerodynamic conductance under vegetation **/
+      ga_veg = 1/Ra_veg[0];
+      /** compute aerodynamic resistance over exposed soil (Ra_bare) **/
+      tmp_wind[0] = wind[0];
+      tmp_wind[1] = ERROR; // unused
+      tmp_wind[2] = ERROR; // unused
+      tmp_height = soil_con->rough/RATIO_RL_HEIGHT_VEG;
+      tmp_displacement[0] = calc_veg_displacement(tmp_height);
+      tmp_roughness[0] = soil_con->rough;
+      tmp_ref_height[0] = WINDH_SOIL; // wind measurement height over bare soil
+      Error = CalcAerodynamic(0,0,0,soil_con->snow_rough,soil_con->rough,0,Ra_bare,tmp_wind,tmp_displacement,tmp_ref_height,tmp_roughness);
+      Ra_bare[0] /= StabilityCorrection(tmp_ref_height[0], tmp_displacement[0], TMean, Tair, tmp_wind[0], tmp_roughness[0]);
+      /** if Ra_bare is non-zero, compute area-weighted average
+          aerodynamic conductance **/
+      if (Ra_bare[0] > 0) {
+        /** aerodynamic conductance over exposed soil **/
+        ga_bare = 1/Ra_bare[0];
+        /** area-weighted average aerodynamic conductance **/
+        ga_average = veg_var->vegcover*ga_veg + (1-veg_var->vegcover)*ga_bare;
+        /** aerodynamic resistance is inverse of conductance **/
+        Ra_used[0] = 1/ga_average;
+      }
+      /** else aerodynamic resistance is zero **/
+      else {
+        Ra_used[0] = 0;
+      }
+    }
+    /** else aerodynamic resistance is zero **/
+    else {
+      Ra_used[0] = 0;
+    }
+  }
+
   /*************************************************
     Compute Evapotranspiration if not snow covered
 
     Should evapotranspiration be active when the 
     ground is only partially covered with snow????
 
-    Use Arno Evap if LAI is set to zero (e.g. no
-    winter crop planted).
+    Use Arno Evap in the exposed soil portion, and/or
+    if LAI is zero.
   *************************************************/
   if ( VEG && !SNOWING && veg_var->vegcover > 0 ) {
     Evap = canopy_evap(layer, veg_var, TRUE, 
@@ -709,20 +769,11 @@ double func_surf_energy_bal(double Ts, va_list ap)
         transp[i] = layer[i].evap;
         layer[i].evap = 0;
       }
-      tmp_wind[0] = wind[0];
-      tmp_wind[1] = -999.;
-      tmp_wind[2] = -999.;
-      tmp_height = soil_con->rough/0.123;
-      tmp_displacement[0] = calc_veg_displacement(tmp_height);
-      tmp_roughness[0] = soil_con->rough;
-      tmp_ref_height[0] = 10;
-      Error = CalcAerodynamic(0,0,0,soil_con->snow_rough,soil_con->rough,0,Ra_bare,tmp_wind,tmp_displacement,tmp_ref_height,tmp_roughness);
-      Ra_bare[0] /= StabilityCorrection(tmp_ref_height[0], tmp_displacement[0], TMean, Tair, tmp_wind[0], tmp_roughness[0]);
       Evap *= veg_var->vegcover;
       Evap += (1-veg_var->vegcover)
 	       * arno_evap(layer, surf_atten*NetBareRad, Tair, vpd, 
 		       depth[0], max_moist * depth[0] * 1000., 
-		       elevation, b_infilt, Ra_bare[0], delta_t, 
+		       elevation, b_infilt, Ra_used[0], delta_t, 
 		       resid_moist[0], frost_fract);
       for (i=0; i<options.Nlayer; i++) {
         layer[i].evap = veg_var->vegcover*transp[i] + (1-veg_var->vegcover)*layer[i].evap;

src/global.h

@@ -22,7 +22,7 @@
   2013-Dec-27 Removed QUICK_FS option.					TJB
   2014-May-20 Added ref_veg_vegcover.					TJB
 **********************************************************************/
-char *version = "4.2.b 2015-January-22";
+char *version = "4.2.c 2015-December-12";
 char *optstring = "g:vo";
 int flag;
 
@@ -65,4 +65,3 @@ param_set_struct param_set;
   double ref_veg_trunk_ratio[] = { 0.0, 0.0, 0.0, 0.0 };
   /* One element for each PET type (non-natural or natural) */
   char ref_veg_ref_crop[] = { FALSE, FALSE, TRUE, TRUE, FALSE, FALSE };
-