Feature/debug cesm driver

docs/Development/ReleaseNotes.md

@@ -197,6 +197,16 @@ 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).
 
+6. 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 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_init_output.c

@@ -224,8 +224,8 @@ initialize_history_file(nc_file_struct *nc,
     double                    *dvar;
 
 
-    // This could be further refined but for now, I've choosen a file naming
-    // Convention that goes like this:
+    // This could be further refined but for now, I've chosen a file naming
+    // convention that goes like this:
     switch (stream->agg_alarm.freq) {
     // If FREQ_NDAYS -- filename = result_dir/prefix.YYYY-MM-DD.nc
     case FREQ_NDAYS:

vic/drivers/shared_image/src/vic_store.c

@@ -48,25 +48,45 @@ 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
     if (mpi_rank == VIC_MPI_ROOT) {
+
+	// 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);
 
@@ -350,7 +370,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;
@@ -376,7 +396,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;
@@ -813,7 +833,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;
@@ -826,7 +846,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;
@@ -972,7 +992,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 */