Skip to content

Commit

Permalink
Standard names model type sediment
Browse files Browse the repository at this point in the history
As the `sediment` model type makes use of hydrological forcing timeseries that are stored differently (internally) than the same timeseries used by model types `sbm` and `sbm_gwf`, the standard name mapping is done separately for the `sediment` model type (different dict). The same hydrological forcing data is used by the `SoilLoss` and `OverlandFlowSediment` components. This data is now only stored in the `SoilLoss` component and shared (memory) with the `OverlandFlowSediment` component.
  • Loading branch information
vers-w committed Dec 18, 2024
1 parent 5858906 commit 8b38fb2
Show file tree
Hide file tree
Showing 17 changed files with 520 additions and 794 deletions.
2 changes: 1 addition & 1 deletion src/Wflow.jl
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,6 @@ function Clock(config, reader)
return Clock(starttime, 0, dt)
end

include("standard_name.jl")
include("io.jl")

"""
Expand Down Expand Up @@ -181,6 +180,7 @@ include("erosion.jl")
include("sediment_flux.jl")
include("sediment_model.jl")
include("sbm_gwf_model.jl")
include("standard_name.jl")
include("utils.jl")
include("bmi.jl")
include("subdomains.jl")
Expand Down
34 changes: 18 additions & 16 deletions src/io.jl
Original file line number Diff line number Diff line change
Expand Up @@ -215,7 +215,7 @@ mover_params = (
)

function load_fixed_forcing!(model)
(; reader, network, config) = model
(; reader, vertical, network, config) = model
(; forcing_parameters) = reader

do_reservoirs = get(config.model, "reservoirs", false)::Bool
Expand All @@ -234,7 +234,7 @@ function load_fixed_forcing!(model)
for (par, ncvar) in forcing_parameters
if ncvar.name === nothing
val = ncvar.value * ncvar.scale + ncvar.offset
lens = standard_name_map[par]
lens = standard_name_map(vertical)[par]
param_vector = lens(model)
param_vector .= val
# set fixed precipitation and evaporation over the lakes and reservoirs and put
Expand All @@ -261,7 +261,7 @@ end

"Get dynamic netCDF input for the given time"
function update_forcing!(model)
(; clock, reader, network, config) = model
(; clock, reader, vertical, network, config) = model
(; dataset, dataset_times, forcing_parameters) = reader

do_reservoirs = get(config.model, "reservoirs", false)::Bool
Expand Down Expand Up @@ -311,7 +311,7 @@ function update_forcing!(model)
end
end
end
lens = standard_name_map[par]
lens = standard_name_map(vertical)[par]
param_vector = lens(model)
sel = active_indices(network, par)
data_sel = data[sel]
Expand Down Expand Up @@ -352,7 +352,7 @@ end

"Get cyclic netCDF input for the given time"
function update_cyclic!(model)
(; clock, reader, network) = model
(; clock, reader, vertical, network) = model
(; cyclic_dataset, cyclic_times, cyclic_parameters) = reader

# pick up the data that is valid for the past model time step
Expand All @@ -368,7 +368,7 @@ function update_cyclic!(model)
error("Could not find applicable cyclic timestep for $month_day")
# load from netCDF into the model according to the mapping
data = get_at(cyclic_dataset, ncvar.name, i)
lens = standard_name_map[par]
lens = standard_name_map(vertical)[par]
param_vector = lens(model)
sel = active_indices(network, par)
param_vector .= data[sel]
Expand Down Expand Up @@ -419,6 +419,7 @@ function setup_scalar_netcdf(
config,
float_type = Float32,
)
(; vertical) = modelmap
ds = create_tracked_netcdf(path)
defDim(ds, "time", Inf) # unlimited
defVar(
Expand All @@ -439,8 +440,8 @@ function setup_scalar_netcdf(
(nc.location_dim,);
attrib = ["cf_role" => "timeseries_id"],
)
v = if haskey(standard_name_map, nc.par)
lens = standard_name_map[nc.par]
v = if haskey(standard_name_map(vertical), nc.par)
lens = standard_name_map(vertical)[nc.par]
lens(modelmap)
else
param(modelmap, nc.par)
Expand Down Expand Up @@ -906,9 +907,10 @@ Create a Dict that maps parameter netCDF names to arrays in the Model.
"""
function out_map(ncnames_dict, modelmap)
output_map = Dict{String, Any}()
(; vertical) = modelmap
for (par, ncname) in ncnames_dict
A = if haskey(standard_name_map, par)
lens = standard_name_map[par]
A = if haskey(standard_name_map(vertical), par)
lens = standard_name_map(vertical)[par]
lens(modelmap)
else
param(modelmap, par)
Expand Down Expand Up @@ -1081,12 +1083,12 @@ end

"Write a new timestep with scalar data to a netCDF file"
function write_netcdf_timestep(model, dataset)
(; writer, clock, config) = model
(; writer, vertical, clock, config) = model

time_index = add_time(dataset, clock.time)
for (nt, nc) in zip(writer.nc_scalar, config.netcdf.variable)
A = if haskey(standard_name_map, nt.parameter)
lens = standard_name_map[nt.parameter]
A = if haskey(standard_name_map(vertical), nt.parameter)
lens = standard_name_map(vertical)[nt.parameter]
lens(model)
else
param(model, nt.parameter)
Expand Down Expand Up @@ -1336,13 +1338,13 @@ function reducer(col, rev_inds, x_nc, y_nc, config, dataset, fileformat)
end

function write_csv_row(model)
(; writer, clock, config) = model
(; writer, vertical, clock, config) = model
isnothing(writer.csv_path) && return nothing
io = writer.csv_io
print(io, string(clock.time))
for (nt, col) in zip(writer.csv_cols, config.csv.column)
A = if haskey(standard_name_map, nt.parameter)
lens = standard_name_map[nt.parameter]
A = if haskey(standard_name_map(vertical), nt.parameter)
lens = standard_name_map(vertical)[nt.parameter]
lens(model)
else
param(model, nt.parameter)
Expand Down
44 changes: 11 additions & 33 deletions src/parameters.jl
Original file line number Diff line number Diff line change
Expand Up @@ -86,16 +86,11 @@ function LandGeometry(nc, config, inds)
sizeinmetres = get(config.model, "sizeinmetres", false)::Bool
xl, yl = cell_lengths(y, cellength, sizeinmetres)
area = xl .* yl
ldd = ncread(nc, config, "ldd"; optional = false, sel = inds, allow_missing = true)
lens = lens_input("ldd")
ldd = ncread(nc, config, lens; optional = false, sel = inds, allow_missing = true)
drain_width = map(get_flow_width, ldd, xl, yl)
landslope = ncread(
nc,
config,
"vertical.land_parameter_set.slope";
optional = false,
sel = inds,
type = Float,
)
lens = lens_input_parameter("land_surface__slope")
landslope = ncread(nc, config, lens; optional = false, sel = inds, type = Float)
clamp!(landslope, 0.00001, Inf)

land_parameter_set =
Expand All @@ -115,30 +110,13 @@ end

"Initialize river geometry parameters"
function RiverGeometry(nc, config, inds)
riverwidth = ncread(
nc,
config,
"lateral.river_parameter_set.width";
optional = false,
sel = inds,
type = Float,
)
riverlength = ncread(
nc,
config,
"lateral.river_parameter_set.length";
optional = false,
sel = inds,
type = Float,
)
riverslope = ncread(
nc,
config,
"lateral.river_parameter_set.slope";
optional = false,
sel = inds,
type = Float,
)
lens = lens_input_parameter("river__width")
riverwidth = ncread(nc, config, lens; optional = false, sel = inds, type = Float)
lens = lens_input_parameter("river__length")
riverlength = ncread(nc, config, lens; optional = false, sel = inds, type = Float)
lens = lens_input_parameter("river__slope")
riverslope = ncread(nc, config, lens; optional = false, sel = inds, type = Float)

minimum(riverlength) > 0 || error("river length must be positive on river cells")
minimum(riverwidth) > 0 || error("river width must be positive on river cells")
clamp!(riverslope, 0.00001, Inf)
Expand Down
31 changes: 7 additions & 24 deletions src/sediment/erosion/overland_flow_erosion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -34,30 +34,13 @@ end

"Initialize ANSWERS overland flow erosion model parameters"
function OverlandFlowErosionAnswersParameters(dataset, config, indices)
usle_k = ncread(
dataset,
config,
"vertical.overland_flow_erosion.parameters.usle_k";
sel = indices,
defaults = 0.1,
type = Float,
)
usle_c = ncread(
dataset,
config,
"vertical.overland_flow_erosion.parameters.usle_c";
sel = indices,
defaults = 0.01,
type = Float,
)
answers_k = ncread(
dataset,
config,
"vertical.overland_flow_erosion.parameters.answers_k";
sel = indices,
defaults = 0.9,
type = Float,
)
lens = lens_input_parameter("soil_erosion__usle_k_factor")
usle_k = ncread(dataset, config, lens; sel = indices, defaults = 0.1, type = Float)
lens = lens_input_parameter("soil_erosion__usle_c_factor")
usle_c = ncread(dataset, config, lens; sel = indices, defaults = 0.01, type = Float)
lens = lens_input_parameter("soil_erosion__answers_overland_flow_factor")
answers_k = ncread(dataset, config, lens; sel = indices, defaults = 0.9, type = Float)

answers_parameters = OverlandFlowErosionAnswersParameters(;
usle_k = usle_k,
usle_c = usle_c,
Expand Down
76 changes: 20 additions & 56 deletions src/sediment/erosion/rainfall_erosion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -51,46 +51,22 @@ end

"Initialize EUROSEM rainfall erosion model parameters"
function RainfallErosionEurosemParameters(dataset, config, indices)
soil_detachability = ncread(
dataset,
config,
"vertical.rainfall_erosion.parameters.soil_detachability";
sel = indices,
defaults = 0.6,
type = Float,
)
eurosem_exponent = ncread(
dataset,
config,
"vertical.rainfall_erosion.parameters.eurosem_exponent";
sel = indices,
defaults = 2.0,
type = Float,
)
canopyheight = ncread(
dataset,
config,
"vertical.rainfall_erosion.parameters.canopyheight";
sel = indices,
defaults = 0.5,
type = Float,
)
canopygapfraction = ncread(
dataset,
config,
"vertical.rainfall_erosion.parameters.canopygapfraction";
sel = indices,
defaults = 0.1,
type = Float,
)
soilcover_fraction = ncread(
dataset,
config,
"vertical.rainfall_erosion.parameters.pathfrac";
sel = indices,
defaults = 0.01,
type = Float,
)
lens = lens_input_parameter("soil_erosion__rainfall_soil_detachability_factor")
soil_detachability =
ncread(dataset, config, lens; sel = indices, defaults = 0.6, type = Float)
lens = lens_input_parameter("soil_erosion__eurosem_exponent")
eurosem_exponent =
ncread(dataset, config, lens; sel = indices, defaults = 2.0, type = Float)
lens = lens_input_parameter("vegetation_canopy__height")
canopyheight =
ncread(dataset, config, lens; sel = indices, defaults = 0.5, type = Float)
lens = lens_input_parameter("vegetation_canopy__gap_fraction")
canopygapfraction =
ncread(dataset, config, lens; sel = indices, defaults = 0.1, type = Float)
lens = lens_input_parameter("soil~compacted__area_fraction")
soilcover_fraction =
ncread(dataset, config, lens; sel = indices, defaults = 0.01, type = Float)

eurosem_parameters = RainfallErosionEurosemParameters(;
soil_detachability = soil_detachability,
eurosem_exponent = eurosem_exponent,
Expand Down Expand Up @@ -184,22 +160,10 @@ end

"Initialize ANSWERS rainfall erosion model parameters"
function RainfallErosionAnswersParameters(dataset, config, indices)
usle_k = ncread(
dataset,
config,
"vertical.rainfall_erosion.parameters.usle_k";
sel = indices,
defaults = 0.1,
type = Float,
)
usle_c = ncread(
dataset,
config,
"vertical.rainfall_erosion.parameters.usle_c";
sel = indices,
defaults = 0.01,
type = Float,
)
lens = lens_input_parameter("soil_erosion__usle_k_factor")
usle_k = ncread(dataset, config, lens; sel = indices, defaults = 0.1, type = Float)
lens = lens_input_parameter("soil_erosion__usle_c_factor")
usle_c = ncread(dataset, config, lens; sel = indices, defaults = 0.01, type = Float)
answers_parameters =
RainfallErosionAnswersParameters(; usle_k = usle_k, usle_c = usle_c)
return answers_parameters
Expand Down
10 changes: 2 additions & 8 deletions src/sediment/erosion/river_erosion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -43,14 +43,8 @@ end

"Initialize Julian and Torres river erosion parameters"
function RiverErosionParameters(dataset, config, indices)
d50 = ncread(
dataset,
config,
"lateral.river.potential_erosion.parameters.d50";
sel = indices,
defaults = 0.1,
type = Float,
)
lens = lens_input_parameter("river_bottom-and-bank_sediment__d50_diameter")
d50 = ncread(dataset, config, lens; sel = indices, defaults = 0.1, type = Float)
river_parameters = RiverErosionParameters(; d50 = d50)

return river_parameters
Expand Down
Loading

0 comments on commit 8b38fb2

Please sign in to comment.