Merge pull request #712 from dgergel/feature/initialize_albedo_cesm_d…

…river_2

Feature/initialize albedo cesm driver

docs/Development/ReleaseNotes.md

@@ -29,36 +29,40 @@ To check which release of VIC you are running:
 
 	[GH#642](https://github.com/UW-Hydro/VIC/pull/642)
 
-	1. Using correct fill value datatypes in MPI Gather steps
-	2. Updated state file name time step to be period-ending rather than period-beginning
-	3. Set the state file name to the RASM case ID
-	4. Removed decimal point for missing values for unsigned integers
-	5. Create dummy forcings when initializing the model (so that there is forcing data for the first time step)
-	6. Changed pressure units from kPa to Pa
-	7. Fixed bug that prevented using the correct local domain grid cells in `cesm_put_data.c`
-	8. Changed reference temperature units from Celsius to Kelvin in `cesm_put_data.c`
+	    1. Using correct fill value datatypes in MPI Gather steps
+	    2. Updated state file name time step to be period-ending rather than period-beginning
+	    3. Set the state file name to the RASM case ID
+	    4. Removed decimal point for missing values for unsigned integers
+	    5. Create dummy forcings when initializing the model (so that there is forcing data for the first time step)
+	    6. Changed pressure units from kPa to Pa
+	    7. Fixed bug that prevented using the correct local domain grid cells in `cesm_put_data.c`
+	    8. Changed reference temperature units from Celsius to Kelvin in `cesm_put_data.c`
 
 	[GH#695](https://github.com/UW-Hydro/VIC/pull/695)
 
-	1. Fix sign for latent heat fluxes passed from VIC to the coupler
-	2. Fix sign for longwave radiation passed from VIC to the coupler
+	    1. Fix sign for latent heat fluxes passed from VIC to the coupler
+	    2. Fix sign for longwave radiation passed from VIC to the coupler
 
 	[GH#696](https://github.com/UW-Hydro/VIC/pull/696)
 
-	1. Changes names of CESM driver functions `trim` and `advance_time` to `trimstr` and `advance_vic_time`, respectively, to avoid conflicts with WRF functions with the same names when compiling RFR case.
+	    Changes names of CESM driver functions `trim` and `advance_time` to `trimstr` and `advance_vic_time`, respectively, to avoid conflicts with WRF functions with the same names when compiling RFR case.
 
-    [GH#702] (https://github.com/UW-Hydro/VIC/pull/702)
+        [GH#702] (https://github.com/UW-Hydro/VIC/pull/702)
 
-    1. Fixes Julian day for the first timestep in the dmy struct for the CESM driver.
+	    Fixes Julian day for the first timestep in the dmy struct for the CESM driver.
 
 3. Speed up NetCDF operations in the image/CESM drivers ([GH#684](https://github.com/UW-Hydro/VIC/pull/684))
 
-	These changes speed up image driver initialization, forcing reads, and history writes by only opening and closing each input netCDF file once.
+    These changes speed up image driver initialization, forcing reads, and history writes by only opening and closing each input netCDF file once.
 
 4. Added two new timers to measure time in I/O operations ([GH#703](https://github.com/UW-Hydro/VIC/pull/703))
 
     These two timers count the CPU and WALL time spent in ``vic_force`` and ``vic_write``. The accumulated time from these timers is printed out at the end of each simulation in the timing table. See also [GH#442](https://github.com/UW-Hydro/VIC/pull/442).
 
+5. Added gridcell-averaged albedo (STATE_AVG_ALBEDO) as a state file variable ([GH#712](https://github.com/UW-Hydro/VIC/pull/712))
+
+    This is for use in the CESM driver for VIC to pass to WRF, but has been implemented in the core structure of VIC (in vic_run) for consistency with the classic and image drivers. Running VIC from a cold start now also includes calculation of gridcell-averaged albedo.
+
 ## VIC 5.0.1
 
 **Release date: (February 1, 2017)**

vic/drivers/cesm/include/vic_driver_cesm.h

@@ -187,5 +187,6 @@ void vic_cesm_run_model(void);
 void vic_cesm_start(vic_clock *vclock, case_metadata *cmeta);
 void vic_populate_model_state(char *runtype_str);
 void write_rpointer_file(char *fname);
+void vic_initialize_albedo(void);
 
 #endif

vic/drivers/cesm/src/cesm_interface_c.c

@@ -102,6 +102,9 @@ vic_cesm_init(vic_clock     *vclock,
     // initialize output structures
     vic_init_output(&dmy_current);
 
+    // initialize albedo
+    vic_initialize_albedo();
+
     // initialization is complete, print settings
     log_info(
         "Initialization is complete, print global param and options structures");

vic/drivers/cesm/src/vic_cesm_init_library.c

@@ -122,3 +122,25 @@ initialize_l2x_data(void)
         l2x_vic[i].l2x_vars_set = true;
     }
 }
+
+/******************************************************************************
+ * @brief      Initialize albedo values in l2x_data_struct.
+ *****************************************************************************/
+void
+vic_initialize_albedo(void)
+{
+    extern l2x_data_struct *l2x_vic;
+    extern domain_struct    local_domain;
+    extern all_vars_struct *all_vars;
+
+    size_t                  i;
+
+    log_info("Initializing albedo values");
+
+    for (i = 0; i < local_domain.ncells_active; i++) {
+        l2x_vic[i].l2x_Sl_avsdr = all_vars[i].gridcell_avg.avg_albedo;
+        l2x_vic[i].l2x_Sl_anidr = all_vars[i].gridcell_avg.avg_albedo;
+        l2x_vic[i].l2x_Sl_avsdf = all_vars[i].gridcell_avg.avg_albedo;
+        l2x_vic[i].l2x_Sl_anidf = all_vars[i].gridcell_avg.avg_albedo;
+    }
+}

vic/drivers/shared_all/include/vic_driver_shared_all.h

@@ -354,6 +354,7 @@ enum
     STATE_LAKE_ICE_SNOW_PACK_WATER,    /**<  lake ice snow pack water: lake_var.pack_water */
     STATE_LAKE_ICE_SNOW_ALBEDO,        /**<  lake ice snow albedo: lake_var.SAlbedo */
     STATE_LAKE_ICE_SNOW_DEPTH,         /**<  lake ice snow depth: lake_var.sdepth */
+    STATE_AVG_ALBEDO,                  /**<  gridcell-averaged albedo: gridcell_avg.avg_albedo */
     // Last value of enum - DO NOT ADD ANYTHING BELOW THIS LINE!!
     // used as a loop counter and must be >= the largest value in this enum
     N_STATE_VARS                       /**< used as a loop counter*/

vic/drivers/shared_all/src/generate_default_state.c

@@ -53,15 +53,21 @@ generate_default_state(all_vars_struct *all_vars,
     };
     double                   Cv;
     double                   tmp;
+    double                   AreaFactor;
+    double                   TreeAdjustFactor = 1.;
+    double                   lakefactor = 1.;
+    double                   albedo_sum;
     double                ***tmpT;
     double                 **tmpZ;
     int                      ErrorFlag;
 
     cell_data_struct       **cell;
     energy_bal_struct      **energy;
+    veg_var_struct         **veg_var;
 
     cell = all_vars->cell;
     energy = all_vars->energy;
+    veg_var = all_vars->veg_var;
     Nveg = veg_con[0].vegetat_type_num;
 
     // allocate memory for tmpT and tmpZ
@@ -125,6 +131,34 @@ generate_default_state(all_vars_struct *all_vars,
         }
     }
 
+
+    /************************************************************************
+       Initialize gridcell-averaged albedo
+    ************************************************************************/
+    // vegetation class-weighted albedo over gridcell
+    albedo_sum = 0;
+    for (veg = 0; veg <= Nveg; veg++) {
+        Cv = veg_con[veg].Cv;
+        if (Cv > 0) {
+            for (band = 0; band < options.SNOW_BAND; band++) {
+                if (soil_con->AreaFract[band] > 0.) {
+                    // TO-DO: account for treeline and lake factors
+                    AreaFactor = (Cv * soil_con->AreaFract[band] *
+                                  TreeAdjustFactor * lakefactor);
+                    // cold start, so assuming bare (free of snow) albedo
+                    if (veg != Nveg) {
+                        albedo_sum += AreaFactor * veg_var[veg][band].albedo;
+                    }
+                    else {
+                        // this is the bare soil class, so use bare soil albedo parameter
+                        albedo_sum += AreaFactor * param.ALBEDO_BARE_SOIL;
+                    }
+                }
+            }
+        }
+    }
+    all_vars->gridcell_avg.avg_albedo = albedo_sum;
+
     /************************************************************************
        Initialize soil layer ice content
     ************************************************************************/

vic/drivers/shared_all/src/vic_history.c

@@ -30,7 +30,7 @@
  * @brief    This routine creates the list of output data.
  *****************************************************************************/
 void
-alloc_out_data(size_t     ngridcells,
+alloc_out_data(size_t    ngridcells,
                double ***out_data)
 {
     extern metadata_struct out_metadata[N_OUTVAR_TYPES];

vic/drivers/shared_image/src/state_metadata.c

@@ -275,6 +275,17 @@ set_state_meta_data_info()
     strcpy(state_metadata[STATE_ENERGY_SNOW_FLUX].description,
            "thermal flux through snowpack");
 
+    // STATE_GRIDCELL_AVG_ALBEDO
+    strcpy(state_metadata[STATE_AVG_ALBEDO].varname,
+           "STATE_AVG_ALBEDO");
+    strcpy(state_metadata[STATE_AVG_ALBEDO].long_name,
+           "state_avg_albedo");
+    strcpy(state_metadata[STATE_AVG_ALBEDO].standard_name,
+           "state_gridcell_avg_albedo");
+    strcpy(state_metadata[STATE_AVG_ALBEDO].units, "fraction");
+    strcpy(state_metadata[STATE_AVG_ALBEDO].description,
+           "gridcell averaged albedo");
+
     if (options.LAKES) {
         // STATE_LAKE_SOIL_MOISTURE
         strcpy(state_metadata[STATE_LAKE_SOIL_MOISTURE].varname,

vic/drivers/shared_image/src/vic_alloc.c

@@ -86,7 +86,7 @@ vic_alloc(void)
     check_alloc_status(out_data, "Memory allocation error.");
 
     // save_data allocation
-    save_data = malloc(local_domain.ncells_active * sizeof(*save_data));
+    save_data = calloc(local_domain.ncells_active, sizeof(*save_data));
     check_alloc_status(save_data, "Memory allocation error.");
 
     // allocate memory for individual grid cells

vic/drivers/shared_image/src/vic_restore.c

@@ -461,6 +461,14 @@ vic_restore(void)
         }
     }
 
+    // grid cell-averaged albedo: gridcell_avg.avg_albedo
+    get_scatter_nc_field_double(&(filenames.init_state),
+                                state_metadata[STATE_AVG_ALBEDO].varname,
+                                d2start, d2count, dvar);
+    for (i = 0; i < local_domain.ncells_active; i++) {
+        all_vars[i].gridcell_avg.avg_albedo = dvar[i];
+    }
+
     // soil node temperatures: energy[veg][band].T[nidx]
     for (m = 0; m < options.NVEGTYPES; m++) {
         d5start[0] = m;

vic/drivers/shared_image/src/vic_store.c

@@ -725,6 +725,19 @@ vic_store(dmy_struct *dmy_state,
         }
     }
 
+    // Grid cell averaged albedo
+    nc_var = &(nc_state_file.nc_vars[STATE_AVG_ALBEDO]);
+    for (i = 0; i < local_domain.ncells_active; i++) {
+        dvar[i] = (double) all_vars[i].gridcell_avg.avg_albedo;
+    }
+    gather_put_nc_field_double(nc_state_file.nc_id,
+                               nc_var->nc_varid,
+                               nc_state_file.d_fillvalue,
+                               d2start, nc_var->nc_counts, dvar);
+    for (i = 0; i < local_domain.ncells_active; i++) {
+        dvar[i] = nc_state_file.d_fillvalue;
+    }
+
 
     if (options.LAKES) {
         // total soil moisture
@@ -1423,6 +1436,15 @@ set_nc_state_var_info(nc_file_struct *nc)
             nc->nc_vars[i].nc_counts[3] = nc->nj_size;
             nc->nc_vars[i].nc_counts[4] = nc->ni_size;
             break;
+        case STATE_AVG_ALBEDO:
+            // 2d vars [j, i]
+            nc->nc_vars[i].nc_dims = 2;
+            nc->nc_vars[i].nc_dimids[0] = nc->nj_dimid;
+            nc->nc_vars[i].nc_dimids[1] = nc->ni_dimid;
+            nc->nc_vars[i].nc_counts[0] = nc->nj_size;
+            nc->nc_vars[i].nc_counts[1] = nc->ni_size;
+            break;
+
         case STATE_LAKE_SOIL_MOISTURE:
             // 3d vars [layer, j, i]
             nc->nc_vars[i].nc_dims = 3;

vic/vic_run/include/vic_def.h

@@ -999,6 +999,14 @@ typedef struct {
                                condensation from snow pack (m) */
 } snow_data_struct;
 
+/******************************************************************************
+ * @brief   This structures stores variables averaged over a grid cell
+ *****************************************************************************/
+typedef struct {
+    // Grid cell averaged variables
+    double avg_albedo;            /**< Average albedo over a grid cell */
+} gridcell_avg_struct;
+
 /******************************************************************************
  * @brief   This structure stores the lake/wetland parameters for a grid cell
  *****************************************************************************/
@@ -1085,6 +1093,7 @@ typedef struct {
     lake_var_struct lake_var;     /**< Stores lake/wetland variables */
     snow_data_struct **snow;      /**< Stores snow variables */
     veg_var_struct **veg_var;     /**< Stores vegetation variables */
+    gridcell_avg_struct gridcell_avg;   /**< Stores gridcell average variables */
 } all_vars_struct;
 
 #endif

vic/vic_run/include/vic_run.h

@@ -49,6 +49,8 @@ double calc_atmos_energy_bal(double, double, double, double, double, double,
                              double, double *, double *, double *, double *,
                              double *, double *, bool *, unsigned int*);
 double calc_density(double);
+void calc_gridcell_avg_albedo(double *, double, size_t, energy_bal_struct **,
+                              veg_con_struct *, soil_con_struct *);
 double calc_latent_heat_of_sublimation(double temp);
 double calc_latent_heat_of_vaporization(double temp);
 int calc_layer_average_thermal_props(energy_bal_struct *, layer_data_struct *,

vic/vic_run/src/calc_gridcell_avg_albedo.c

@@ -0,0 +1,75 @@
+/******************************************************************************
+ * @section DESCRIPTION
+ *
+ * This routine computes gridcell-averaged albedo.
+ *
+ * Note: this routine is specifically designed for the CESM driver, for WRF,
+ * but has been implemented in other drivers as well for the sake of consistency.
+ *
+ * @section LICENSE
+ *
+ * The Variable Infiltration Capacity (VIC) macroscale hydrological model
+ * Copyright (C) 2016 The Computational Hydrology Group, Department of Civil
+ * and Environmental Engineering, University of Washington.
+ *
+ * The VIC model is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *****************************************************************************/
+
+#include <vic_run.h>
+
+/******************************************************************************
+ * @brief    Compute gridcell-averaged albedo.
+ *****************************************************************************/
+void
+calc_gridcell_avg_albedo(double             *albedo,
+                         double              shortwave,
+                         size_t              Nveg,
+                         energy_bal_struct **energy,
+                         veg_con_struct     *veg_con,
+                         soil_con_struct    *soil_con)
+{
+    extern option_struct options;
+    size_t               veg;
+    size_t               band;
+    double               Cv;
+    double               AreaFactor;
+    double               TreeAdjustFactor = 1.;
+    double               lakefactor = 1;
+    double               swnet;
+
+    swnet = 0;
+
+    for (veg = 0; veg <= Nveg; veg++) {
+        Cv = veg_con[veg].Cv;
+        if (Cv > 0) {
+            for (band = 0; band < options.SNOW_BAND; band++) {
+                if (soil_con->AreaFract[band] > 0.) {
+                    // TO-DO: account for treeline and lake factors
+                    AreaFactor = (Cv * soil_con->AreaFract[band] *
+                                  TreeAdjustFactor * lakefactor);
+                    swnet += energy[veg][band].NetShortAtmos * AreaFactor;
+                }
+            }
+        }
+    }
+
+    // compute gridcell-averaged albedo using average shortwave
+    if (shortwave > 0) {
+        *albedo = 1. - (swnet / shortwave);
+    }
+    else {
+        *albedo = 0.0;
+    }
+}

vic/vic_run/src/vic_run.c

@@ -382,6 +382,11 @@ vic_run(force_data_struct   *force,
         }
     }
 
+    // Compute gridcell-averaged albedo
+    calc_gridcell_avg_albedo(&all_vars->gridcell_avg.avg_albedo,
+                             force->shortwave[NR], Nveg, energy,
+                             veg_con, soil_con);
+
     /****************************
        Run Lake Model
     ****************************/