Change basic time unit to seconds about vic

[jhamman]

I'm in the middle of working through this feature and could use a bit of input on a few things.

First, I've already completed changing the timestep inside vic_run. That process was fairly straight forward and the dt variable passed throughout the run step is now a double precision float with units of seconds.

The real fun comes in the drivers. What I'd like to bring up for discussion is the constraints we want to put on a sub hourly timestep. In other words, what are acceptable values for dt? Below are the minimum constraints that I have in place for the timestep:

// Timestep Assertions:
assert ( global.dt <= SEC_PER_DAY );  // Maximum timestep is 1day
if ( global.dt < SEC_PER_DAY ) {  // Sub-daily timestep
    assert ( SEC_PER_DAY % global.dt == 0 );  // Timestep is evenly divisible into a day
    assert ( global.dt <= 6*SEC_PER_HOUR );  // Timestep is less than 6hrs
}
assert ( global.dt >= global.snow_dt == 0 );  // Snow timestep is less than or equal to the model timestep
assert ( global.dt % global.snow_dt == 0 );  // Snow timestep is evenly divisible into the model timestep
assert ( global.snow_dt <= 6*SEC_PER_HOUR );  // Snow timestep is less than 6hrs

assert ( global.dt >= MINIMUM_DT );  // Minimum timestep (e.g. 60 seconds)
assert ( global.snow_dt >= MINIMUM_DT );  // Minimum snow timestep (e.g. 60 seconds)
if ( global.dt < SEC_PER_HOUR ) { // Sub-hourly timestep
    assert ( SEC_PER_HOUR % global.dt == 0 );  // Timestep is evenly divisible into an hour
}

These last three are not entirely necessary but would allow for some assumptions to be made about output structures, etc.

Am I missing anything here that we think would be useful?

[bartnijssen]

This will depend on the driver and so we would have different checks for each driver (so isolating this into a check_timestep(double dt) function would be useful. This would be driver specific. For example, the CESM driver won't have separate global.dt and global.snow_dt variables (or the variables will always have the same value) and there should be few restrictions on what dt is in that case.

In general, you need to be careful because if global.dt is a double, operations like (SEC_PER_HOUR % global.dt) may not work the way you anticipate, nor can you really do an '==' operation. You'll need to typecast it to a size_t or something like that (unless we want to handle sub-second timesteps. I don't think we need to).

In addition, we should not hardcode the 6 as the maximum snow time step in hours (just put a maximum step in seconds in a header file as MAXIMUM_SUBDAILY_DT or something like it).

I'd also group similar checks (dt within certain bounds, e.g.

    assert ( global.dt <= SEC_PER_DAY );  // Maximum timestep is 1day
    assert ( global.dt >= MINIMUM_DT );  // Minimum timestep (e.g. 60 seconds)

[jhamman]

@bartnijssen, thanks for your above comments, that helps keep me moving on this. In get_global_param.c, I have a temporary unsigned short that I use to validate the timesteps, I then type cast it to a double when populating the global.dt field.

One question has recently come up in vic_run/src/runoff.c#L192. In this portion of the code, there is currently a fixed hourly timestep for calculating the flow between the soil layers. How do we want to handle this with a subhourly timestep? One option is to only use a "sub-dt" timestep if the model is run at a daily or > 1hr timestep.

[bartnijssen]

@jhamman : I'd be careful casting a double to an unsigned short, since they are basically incompatible (size wise). Even though by the time you do this, you probably already tested that dt <= 86400, so the unsigned short should be large enough to hold the dt value, it would be more consistent to use a size_t.

The runoff part is interesting, because of the explicit assumption of a time step that is an integer number of hours. We must already handle this in the RASM code, which probably has a much older version of runoff.c. I wonder whether this occurs in other routines in vic_run() as well. If nothing else, it would be good to isolate all the dt related changes as a separate branch, so that we can check that we are not changing the answer by switching (or at least we know why we are changing the answer).

A quick grep -i hour vic_run/*.c shows that there are a fair number of cases where hour is assumed or implied. Meaning that a change from int to double for the timestep is only the first step and making this a separate branch is important.

I know that this does not answer the original question of how to handle the sub-hourly timestep in runoff.c. @tbohn : Could you comment on why this method was adopted (you may not have worked on this either).

[jhamman]

@bartnijssen : I have all my timestep changes on my feature/subhourly_timestep branch. Once I push my changes to github, I'll be asking for a code review from you and @tbohn.

Looking at the RASM version of runoff.c, we basically are just setting the "runoff" dt to the "model" dt:

if (dt < 1.0) {
    dt = 1.0;
}

dt_inflow = inflow / (double) dt;

for (time_step = 0; time_step < dt; time_step++) {
    inflow = dt_inflow;

...

So one option would be to do basically the same thing when global.dt < SEC_PER_HOUR. The more explicit option would be to set a model parameter global.runoff_dt where (global.dt % global.runoff_dt == 0 && global.runoff_dt <= SEC_PER_HOUR).

[bartnijssen]

That looks extremely dodgy to me.

First I don't like the

    if (dt < 1.0) {
        dt = 1.0;
    }

Since dt really shouldn't be set locally but is the timestep of the model. We should copy it to a more descriptive variable in that case. Also, using a non-integer as the upper-bound for a loop counter is not a good idea. It would read a lot better if it said something along the lines of

    if (dt < SECS_PER_HOUR) {
        n_local_time_steps = 1;
    }

Although this is still rather confusing.

Perhaps more importantly, I am concerned about the numerical implementation, i.e. what is so special about an hourly timestep (note that this used to be the smallest possible time step in VIC and that may be why it was selected). But maybe we need to shelve that for future discussion (and open a new issue for that) rather than make that part of the change in the time step units from hours to seconds.

I understand what the intention was: If you run at a 24 hour time step, then you could saturate the top layer if you put all the moisture in at once. In that case the amount of fast response runoff (among other things) would become a function of the model time step. By running runoff() at a fixed (hourly) timestep, you avoid this to some extent.

I think there may be a few other instances like this in the code.

[jhamman]

closed via #188.