Merge pull request #262 from Deltares/starttime

Starttime

Author: vers-w

docs/src/changelog.md

@@ -17,6 +17,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 - The time values returned in the BMI interface are no longer in seconds since 1970, but in
   seconds since the model start time. This is more in line with standard BMI practices.
+- The `starttime` was defined one model timestep `Δt` ahead of the actual model time (the
+  initial conditions timestamp (state time)). As a consequence this was also the case for
+  the current model time. To allow for an easier interpretation of Wflow time handling,
+  either through BMI or directly, the `starttime` is now equal to the state time, resulting
+  in current model times without an offset.
 
 ### Added
 - For (regulated) lakes with rating curve of type 1 (H-Q table), lake `storage` above the

docs/src/user_guide/step2_settings_file.md

@@ -6,15 +6,19 @@ The filepaths that are provided in this file are relative to the location of the
 or to `dir_input` and `dir_output` if they are given.
 
 ## General time info
-Time information is optional. When left out, each time step in the forcing NetCDF will be
-calculated. If you wish to calculate a subset of this time range, or a different timestep,
-you can specify a `starttime`, `endtime` and `timestepsecs` yourself. The `time_units`
-optional information is used by the `writer` of the model, for model output in netCDF
-format. The `calendar` option allows you to calculate in one of the different [CF
-conventions calendars](http://cfconventions.org/cf-conventions/cf-conventions.html#calendar)
-provided by the [CFTime.jl package](https://juliageo.org/CFTime.jl/latest/), such as
-`"360_day"`. This is useful if you want to calculate climate scenarios which are sometimes
-provided in these alternative calendars.
+Time information is optional. When left out, for each timestamp in the forcing NetCDF Wflow
+will do computations, except for the first forcing timestamp that is considered equal to the
+initial conditions of the Wflow model (state time). If you wish to calculate a subset of
+this time range, or a different timestep, you can specify a `starttime`, `endtime` and
+`timestepsecs` yourself. The `starttime` is defined as the model state time. In the TOML
+file settings below the `starttime` is 2000-01-01T00:00:00 (state time) and the first update
+(and output) of the Wflow model is at 2000-01-02T00:00:00. The `time_units` optional
+information is used by the `writer` of the model, for model output in netCDF format. The
+`calendar` option allows you to calculate in one of the different [CF conventions
+calendars](http://cfconventions.org/cf-conventions/cf-conventions.html#calendar) provided by
+the [CFTime.jl package](https://juliageo.org/CFTime.jl/latest/), such as `"360_day"`. This
+is useful if you want to calculate climate scenarios which are sometimes provided in these
+alternative calendars.
 
 ```toml
 calendar = "standard"                           # optional, this is default value

src/Wflow.jl

@@ -40,40 +40,41 @@ function Clock(config)
     calendar = get(config, "calendar", "standard")::String
     starttime = cftime(config.starttime, calendar)
     Δt = Second(config.timestepsecs)
-    Clock(starttime, 1, Δt)
+    Clock(starttime, 0, Δt)
 end
 
 function Clock(config, reader)
     nctimes = reader.dataset["time"][:]
-    # if the config file does not have a start or endtime, folow the NetCDF times
-    # and add them to the config
+    
     # if the timestep is not given, use the difference between NetCDF time 1 and 2
+    timestepsecs = get(config, "timestepsecs", nothing)
+    if timestepsecs === nothing
+        timestepsecs = Dates.value(Second(nctimes[2] - nctimes[1]))
+        config.timestepsecs = timestepsecs
+    end
+    Δt = Second(timestepsecs)
+    
+    # if the config file does not have a start or endtime, follow the NetCDF times
+    # and add them to the config
     starttime = get(config, "starttime", nothing)
     if starttime === nothing
-        starttime = first(nctimes)
+        starttime = first(nctimes) - Δt
         config.starttime = starttime
     end
     endtime = get(config, "endtime", nothing)
     if endtime === nothing
         endtime = last(nctimes)
         config.endtime = endtime
     end
-    timestepsecs = get(config, "timestepsecs", nothing)
-    if timestepsecs === nothing
-        timestepsecs = Dates.value(Second(nctimes[2] - nctimes[1]))
-        config.timestepsecs = timestepsecs
-    end
 
     calendar = get(config, "calendar", "standard")::String
-    starttime = cftime(config.starttime, calendar)
-    Δt = Second(timestepsecs)
-
     fews_run = get(config, "fews_run", false)::Bool
     if fews_run
-        starttime = starttime + Δt
+        config.starttime = starttime + Δt
     end
-
-    Clock(starttime, 1, Δt)
+    starttime = cftime(config.starttime, calendar)
+    
+    Clock(starttime, 0, Δt)
 end
 
 include("io.jl")
@@ -192,31 +193,42 @@ function run(config::Config)
     run(model)
 end
 
+function run_timestep(model::Model; update_func = update, write_model_output = true)
+    advance!(model.clock)
+    load_dynamic_input!(model)
+    model = update_func(model)
+    if write_model_output
+        write_output(model)
+    end
+    return model
+end
+
 function run(model::Model; close_files = true)
     @unpack network, config, writer, clock = model
 
     model_type = config.model.type::String
 
     # determine timesteps to run
     calendar = get(config, "calendar", "standard")::String
+    @warn string(
+        "The definition of `starttime` has changed (equal to model state time).\n Please", 
+        " update your settings TOML file by subtracting one model timestep Δt from the", 
+        " `starttime`, if it was used with a Wflow version up to v0.6.3.",
+    )
     starttime = clock.time
     Δt = clock.Δt
     endtime = cftime(config.endtime, calendar)
-    times = range(starttime, endtime, step = Δt)
+    times = range(starttime + Δt, endtime, step = Δt)
 
     @info "Run information" model_type starttime Δt endtime nthreads()
     runstart_time = now()
     @progress for (i, time) in enumerate(times)
         @debug "Starting timestep." time i now()
-        load_dynamic_input!(model)
-        model = update(model)
+        model = run_timestep(model)
     end
     @info "Simulation duration: $(canonicalize(now() - runstart_time))"
 
     # write output state NetCDF
-    # undo the clock advance at the end of the last iteration, since there won't
-    # be a next step, and then the output state falls on the correct time
-    rewind!(clock)
     if haskey(config, "state") && haskey(config.state, "path_output")
         @info "Write output states to NetCDF file `$(model.writer.state_nc_path)`."
     end

src/bmi.jl

@@ -48,35 +48,34 @@ Update the model for a single timestep.
 - `run = nothing`: to update a model partially.
 """
 function BMI.update(model::Model; run = nothing)
-    @unpack network, config = model
+    @unpack clock, network, config = model
     if isnothing(run)
-        update_func = update
+        model = run_timestep(model)
     elseif run == "sbm_until_recharge"
-        update_func = update_until_recharge
+        model = run_timestep(
+            model,
+            update_func = update_until_recharge,
+            write_model_output = false,
+        )
     elseif run == "sbm_after_subsurfaceflow"
-        update_func = update_after_subsurfaceflow
+        model = run_timestep(model, update_func = update_after_subsurfaceflow)
     end
-    load_dynamic_input!(model)
-    return update_func(model)
+    return model
 end
 
 function BMI.update_until(model::Model, time::Float64)
-    @unpack network, config = model
+    @unpack clock, network, config = model
     curtime = BMI.get_current_time(model)
-    n_iter = Int(max(0, (time - curtime) / model.clock.Δt.value))
-    end_time = curtime + n_iter * config.timestepsecs
-    @info "Updating model until $end_time."
-    for _ = 1:n_iter
-        load_dynamic_input!(model)
-        update(model)
+    n = Int(max(0, (time - curtime) / model.clock.Δt.value))
+    for _ = 1:n
+        model = run_timestep(model)
     end
     return model
 end
 
 "Write state output to netCDF and close files."
 function BMI.finalize(model::Model)
     @unpack config, writer, clock = model
-    rewind!(clock)
     write_netcdf_timestep(model, writer.state_dataset, writer.state_parameters)
     reset_clock!(model.clock, config)
     close_files(model, delete_output = false)

src/flextopo_model.jl

@@ -720,10 +720,5 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:FlextopoModel}
 
     surface_routing(model)
 
-    write_output(model)
-
-    # update the clock
-    advance!(clock)
-
     return model
 end

src/hbv_model.jl

@@ -426,10 +426,5 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:HbvModel}
 
     surface_routing(model)
 
-    write_output(model)
-
-    # update the clock
-    advance!(clock)
-
     return model
 end

src/sbm_gwf_model.jl

@@ -476,10 +476,5 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
     ssf_toriver[inds_riv] = -lateral.subsurface.river.flux ./ lateral.river.Δt
     surface_routing(model, ssf_toriver = ssf_toriver)
 
-    write_output(model)
-
-    # update the clock
-    advance!(clock)
-
     return model
 end

src/sbm_model.jl

@@ -496,10 +496,5 @@ function update_after_subsurfaceflow(
     ssf_toriver = lateral.subsurface.to_river ./ tosecond(basetimestep)
     surface_routing(model, ssf_toriver = ssf_toriver)
 
-    write_output(model)
-
-    # update the clock
-    advance!(clock)
-
     return model
 end

src/sediment_model.jl

@@ -199,10 +199,5 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SedimentModel}
         update(lateral.river, network.river, config)
     end
 
-    write_output(model)
-
-    # update the clock
-    advance!(clock)
-
     return model
 end

test/bmi.jl

@@ -13,7 +13,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             @test BMI.get_time_step(model) == 86400.0
             @test BMI.get_start_time(model) == 0.0
             @test BMI.get_current_time(model) == 0.0
-            @test BMI.get_end_time(model) == 30 * 86400.0
+            @test BMI.get_end_time(model) == 31 * 86400.0
         end
 
         @testset "model information functions" begin

test/flextopo_config.toml

@@ -1,6 +1,6 @@
 casename = "wflow_meuse"
 calendar = "proleptic_gregorian"
-starttime = "2010-01-01T00:00:00"
+starttime = "2009-12-31T00:00:00"
 endtime = "2010-07-01T00:00:00"
 time_units = "days since 1900-01-01 00:00:00"
 timestepsecs = 86400

test/hbv_config.toml

@@ -5,7 +5,7 @@
 
 calendar = "proleptic_gregorian"
 endtime = 2000-02-01T00:00:00
-starttime = 2000-01-01T00:00:00
+starttime = 1999-12-31T00:00:00
 time_units = "days since 1900-01-01 00:00:00"
 timestepsecs = 86400
 

test/io.jl

@@ -18,7 +18,7 @@ config = Wflow.Config(tomlpath)
     @test dirname(config) == dirname(tomlpath)
 
     # test if the values are parsed as expected
-    @test config.starttime === DateTime(2000, 1, 2)
+    @test config.starttime === DateTime(2000, 1, 1)
     @test config.endtime === DateTime(2000, 2)
     @test config.output.path == "output_moselle.nc"
     @test config.output isa Wflow.Config
@@ -71,11 +71,11 @@ end
     pop!(Dict(config), "timestepsecs")
     clock = Wflow.Clock(config, reader)
 
-    @test clock.time == DateTimeProlepticGregorian(2000, 1, 2)
-    @test clock.iteration == 1
+    @test clock.time == DateTimeProlepticGregorian(2000, 1, 1)
+    @test clock.iteration == 0
     @test clock.Δt == Second(Day(1))
     # test that the missing keys have been added to the config
-    @test config.starttime == DateTime(2000, 1, 2)
+    @test config.starttime == DateTime(2000, 1, 1)
     @test config.endtime == DateTime(2001, 1, 1)
     @test config.timestepsecs == 86400
 
@@ -87,7 +87,7 @@ end
 
     clock = Wflow.Clock(config, reader)
     @test clock.time == DateTimeStandard(2003, 4, 5)
-    @test clock.iteration == 1
+    @test clock.iteration == 0
     @test clock.Δt == Second(Hour(1))
 
     close(ds)
@@ -98,43 +98,43 @@ end
     # 29 days in this February due to leap year
     starttime = DateTimeStandard(2000, 2, 28)
     Δt = Day(1)
-    clock = Wflow.Clock(starttime, 1, Second(Δt))
+    clock = Wflow.Clock(starttime, 0, Second(Δt))
 
     Wflow.advance!(clock)
     Wflow.advance!(clock)
     @test clock.time == DateTimeStandard(2000, 3, 1)
-    @test clock.iteration == 3
+    @test clock.iteration == 2
     @test clock.Δt == Δt
 
     Wflow.rewind!(clock)
     @test clock.time == DateTimeStandard(2000, 2, 29)
-    @test clock.iteration == 2
+    @test clock.iteration == 1
     @test clock.Δt == Δt
 
     config = Wflow.Config(
         Dict("starttime" => starttime, "timestepsecs" => Dates.value(Second(Δt))),
     )
     Wflow.reset_clock!(clock, config)
     @test clock.time == starttime
-    @test clock.iteration == 1
+    @test clock.iteration == 0
     @test clock.Δt == Δt
 end
 
 @testset "Clock{DateTime360Day}" begin
     # 30 days in each month
     starttime = DateTime360Day(2000, 2, 29)
     Δt = Day(1)
-    clock = Wflow.Clock(starttime, 1, Second(Δt))
+    clock = Wflow.Clock(starttime, 0, Second(Δt))
 
     Wflow.advance!(clock)
     Wflow.advance!(clock)
     @test clock.time == DateTime360Day(2000, 3, 1)
-    @test clock.iteration == 3
+    @test clock.iteration == 2
     @test clock.Δt == Δt
 
     Wflow.rewind!(clock)
     @test clock.time == DateTime360Day(2000, 2, 30)
-    @test clock.iteration == 2
+    @test clock.iteration == 1
     @test clock.Δt == Δt
 
     config = Wflow.Config(
@@ -147,7 +147,7 @@ end
     Wflow.reset_clock!(clock, config)
     @test clock.time isa DateTime360Day
     @test string(clock.time) == "2020-02-29T00:00:00"
-    @test clock.iteration == 1
+    @test clock.iteration == 0
     @test clock.Δt == Δt
 end
 
@@ -193,6 +193,7 @@ config.input["path_forcing"] = abs_path_forcing
 @test isabspath(config.input.path_forcing)
 
 model = Wflow.initialize_sbm_model(config)
+Wflow.advance!(model.clock)
 Wflow.load_dynamic_input!(model)
 
 @unpack vertical, clock, reader, writer = model
@@ -273,6 +274,7 @@ config.input.vertical.c = Dict(
 )
 
 model = Wflow.initialize_sbm_model(config)
+Wflow.advance!(model.clock)
 Wflow.load_dynamic_input!(model)
 
 @testset "changed parameter values" begin
@@ -429,11 +431,11 @@ end
     # test Clock{DateTimeNoLeap}
     clock = Wflow.Clock(config, reader)
     @test clock.time isa DateTimeNoLeap
-    @test clock.time == DateTimeNoLeap(2000, 1, 2)
+    @test clock.time == DateTimeNoLeap(2000, 1, 1)
 
     starttime = DateTimeNoLeap(2000, 2, 28)
     Δt = Day(1)
-    clock = Wflow.Clock(starttime, 1, Second(Δt))
+    clock = Wflow.Clock(starttime, 0, Second(Δt))
     Wflow.advance!(clock)
     @test clock.time == DateTimeNoLeap(2000, 3, 1)
 end