Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Timestep and Calendar Modifications #188

Merged
merged 13 commits into from
Jan 21, 2015
Merged

Timestep and Calendar Modifications #188

merged 13 commits into from
Jan 21, 2015

Conversation

jhamman
Copy link
Member

@jhamman jhamman commented Jan 12, 2015

Initial mods for subhourly timestep and alternative calendar support

  • Timestep units changed from hours to seconds,
  • Ported netcdftime.py to vic_time.c supporting all netCDF calendars.

Addresses #42 and #70 and contributes to #179.

To do (in upcoming commits):

  • Port to image driver
  • Test for all calendar and time units options
  • Test wide range of variable combinations in initialize_atmos.c
  • Code review from @bartnijssen
  • Code review from @tbohn

- Timestep units changed from hours to seconds,
- Ported netcdftime.py to vic_time.c supporting all netCDF calendars.

Not stable.  Timestep/Calendar commit num. 1.
@jhamman
Copy link
Member Author

jhamman commented Jan 13, 2015

Note: the Travis build failed only for driver/image which I have not made mods to yet.

@bartnijssen
Copy link
Member

@jhamman
In
https://github.com/jhamman/VIC/blob/feature/subhourly_timestep/vic_run/include/vic_physical_constants.h

put parentheses around every mathematical operation that is done in a #define statement.

For example,

#define MIN_PER_DAY MIN_PER_HOUR * HOURS_PER_DAY

should read

#define MIN_PER_DAY (MIN_PER_HOUR * HOURS_PER_DAY)

Similarly:

#define MIN_PER_DAY (MIN_PER_HOUR * HOURS_PER_DAY)  /**< hours per day */
#define SEC_PER_HOUR (SEC_PER_MIN * MIN_PER_HOUR)  /**< seconds per hour */
#define SEC_PER_DAY (SEC_PER_HOUR * HOURS_PER_DAY)  /**< hours per day */
#define CONST_OMEGA (2.0 * CONST_PI / CONST_SDAY)
#define CONST_RGAS (CONST_AVOGAD * CONST_BOLTZ) /**< Universal gas constant ~ J/K/kmole */
#define CONST_RDAIR (CONST_RGAS / CONST_MWDAIR)  /**< Dry air gas constant ~ J/K/kg */
#define CONST_RWV (CONST_RGAS / CONST_MWWV)  /**< Water vapor gas constant ~ J/K/kg */
#define CONST_EPS (CONST_MWWV / CONST_MWAIR)  /**< Ratio of molecular weights */
#define CONST_TSTD (CONST_TKFRZ + 15.0)  /**< (K) standard temp at 0.0 m elevation */
#define CONST_RHODAIR (CONST_PSTD / (CONST_RDAIR * CONST_TKFRZ))  /**< density of dry air at STP ~ kg/m^3 */

We should similarly check the other header files.

Keep in mind that the precompiler does textual substitution for #define. So

day = (double)date->day + (double)date->dayseconds / (double)SEC_PER_DAY;

becomes

day = (double)date->day + (double)date->dayseconds / (double) SEC_PER_HOUR * HOURS_PER_DAY;

which becomes

day = (double)date->day + (double)date->dayseconds / (double) SEC_PER_MIN * MIN_PER_HOUR * HOURS_PER_DAY;

which becomes

day = (double)date->day + (double)date->dayseconds / (double) 60 * 60 * 24;

which becomes

day = (double)date->day + (double)date->dayseconds * 24;

which is a big number for any dayseconds > 0. You just have to be careful with the define statements.

Note that if you put the parentheses around the statements, you will have

day = (double)date->day + (double)date->dayseconds / (double) ((SEC_PER_MIN * MIN_PER_HOUR) * HOURS_PER_DAY);

which will give you the right answer.

int i;
size_t stepsperday;

stepsperday = SEC_PER_DAY / stepsize;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this require extra checking to make sure that stepsize evenly divides SEC_PER_DAY. If stepsperday is a quantity that is regularly used, it may make sense to add it to the time structure for VIC.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a good point. I'm going to pass a stepsperday variable into the function and then calculate the stepsize inside the function.

In the global_param_struct, I store these values:

    double dt;                     /**< Time step in seconds */
    double snow_dt;                /**< Snow model time step in seconds */
    double runoff_dt;              /**< Runoff time step in seconds */
    double out_dt;                 /**< Output time step in seconds */
    size_t model_steps_per_day;    /**< Number of model timesteps per day */
    size_t snow_steps_per_day;     /**< Number of snow timesteps per day */
    size_t runoff_steps_per_day;   /**< Number of runoff timesteps per day */

but those do not necessarily correspond to the stepsize of this function.

@jhamman
Copy link
Member Author

jhamman commented Jan 13, 2015

I see. I see. That is much clearer now and I will implement that change shortly.

31, /* DECEMBER */
};
unsigned tmpstartdate;
unsigned tmpenddate;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Be explicit about the size of unsigned. I know it defaults to unsigned int, but since the rest of VIC is explicit about that, stick with it.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done. will extend to other locations as well.

@bartnijssen
Copy link
Member

There are a lot of whitespace-only changes in the commit that should not be included. Are you using the same uncrustify configuration file as we used before? You should rollback most of the whitespace only changes, cause they just clutter the commit and since uncrustify was already run on the original branch we should not include whitespace only changes.

short local_starthour;
unsigned short day_in_year, year, month, days_in_month;
int sec_offset_int;
unsigned tmp_endsec;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

unsigned int

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

will be done via #194

@jhamman jhamman added this to the 5.0 milestone Jan 15, 2015
@jhamman jhamman self-assigned this Jan 15, 2015
@bartnijssen
Copy link
Member

The history file is opened (and time dimension defined) in:
initialize_history_file()
https://github.com/bartnijssen/VIC/blob/develop/drivers/image/src/vic_init_output.c

For the state file in (note this does not currently have a time dimension):
initialize_state_file()
https://github.com/bartnijssen/VIC/blob/develop/drivers/image/src/vic_store.c

For the reads, the time dimension is currently not treated explicitly (just read one slice at a time from the beginning of the file, which obviously needs work). I just assume that the time slice that is being read in vic_force() in https://github.com/bartnijssen/VIC/blob/develop/drivers/image/src/vic_force.c increments by 1 every time (dstart[0] increases, the rest is constant). Proper treatment of the calendar would mean that dstart[0] is calculated in a better way (I just start at 0 and increment).

On Jan 15, 2015, at 12:29 PM, Joe Hamman notifications@github.com wrote:

@bartnijssen , I added the initial mods for the image driver (5d2ad84). I think we will still want to update the netCDF time dimensions to use the new calendar functionality in vic_time(). I had trouble tracking down what you had put into the time variable in the history files so I left it as is for now. If you can point me to the right place, I'll add that as another commit.


Reply to this email directly or view it on GitHub.

…rly_timestep

Conflicts:
	drivers/image/src/get_global_param.c
	drivers/shared/include/vic_driver_shared.h
@jhamman
Copy link
Member Author

jhamman commented Jan 18, 2015

Thanks @bartnijssen. In order to to keep this PR to a manageable size, I'm going to leave the updating of the image driver I/O time handling to another issue (#197) and PR .

@jhamman
Copy link
Member Author

jhamman commented Jan 19, 2015

I've completed my modifications for this branch. I've tested on a few grid cells with the classic driver only but it is running stably for all calendars. Remaining issues for the image driver have been opened so I think this can be merged then closed.

@bartnijssen
Copy link
Member

it is running stably for all calendars

Does that mean all the mtclim operations work correctly in all the different calendars (skipping leap days, etc) or simply that the model does not crash if we specify different calendars? The main concern with the different calendars is that this creates a lot of new options (and hence test pathways). If this has not been fully tested with mtclim (e.g. would that handle all the different calendars correctly - which is not the same as not crashing), then we need to disallow the calendars that have not been tested or only allow the calendars when all the forcings are fully specified.

@jhamman
Copy link
Member Author

jhamman commented Jan 20, 2015

@bartnijssen:

That is a good point. The model does not crash with the different calendars.

Maybe @tbohn has some ideas on the applicability of other calendars to mtclim. I know there are a number of places where the code uses a hard coded year length but I seem to remember @tbohn thinking it would be ok for the standard, all_leap, and no_leap calendars. I think it would be an issue for the 360_day calendar.

I'm fine with restricting the use of other calendars to cases when mtclim is not used.

@bartnijssen
Copy link
Member

Even a single leap day difference can add up to a month in a 100 year simulation, which is why this can come back to haunt us.

We'll want to include at the very least a warning with non-standard calendars that they have not been tested. I'd rather not rely on someone thinking it'd probably be OK, but prefer to do real testing. In the absence of that (I do not suggest making this a major push right now), keep the code, but warn the user and open a new issue for more thorough testing of this feature.

Bart

On Jan 20, 2015, at 10:46 AM, Joe Hamman notifications@github.com wrote:

@bartnijssen:

That is a good point. The model does not crash with the different calendars.

Maybe @tbohn has some ideas on the applicability of other calendars to mtclim. I know there are a number of places where the code uses a hard coded year length but I seem to remember @tbohn thinking it would be ok for the standard, all_leap, and no_leap calendars. I think it would be an issue for the 360_day calendar.

I'm fine with restricting the use of other calendars to cases when mtclim is not used.


Reply to this email directly or view it on GitHub.

@tbohn
Copy link
Contributor

tbohn commented Jan 20, 2015

Yes, I think the standard, no_leap and all_leap calendars would be fine
with MTCLIM - I believe MTCLIM uses a no_leap calendar and we just
duplicate a day for leap years.

I think we could even find a way to make the 360-day calendar work, if we
realy wanted to. But I am guessing that only the CCSM system uses it? And
if so, the forcings would already be supplied to VIC at sub-daily
timescales, so MTCLIM would be bypassed.

Ted

On Tue, Jan 20, 2015 at 10:46 AM, Joe Hamman notifications@github.com
wrote:

@bartnijssen https://github.com/bartnijssen:

That is a good point. The model does not crash with the different
calendars.

Maybe @tbohn https://github.com/tbohn has some ideas on the
applicability of other calendars to mtclim. I know there are a number of
places where the code uses a hard coded year length but I seem to remember
@tbohn https://github.com/tbohn thinking it would be ok for the standard,
all_leap, and no_leap calendars. I think it would be an issue for the
360_day calendar.

I'm fine with restricting the use of other calendars to cases when mtclim
is not used.


Reply to this email directly or view it on GitHub
#188 (comment).

@jhamman
Copy link
Member Author

jhamman commented Jan 20, 2015

A cumulative leap day offset should not be a problem since mtclim only knows which day in the year [1, 366]. The only problem I see in mtclim is when there is a calculation that uses something like day_in_year / total_days_in_year. For example in calc_srad_humidity_iterative():

 /* begin loop through yeardays */
for (i = 0; i < DAYS_PER_YEAR; i++) {
    ...
     /* solar constant as a function of yearday (W/m^2) */
    sc = param.MTCLIM_SOLAR_CONSTANT + 45.5 *
    sin((2.0 * CONST_PI * (double)i / CONST_DDAYS_PER_YEAR) + 1.7);
    ....
}

In previous versions of mtclim, this was insensitive to whether or not this was a leap year.

@bartnijssen
Copy link
Member

@jhamman : If you add the small changes to the MTCLIM code that we discussed (pass both day of year and the number of days in the year), then I'll go ahead and merge

@bartnijssen
Copy link
Member

@jhamman : BTW - also note that the pull request cannot be merged automatically, because I merged another pull request first. If you can resolve the conflict on your end that would be great.

@jhamman
Copy link
Member Author

jhamman commented Jan 20, 2015

No problem, I'll bring this branch up to date before committing these last changes.

Joe Hamman added 3 commits January 20, 2015 17:26
…s a parameter "days_per_year" to determine if leap days should be included.
…rly_timestep

Conflicts:
	drivers/shared/include/vic_driver_shared.h
@jhamman
Copy link
Member Author

jhamman commented Jan 21, 2015

Ok, with this last set of commits I

  • brought the branch up to date with the develop branch
  • applied the mtclim calendar fix we discussed today. This results in a no-change outcome for the standard, gregorian, and julian calendars where as the no_leap, all_leap and 360_day calendars are now handled in a better way.

The build passed and I verified in a few test cases that this is doing what we want.

bartnijssen added a commit that referenced this pull request Jan 21, 2015
@bartnijssen bartnijssen merged commit 6849a5b into UW-Hydro:develop Jan 21, 2015
This was referenced Jan 21, 2015
@jhamman jhamman deleted the feature/subhourly_timestep branch January 22, 2015 18:15
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants