Address review comment and bug fix

-Use consistently `!` for functions that are mutating, and return always  `nothing` for these functions.
- Fix missing `canonicalize` from `Dates`

src/Wflow.jl

@@ -19,7 +19,8 @@ using Dates:
     DateTime,
     now,
     isleapyear,
-    datetime2unix
+    datetime2unix,
+    canonicalize
 using DelimitedFiles: readdlm
 using FieldMetadata: @metadata
 using Glob: glob
@@ -170,7 +171,7 @@ include("states.jl")
 """
     run(tomlpath::AbstractString; silent=false)
     run(config::Config)
-    run(model::Model)
+    run!(model::Model)
     run()
 
 Run an entire simulation starting either from a path to a TOML settings file,
@@ -211,6 +212,7 @@ function run(tomlpath::AbstractString; silent = nothing)
             close(logfile)
         end
     end
+    return nothing
 end
 
 function run(config::Config)
@@ -225,21 +227,22 @@ function run(config::Config)
     else
         error("unknown model type")
     end
-    load_fixed_forcing(model)
-    return run(model)
+    load_fixed_forcing!(model)
+    run!(model)
+    return model
 end
 
-function run_timestep(model::Model; update_func = update, write_model_output = true)
+function run_timestep!(model::Model; update_func = update!, write_model_output = true)
     advance!(model.clock)
     load_dynamic_input!(model)
-    model = update_func(model)
+    update_func(model)
     if write_model_output
         write_output(model)
     end
-    return model
+    return nothing
 end
 
-function run(model::Model; close_files = true)
+function run!(model::Model; close_files = true)
     (; config, writer, clock) = model
 
     model_type = config.model.type::String
@@ -260,7 +263,7 @@ function run(model::Model; close_files = true)
     runstart_time = now()
     @progress for (i, time) in enumerate(times)
         @debug "Starting timestep." time i now()
-        model = run_timestep(model)
+        run_timestep!(model)
     end
     @info "Simulation duration: $(canonicalize(now() - runstart_time))"
 
@@ -287,7 +290,7 @@ function run(model::Model; close_files = true)
             cp(src, dst; force = true)
         end
     end
-    return model
+    return nothing
 end
 
 function run()

src/bmi.jl

@@ -32,7 +32,7 @@ function BMI.initialize(::Type{<:Model}, config_file)
     else
         error("unknown model type")
     end
-    load_fixed_forcing(model)
+    load_fixed_forcing!(model)
     return model
 end
 
@@ -45,15 +45,15 @@ Update the model for a single timestep.
 """
 function BMI.update(model::Model; run = nothing)
     if isnothing(run)
-        model = run_timestep(model)
+        run_timestep!(model)
     elseif run == "sbm_until_recharge"
-        model = run_timestep(
+        run_timestep!(
             model;
-            update_func = update_until_recharge,
+            update_func = update_until_recharge!,
             write_model_output = false,
         )
     elseif run == "sbm_after_subsurfaceflow"
-        model = run_timestep(model; update_func = update_after_subsurfaceflow)
+        run_timestep!(model; update_func = update_after_subsurfaceflow!)
     end
     return model
 end
@@ -72,7 +72,7 @@ function BMI.update_until(model::Model, time::Float64)
         error(error_message)
     end
     for _ in 1:steps
-        model = run_timestep(model)
+        run_timestep!(model)
     end
     return model
 end
@@ -400,7 +400,7 @@ end
 # Extension of BMI functions (state handling and start time), required for OpenDA coupling.
 # May also be useful for other external software packages.
 function load_state(model::Model)
-    model = set_states(model)
+    set_states!(model)
     return model
 end
 

src/demand/water_demand.jl

@@ -199,6 +199,7 @@ function update_demand_gross!(nonpaddy::NonPaddy, soil::SbmSoilModel)
         end
         nonpaddy.variables.demand_gross[i] = irri_dem_gross
     end
+    return nothing
 end
 update_demand_gross!(nonpaddy::NoIrrigationNonPaddy, soil::SbmSoilModel) = nothing
 
@@ -324,6 +325,7 @@ function evaporation!(model::Paddy, potential_evaporation)
             model.variables.evaporation[i] = evaporation
         end
     end
+    return nothing
 end
 evaporation!(model::NoIrrigationPaddy, potential_evaporation) = nothing
 
@@ -345,6 +347,7 @@ function update_runoff!(model::Paddy, runoff)
             runoff[i] = paddy_runoff
         end
     end
+    return nothing
 end
 update_runoff!(model::NoIrrigationPaddy, runoff) = nothing
 
@@ -389,6 +392,7 @@ function update_demand_gross!(model::Paddy)
         end
         demand_gross[i] = irri_dem_gross
     end
+    return nothing
 end
 
 update_demand_gross!(paddy::NoIrrigationPaddy) = nothing
@@ -578,13 +582,14 @@ function return_flow_fraction!(model::NonIrrigationDemand)
     (; returnflow_fraction) = model.variables
     (; demand_gross, demand_net) = model.demand
     @. returnflow_fraction = return_flow_fraction(demand_gross, demand_net)
+    return nothing
 end
 
 # return zero (gross water demand) if non-irrigation water demand sector is not defined
 return_flow_fraction!(model::NoNonIrrigationDemand) = nothing
 
 "Update water allocation for river and land domains based on local surface water (river) availability."
-function surface_water_allocation_local(land_allocation, demand, river, network, dt)
+function surface_water_allocation_local!(land_allocation, demand, river, network, dt)
     (; surfacewater_alloc) = land_allocation.variables
     (; surfacewater_demand) = demand.variables
     (; act_surfacewater_abst_vol, act_surfacewater_abst, available_surfacewater) =
@@ -616,10 +621,11 @@ function surface_water_allocation_local(land_allocation, demand, river, network,
             surfacewater_alloc[i] = abstraction
         end
     end
+    return nothing
 end
 
 "Update water allocation for river and land domains based on surface water (river) availability for allocation areas."
-function surface_water_allocation_area(land_allocation, demand, river, network)
+function surface_water_allocation_area!(land_allocation, demand, river, network)
     inds_river = network.river.indices_allocation_areas
     inds_land = network.land.indices_allocation_areas
     res_index = network.river.reservoir_index
@@ -680,10 +686,11 @@ function surface_water_allocation_area(land_allocation, demand, river, network)
             surfacewater_alloc[j] += frac_allocate_sw * surfacewater_demand[j]
         end
     end
+    return nothing
 end
 
 "Update water allocation for land domain based on local groundwater availability."
-function groundwater_allocation_local(land_allocation, demand, groundwater_volume, network)
+function groundwater_allocation_local!(land_allocation, demand, groundwater_volume, network)
     (;
         surfacewater_alloc,
         act_groundwater_abst_vol,
@@ -712,10 +719,11 @@ function groundwater_allocation_local(land_allocation, demand, groundwater_volum
             groundwater_alloc[i] = abstraction
         end
     end
+    return nothing
 end
 
 "Update water allocation for land domain based on groundwater availability for allocation areas."
-function groundwater_allocation_area(land_allocation, demand, network)
+function groundwater_allocation_area!(land_allocation, demand, network)
     inds_river = network.river.indices_allocation_areas
     inds_land = network.land.indices_allocation_areas
     (;
@@ -754,6 +762,7 @@ function groundwater_allocation_area(land_allocation, demand, network)
             groundwater_alloc[j] += frac_allocate_gw * groundwater_demand[j]
         end
     end
+    return nothing
 end
 
 "Return and update non-irrigation sector (domestic, livestock, industry) return flow"
@@ -815,9 +824,9 @@ function update_water_allocation!(land_allocation, demand::Demand, lateral, netw
         frac_sw_used * nonirri_demand_gross + frac_sw_used * irri_demand_gross
 
     # local surface water demand and allocation (river, excluding reservoirs and lakes)
-    surface_water_allocation_local(land_allocation, demand, river, network, dt)
+    surface_water_allocation_local!(land_allocation, demand, river, network, dt)
     # surface water demand and allocation for areas
-    surface_water_allocation_area(land_allocation, demand, river, network)
+    surface_water_allocation_area!(land_allocation, demand, river, network)
 
     @. river.abstraction = act_surfacewater_abst_vol / dt
 
@@ -838,14 +847,14 @@ function update_water_allocation!(land_allocation, demand::Demand, lateral, netw
     act_groundwater_abst_vol .= 0.0
     act_groundwater_abst .= 0.0
     # local groundwater demand and allocation
-    groundwater_allocation_local(
+    groundwater_allocation_local!(
         land_allocation,
         demand,
         groundwater_volume(lateral.subsurface),
         network.land,
     )
     # groundwater demand and allocation for areas
-    groundwater_allocation_area(land_allocation, demand, network)
+    groundwater_allocation_area!(land_allocation, demand, network)
 
     # irrigation allocation
     for i in eachindex(total_alloc)
@@ -895,6 +904,8 @@ function update_demand_gross!(demand::Demand)
     @. nonirri_demand_gross = industry_dem + domestic_dem + livestock_dem
     @. total_gross_demand =
         nonpaddy_dem_gross + paddy_dem_gross + industry_dem + domestic_dem + livestock_dem
+
+    return nothing
 end
 
 update_demand_gross!(demand::NoDemand) = nothing
@@ -917,5 +928,7 @@ function update_water_demand!(demand::Demand, soil)
     update_demand_gross!(nonpaddy, soil)
     update_demand_gross!(paddy)
     update_demand_gross!(demand)
+
+    return nothing
 end
 update_water_demand!(demand::NoDemand, soil) = nothing