Merge pull request #677 from yixinmao/develop

Merge VIC5.0.1 release into develop branch

docs/Development/ReleaseNotes.md

@@ -14,6 +14,30 @@ To check which release of VIC you are running:
 - For VIC 5 and later, type `vic_{classic,image}.exe -v`
 
 ------------------------------
+
+## VIC 5.1.0
+
+**Release date: (Unreleased)**
+
+#### Model enhancement:
+
+1. Improved calculation of drainage between soil layers ([GH#656](https://github.com/UW-Hydro/VIC/pull/656))
+
+	Drainage from upper layer to adjacent lower layer is calculated according to Brook & Corey curve (where drainage rate is a function of upper-layer soil moisture). In previous versions, a simple numerical solution is applied which uses the timestep-beginning upper-layer soil moisture to calculate drainage rate, and assume this constant rate over the entire timestep. This can cause unreasonably large drainage if the curve has a steep shape and when soil moisture is high. Now, the current version uses exact integral (instead of numerical solution) for layer drainage calculation.
+
+2. Fixes for the CESM driver ([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`
+
+------------------------------
+
 ## VIC 5.0.1
 
 **Release date: (February 1, 2017)**
@@ -61,7 +85,7 @@ To check which release of VIC you are running:
 
 ------------------------------
 
-## VIC 5.0.0 [![DOI](https://zenodo.org/badge/7766/UW-Hydro/VIC.svg)](https://zenodo.org/badge/latestdoi/7766/UW-Hydro/VIC)
+## VIC 5.0.0 [![DOI](https://zenodo.org/badge/doi/10.5281/zenodo.61422.svg)](http://dx.doi.org/10.5281/zenodo.61422)
 
 **Release date: (September 2, 2016)**
 
@@ -236,6 +260,7 @@ This is a major update from VIC 4. The VIC 5.0.0 release aims to have nearly ide
 
 	Fixed a bug where volumetric heat capacity of water should be used in `func_canopy_energy_bal` (previously specific heat capacity was used).
 
+
 ------------------------------
 
 ## VIC 4.2.d [![DOI](https://zenodo.org/badge/doi/10.5281/zenodo.56058.svg)](http://dx.doi.org/10.5281/zenodo.56058)

vic/drivers/cesm/src/cesm_interface_c.c

@@ -96,6 +96,9 @@ vic_cesm_init(vic_clock     *vclock,
     // populate model state, either using a cold start or from a restart file
     vic_populate_model_state(trim(cmeta->starttype));
 
+    // initialize forcings
+    vic_force();
+
     // initialize output structures
     vic_init_output(&dmy_current);
 
@@ -137,6 +140,7 @@ vic_cesm_run(vic_clock *vclock)
 
     // if save:
     if (vclock->state_flag) {
+        // write state file
         vic_store(&dmy_current, state_filename);
         write_rpointer_file(state_filename);
     }

vic/drivers/cesm/src/cesm_put_data.c

@@ -140,12 +140,12 @@ vic_cesm_put_data()
 
                 // 2m reference temperature
                 // CESM units: K
-                l2x_vic[i].l2x_Sl_tref += AreaFactor * force->air_temp[NR];
+                l2x_vic[i].l2x_Sl_tref += AreaFactor * (force[i].air_temp[NR] + CONST_TKFRZ);
 
                 // 2m reference specific humidity
                 // CESM units: g/g
                 l2x_vic[i].l2x_Sl_qref += AreaFactor * CONST_EPS *
-                                          force->vp[NR] / force->pressure[NR];
+                                          force[i].vp[NR] / force[i].pressure[NR];
 
                 // Albedo Note: VIC does not partition its albedo, all returned
                 // values will be the same
@@ -156,10 +156,10 @@ vic_cesm_put_data()
                 // force->NetShortAtmos net shortwave flux (+ down)
                 // SWup = force->shortwave[NR] - energy.NetShortAtmos
                 // Set the albedo to zero for the case where there is no shortwave down
-                if (force->shortwave[NR] > 0.) {
+                if (force[i].shortwave[NR] > 0.) {
                     albedo = AreaFactor *
-                             (force->shortwave[NR] - energy.NetShortAtmos) /
-                             force->shortwave[NR];
+                             (force[i].shortwave[NR] - energy.NetShortAtmos) /
+                             force[i].shortwave[NR];
                 }
                 else {
                     albedo = 0.;
@@ -184,7 +184,7 @@ vic_cesm_put_data()
 
                 // 10m wind
                 // CESM units: m/s
-                l2x_vic[i].l2x_Sl_u10 += AreaFactor * force->wind[NR];
+                l2x_vic[i].l2x_Sl_u10 += AreaFactor * force[i].wind[NR];
 
                 // dry deposition velocities (optional)
                 // CESM units: ?
@@ -259,7 +259,7 @@ vic_cesm_put_data()
                 // upward longwave heat flux
                 // CESM units: W m-2
                 l2x_vic[i].l2x_Fall_lwup += AreaFactor *
-                                            (force->longwave[NR] -
+                                            (force[i].longwave[NR] -
                                              energy.NetLongAtmos);
 
                 // evaporation water flux
@@ -278,7 +278,7 @@ vic_cesm_put_data()
 
                 // heat flux shortwave net
                 l2x_vic[i].l2x_Fall_swnet += AreaFactor *
-                                             (force->shortwave[NR] -
+                                             (force[i].shortwave[NR] -
                                               energy.NetShortAtmos);
 
                 // co2 flux **For testing set to 0

vic/drivers/cesm/src/get_global_param.c

@@ -385,6 +385,12 @@ get_global_param(FILE *gp)
                             "control file.");
                 }
             }
+            else if (strcasecmp("SNOW_BAND", optstr) == 0) {
+	        sscanf(cmdstr, "%*s %s", flgstr);
+		if (str_to_bool(flgstr)) {
+                    options.SNOW_BAND = SNOW_BAND_TRUE_BUT_UNSET;
+                }
+            }
             else if (strcasecmp("LAKES", optstr) == 0) {
                 sscanf(cmdstr, "%*s %s", flgstr);
                 if (strcasecmp("FALSE", flgstr) == 0) {
@@ -415,6 +421,9 @@ get_global_param(FILE *gp)
             else if (strcasecmp("OUTVAR", optstr) == 0) {
                 ; // do nothing
             }
+            else if (strcasecmp("AGGFREQ", optstr) == 0) {
+		; // do nothing
+            }
             else if (strcasecmp("OUTPUT_STEPS_PER_DAY", optstr) == 0) {
                 ; // do nothing
             }

vic/drivers/cesm/src/vic_cesm_start.c

@@ -46,6 +46,11 @@ vic_cesm_start(vic_clock     *vclock,
     // Driver specific settings
     if (mpi_rank == VIC_MPI_ROOT) {
         strcpy(filenames.global, GLOBALPARAM);
+
+	// assign case name to state file name
+    	strncpy(filenames.statefile, trim(cmeta->caseid),
+                sizeof(filenames.statefile));
+
         // read global settings
         filep.globalparam = open_file(filenames.global, "r");
         get_global_param(filep.globalparam);

vic/drivers/cesm/src/vic_force.c

@@ -127,8 +127,8 @@ vic_force(void)
     for (j = 0; j < NF; j++) {
         for (i = 0; i < local_domain.ncells_active; i++) {
             // CESM units: Pa
-            // VIC units: kPa
-            force[i].pressure[j] = x2l_vic[i].x2l_Sa_pbot / PA_PER_KPA;
+            // VIC units: Pa
+            force[i].pressure[j] = x2l_vic[i].x2l_Sa_pbot;
         }
     }
 

vic/drivers/shared_image/src/vic_store.c

@@ -48,25 +48,44 @@ vic_store(dmy_struct *dmy_current,
     size_t                     k;
     size_t                     m;
     size_t                     p;
+    double                     offset;
+    double                     time_num;
+    double                     end_time_num;
     int                       *ivar = NULL;
     double                    *dvar = NULL;
     size_t                     d2start[2];
     size_t                     d3start[3];
     size_t                     d4start[4];
     size_t                     d5start[5];
     size_t                     d6start[6];
+    dmy_struct                 end_time_date;
     nc_file_struct             nc_state_file;
     nc_var_struct             *nc_var;
 
     set_nc_state_file_info(&nc_state_file);
 
-    // only open and initialize the netcdf file on the first thread
+    // advance dmy_current by one timestep since dmy_current is the 
+    // timestep-beginning timestamp, and state file date should be 
+    // the end of the current time step 
+    dt_seconds_to_time_units(global_param.time_units, global_param.dt,
+                             &offset);
+    time_num = date2num(global_param.time_origin_num, dmy_current, 0,
+                        global_param.calendar, global_param.time_units);
+    end_time_num = time_num + offset;
+
+    // allocate dmy struct for end of current time step 
+    num2date(global_param.time_origin_num, end_time_num, 0.,
+                             global_param.calendar, global_param.time_units,
+                             &end_time_date);
+
     // create netcdf file for storing model state
     sprintf(filename, "%s.%04i%02i%02i_%05u.nc",
-            filenames.statefile, global_param.stateyear,
-            global_param.statemonth, global_param.stateday,
-            global_param.statesec);
+            filenames.statefile, end_time_date.year,
+            end_time_date.month, end_time_date.day,
+            end_time_date.dayseconds);
+
     initialize_state_file(filename, &nc_state_file, dmy_current);
+
     if (mpi_rank == VIC_MPI_ROOT) {
         debug("writing state file: %s", filename);
     }
@@ -348,7 +367,7 @@ vic_store(dmy_struct *dmy_current,
             }
             gather_put_nc_field_int(nc_state_file.nc_id,
                                     nc_var->nc_varid,
-                                    nc_state_file.d_fillvalue,
+                                    nc_state_file.i_fillvalue,
                                     d4start, nc_var->nc_counts, ivar);
             for (i = 0; i < local_domain.ncells_active; i++) {
                 ivar[i] = nc_state_file.i_fillvalue;
@@ -374,7 +393,7 @@ vic_store(dmy_struct *dmy_current,
             }
             gather_put_nc_field_int(nc_state_file.nc_id,
                                     nc_var->nc_varid,
-                                    nc_state_file.d_fillvalue,
+                                    nc_state_file.i_fillvalue,
                                     d4start, nc_var->nc_counts, ivar);
             for (i = 0; i < local_domain.ncells_active; i++) {
                 ivar[i] = nc_state_file.i_fillvalue;
@@ -811,7 +830,7 @@ vic_store(dmy_struct *dmy_current,
         }
         gather_put_nc_field_int(nc_state_file.nc_id,
                                 nc_var->nc_varid,
-                                nc_state_file.d_fillvalue,
+                                nc_state_file.i_fillvalue,
                                 d2start, nc_var->nc_counts, ivar);
         for (i = 0; i < local_domain.ncells_active; i++) {
             ivar[i] = nc_state_file.i_fillvalue;
@@ -824,7 +843,7 @@ vic_store(dmy_struct *dmy_current,
         }
         gather_put_nc_field_int(nc_state_file.nc_id,
                                 nc_var->nc_varid,
-                                nc_state_file.d_fillvalue,
+                                nc_state_file.i_fillvalue,
                                 d2start, nc_var->nc_counts, ivar);
         for (i = 0; i < local_domain.ncells_active; i++) {
             ivar[i] = nc_state_file.i_fillvalue;
@@ -970,7 +989,7 @@ vic_store(dmy_struct *dmy_current,
         }
         gather_put_nc_field_int(nc_state_file.nc_id,
                                 nc_var->nc_varid,
-                                nc_state_file.d_fillvalue,
+                                nc_state_file.i_fillvalue,
                                 d2start, nc_var->nc_counts, ivar);
         for (i = 0; i < local_domain.ncells_active; i++) {
             ivar[i] = nc_state_file.i_fillvalue;

vic/drivers/shared_image/src/vic_write.c

@@ -194,7 +194,7 @@ vic_write(stream_struct  *stream,
                 }
                 gather_put_nc_field_schar(nc_hist_file->nc_id,
                                           nc_hist_file->nc_vars[k].nc_varid,
-                                          nc_hist_file->d_fillvalue,
+                                          nc_hist_file->c_fillvalue,
                                           dstart, dcount, cvar);
             }
             else {

vic/vic_run/include/vic_def.h

@@ -52,7 +52,7 @@
 /***** Model Constants *****/
 #define MAXSTRING    2048
 #define MISSING      -99999.   /**< missing value */
-#define MISSING_USI  99999.    /**< missing value for unsigned ints */
+#define MISSING_USI  99999    /**< missing value for unsigned ints */
 #define MISSING_S    "MISSING"    /**< missing value for strings */
 #define NODATA_VH    -1        /**< missing value for veg_hist inputs */
 #define NODATA_VEG   -1        /**< flag for veg types not in grid cell */

vic/vic_run/include/vic_run.h

@@ -113,6 +113,7 @@ void compute_pot_evap(size_t, double, double, double, double, double, double,
                       double, double *);
 void compute_runoff_and_asat(soil_con_struct *, double *, double, double *,
                              double *);
+double calc_Q12(double, double, double, double, double);
 void compute_soil_resp(int, double *, double, double, double *, double *,
                        double, double, double, double *, double *, double *);
 void compute_soil_layer_thermal_properties(layer_data_struct *, double *,

vic/vic_run/src/runoff.c

@@ -200,12 +200,10 @@ runoff(cell_data_struct  *cell,
                 }
 
                 if (liq[lindex] > resid_moist[lindex]) {
-                    Q12[lindex] = Ksat[lindex] *
-                                  pow(((tmp_liq -
-                                        resid_moist[lindex]) /
-                                       (soil_con->max_moist[lindex] -
-                                        resid_moist[lindex])),
-                                      soil_con->expt[lindex]);
+                    Q12[lindex] = calc_Q12(Ksat[lindex], tmp_liq,
+                                           resid_moist[lindex],
+                                           soil_con->max_moist[lindex],
+                                           soil_con->expt[lindex]);
                 }
                 else {
                     Q12[lindex] = 0.;
@@ -487,3 +485,20 @@ compute_runoff_and_asat(soil_con_struct *soil_con,
         *runoff = 0.;
     }
 }
+
+/******************************************************************************
+* @brief    Calculate drainage between two layers
+******************************************************************************/
+double
+calc_Q12(double Ksat, double init_moist, double resid_moist,
+         double max_moist, double expt)
+{
+    double Q12;
+
+    Q12 = init_moist - pow(pow(init_moist - resid_moist, 1.0 - expt) -
+          Ksat / pow(max_moist - resid_moist, expt) * (1.0 - expt),
+          1.0 / (1.0 - expt)) - resid_moist;
+
+    return Q12;
+}
+