Update to code for Esoil in sparse vegetation.

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/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) {

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/vicNl_def.h

@@ -314,11 +314,19 @@ extern char ref_veg_ref_crop[];
 #define B_SVP 17.269
 #define C_SVP 237.3
 
+/* define constants for vegetation turbulence parameters */
+#define RATIO_RL_HEIGHT_VEG 0.123 /* Ratio of roughness length (m) to vegetation height (m) */
+#define RATIO_DH_HEIGHT_VEG 0.67  /* Ratio of displacement height (m) to vegetation height (m) */
+
 /* define constants for canopy resistance */
 #define CLOSURE 4000		/* Threshold vapor pressure deficit for stomatal closure (Pa) */
 #define RSMAX 5000              /* Maximum allowable resistance (s/m) */
 #define VPDMINFACTOR 0.1        /* Minimum allowable vapor pressure deficit factor */
 
+/* define constants for soil evaporation */
+#define RARC_SOIL 100.0         /* Architectural resistance (s/m) of soil when computing soil evaporation via Penman-Monteith eqn */
+#define WINDH_SOIL 10.0         /* Default wind measurement height over soil (m) */
+
 /* define constants for penman evaporation */
 #define CP_PM 1013		/* specific heat of moist air at constant pressure (J/kg/C)
 				   (Handbook of Hydrology) */