From 16e5259384e8871a9023effd26e0c2b9c9741b13 Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Tue, 26 Nov 2024 15:54:15 +0100
Subject: [PATCH] Rename variables `nc` and `inds`

---
 src/demand/water_demand.jl             |  90 +++++++--------
 src/flow.jl                            |   6 +-
 src/glacier/glacier.jl                 |  30 ++---
 src/groundwater/aquifer.jl             |  35 ++++--
 src/groundwater/boundary_conditions.jl |  30 +++--
 src/horizontal_process.jl              |   8 +-
 src/parameters.jl                      |  32 +++---
 src/reservoir_lake.jl                  |  22 ++--
 src/sbm.jl                             |  21 ++--
 src/snow/snow.jl                       |  28 ++---
 src/soil/soil.jl                       | 149 ++++++++++++++-----------
 src/surfacewater/runoff.jl             |  10 +-
 src/vegetation/canopy.jl               |   8 +-
 13 files changed, 253 insertions(+), 216 deletions(-)

diff --git a/src/demand/water_demand.jl b/src/demand/water_demand.jl
index a9c09f451..2dd5dc8b8 100644
--- a/src/demand/water_demand.jl
+++ b/src/demand/water_demand.jl
@@ -31,26 +31,26 @@ get_demand_gross(model::NonIrrigationDemand) = model.demand.demand_gross
 get_demand_gross(model::NoNonIrrigationDemand) = 0.0
 
 "Initialize non-irrigation water demand model for a water use `sector`"
-function NonIrrigationDemand(nc, config, inds, dt, sector)
+function NonIrrigationDemand(dataset, config, indices, dt, sector)
     demand_gross =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.demand.$(sector).demand_gross";
-            sel = inds,
+            sel = indices,
             defaults = 0.0,
             type = Float,
         ) .* (dt / basetimestep)
     demand_net =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.demand.$(sector).demand_net";
-            sel = inds,
+            sel = indices,
             defaults = 0.0,
             type = Float,
         ) .* (dt / basetimestep)
-    n = length(inds)
+    n = length(indices)
     returnflow_f = return_flow_fraction.(demand_gross, demand_net)
 
     demand = PrescibedDemand{Float}(; demand_gross, demand_net)
@@ -83,39 +83,39 @@ end
 end
 
 "Initialize non-paddy irrigation model"
-function NonPaddy(nc, config, inds, dt)
+function NonPaddy(dataset, config, indices, dt)
     efficiency = ncread(
-        nc,
+        dataset,
         config,
         "vertical.demand.nonpaddy.parameters.irrigation_efficiency";
-        sel = inds,
+        sel = indices,
         defaults = 1.0,
         type = Float,
     )
     areas = ncread(
-        nc,
+        dataset,
         config,
         "vertical.demand.nonpaddy.parameters.irrigation_areas";
-        sel = inds,
+        sel = indices,
         defaults = 1,
         optional = false,
         type = Int,
     )
     irrigation_trigger = ncread(
-        nc,
+        dataset,
         config,
         "vertical.demand.nonpaddy.parameters.irrigation_trigger";
-        sel = inds,
+        sel = indices,
         defaults = 1,
         optional = false,
         type = Bool,
     )
     max_irri_rate =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.demand.nonpaddy.parameters.maximum_irrigation_rate";
-            sel = inds,
+            sel = indices,
             defaults = 25.0,
             type = Float,
         ) .* (dt / basetimestep)
@@ -126,7 +126,7 @@ function NonPaddy(nc, config, inds, dt)
         irrigation_areas = areas,
         irrigation_trigger,
     )
-    vars = NonPaddyVariables{Float}(; demand_gross = fill(mv, length(inds)))
+    vars = NonPaddyVariables{Float}(; demand_gross = fill(mv, length(indices)))
 
     nonpaddy = NonPaddy{Float}(; variables = vars, parameters = params)
 
@@ -229,65 +229,65 @@ end
 end
 
 "Initialize paddy irrigation model"
-function Paddy(nc, config, inds, dt)
+function Paddy(dataset, config, indices, dt)
     h_min = ncread(
-        nc,
+        dataset,
         config,
         "vertical.demand.paddy.parameters.h_min";
-        sel = inds,
+        sel = indices,
         defaults = 20.0,
         type = Float,
     )
     h_opt = ncread(
-        nc,
+        dataset,
         config,
         "vertical.demand.paddy.parameters.h_opt";
-        sel = inds,
+        sel = indices,
         defaults = 50.0,
         type = Float,
     )
     h_max = ncread(
-        nc,
+        dataset,
         config,
         "vertical.demand.paddy.parameters.h_max";
-        sel = inds,
+        sel = indices,
         defaults = 80.0,
         type = Float,
     )
     efficiency = ncread(
-        nc,
+        dataset,
         config,
         "vertical.demand.paddy.parameters.irrigation_efficiency";
-        sel = inds,
+        sel = indices,
         defaults = 1.0,
         type = Float,
     )
     areas = ncread(
-        nc,
+        dataset,
         config,
         "vertical.demand.paddy.parameters.irrigation_areas";
-        sel = inds,
+        sel = indices,
         optional = false,
         type = Bool,
     )
     irrigation_trigger = ncread(
-        nc,
+        dataset,
         config,
         "vertical.demand.paddy.parameters.irrigation_trigger";
-        sel = inds,
+        sel = indices,
         optional = false,
         type = Bool,
     )
     max_irri_rate =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.demand.paddy.parameters.maximum_irrigation_rate";
-            sel = inds,
+            sel = indices,
             defaults = 25.0,
             type = Float,
         ) .* (dt / basetimestep)
-    n = length(inds)
+    n = length(indices)
     params = PaddyParameters{Float}(;
         irrigation_efficiency = efficiency,
         maximum_irrigation_rate = max_irri_rate,
@@ -437,34 +437,34 @@ end
 end
 
 "Initialize water demand model"
-function Demand(nc, config, inds, dt)
+function Demand(dataset, config, indices, dt)
     domestic = if get(config.model.water_demand, "domestic", false)
-        NonIrrigationDemand(nc, config, inds, dt, "domestic")
+        NonIrrigationDemand(dataset, config, indices, dt, "domestic")
     else
         NoNonIrrigationDemand()
     end
     industry = if get(config.model.water_demand, "industry", false)
-        NonIrrigationDemand(nc, config, inds, dt, "industry")
+        NonIrrigationDemand(dataset, config, indices, dt, "industry")
     else
         NoNonIrrigationDemand()
     end
     livestock = if get(config.model.water_demand, "livestock", false)
-        NonIrrigationDemand(nc, config, inds, dt, "livestock")
+        NonIrrigationDemand(dataset, config, indices, dt, "livestock")
     else
         NoNonIrrigationDemand()
     end
     paddy = if get(config.model.water_demand, "paddy", false)
-        Paddy(nc, config, inds, dt)
+        Paddy(dataset, config, indices, dt)
     else
         NoIrrigationPaddy{Float}()
     end
     nonpaddy = if get(config.model.water_demand, "nonpaddy", false)
-        NonPaddy(nc, config, inds, dt)
+        NonPaddy(dataset, config, indices, dt)
     else
         NoIrrigationNonPaddy{Float}()
     end
 
-    n = length(inds)
+    n = length(indices)
     vars = DemandVariables(Float, n)
     demand = Demand(; domestic, industry, livestock, paddy, nonpaddy, variables = vars)
     return demand
@@ -544,25 +544,25 @@ end
 end
 
 "Initialize water allocation for the land domain"
-function AllocationLand(nc, config, inds)
+function AllocationLand(dataset, config, indices)
     frac_sw_used = ncread(
-        nc,
+        dataset,
         config,
         "vertical.allocation.parameters.frac_sw_used";
-        sel = inds,
+        sel = indices,
         defaults = 1,
         type = Float,
     )
     areas = ncread(
-        nc,
+        dataset,
         config,
         "vertical.allocation.parameters.areas";
-        sel = inds,
+        sel = indices,
         defaults = 1,
         type = Int,
     )
 
-    n = length(inds)
+    n = length(indices)
 
     params = AllocationLandParameters(; areas = areas, frac_sw_used = frac_sw_used)
     vars = AllocationLandVariables(Float, n)
diff --git a/src/flow.jl b/src/flow.jl
index 4c2173c54..56e9a2ee7 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -885,10 +885,10 @@ end
     error::Vector{T} | "m3"                             # error volume    
 end
 
-function ShallowWaterRiverVariables(nc, config, inds, n_edges, inds_pit)
+function ShallowWaterRiverVariables(dataset, config, indices, n_edges, inds_pit)
     floodplain_1d = get(config.model, "floodplain_1d", false)::Bool
     riverdepth_bc = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.riverdepth_bc";
         sel = inds_pit,
@@ -896,7 +896,7 @@ function ShallowWaterRiverVariables(nc, config, inds, n_edges, inds_pit)
         type = Float,
     )
 
-    n = length(inds)
+    n = length(indices)
     # set river depth h to zero (including reservoir and lake locations)
     h = zeros(n)
     q_av = zeros(n_edges)
diff --git a/src/glacier/glacier.jl b/src/glacier/glacier.jl
index 8f5be9553..d61087be1 100644
--- a/src/glacier/glacier.jl
+++ b/src/glacier/glacier.jl
@@ -9,12 +9,12 @@ abstract type AbstractGlacierModel{T} end
 end
 
 "Initialize glacier model variables"
-function GlacierVariables(nc, config, inds)
+function GlacierVariables(dataset, config, indices)
     glacier_store = ncread(
-        nc,
+        dataset,
         config,
         "vertical.glacier.variables.glacier_store";
-        sel = inds,
+        sel = indices,
         defaults = 5500.0,
         type = Float,
         fill = 0.0,
@@ -54,41 +54,41 @@ end
 struct NoGlacierModel{T} <: AbstractGlacierModel{T} end
 
 "Initialize glacier HBV model parameters"
-function GlacierHbvParameters(nc, config, inds, dt)
+function GlacierHbvParameters(dataset, config, indices, dt)
     g_tt = ncread(
-        nc,
+        dataset,
         config,
         "vertical.glacier.parameters.g_tt";
-        sel = inds,
+        sel = indices,
         defaults = 0.0,
         type = Float,
         fill = 0.0,
     )
     g_cfmax =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.glacier.parameters.g_cfmax";
-            sel = inds,
+            sel = indices,
             defaults = 3.0,
             type = Float,
             fill = 0.0,
         ) .* (dt / basetimestep)
     g_sifrac =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.glacier.parameters.g_sifrac";
-            sel = inds,
+            sel = indices,
             defaults = 0.001,
             type = Float,
             fill = 0.0,
         ) .* (dt / basetimestep)
     glacier_frac = ncread(
-        nc,
+        dataset,
         config,
         "vertical.glacier.parameters.glacier_frac";
-        sel = inds,
+        sel = indices,
         defaults = 0.0,
         type = Float,
         fill = 0.0,
@@ -105,9 +105,9 @@ function GlacierHbvParameters(nc, config, inds, dt)
 end
 
 "Initialize glacier HBV model"
-function GlacierHbvModel(nc, config, inds, dt, bc)
-    params = GlacierHbvParameters(nc, config, inds, dt)
-    vars = GlacierVariables(nc, config, inds)
+function GlacierHbvModel(dataset, config, indices, dt, bc)
+    params = GlacierHbvParameters(dataset, config, indices, dt)
+    vars = GlacierVariables(dataset, config, indices)
     model =
         GlacierHbvModel(; boundary_conditions = bc, parameters = params, variables = vars)
     return model
diff --git a/src/groundwater/aquifer.jl b/src/groundwater/aquifer.jl
index 32b901100..7b6689f4c 100644
--- a/src/groundwater/aquifer.jl
+++ b/src/groundwater/aquifer.jl
@@ -126,15 +126,26 @@ instead. Specific yield will vary roughly between 0.05 (clay) and 0.45 (peat)
     # conductivity_profile is set to "exponential")
 end
 
-function UnconfinedAquiferParameters(nc, config, inds, top, bottom, area)
-    k = ncread(nc, config, "lateral.subsurface.conductivity"; sel = inds, type = Float)
-    specific_yield =
-        ncread(nc, config, "lateral.subsurface.specific_yield"; sel = inds, type = Float)
+function UnconfinedAquiferParameters(dataset, config, indices, top, bottom, area)
+    k = ncread(
+        dataset,
+        config,
+        "lateral.subsurface.conductivity";
+        sel = indices,
+        type = Float,
+    )
+    specific_yield = ncread(
+        dataset,
+        config,
+        "lateral.subsurface.specific_yield";
+        sel = indices,
+        type = Float,
+    )
     f = ncread(
-        nc,
+        dataset,
         config,
         "lateral.subsurface.gwf_f";
-        sel = inds,
+        sel = indices,
         type = Float,
         defaults = 3.0,
     )
@@ -155,8 +166,8 @@ end
     variables::UnconfinedAquiferVariables{T}
 end
 
-function UnconfinedAquifer(nc, config, inds, top, bottom, area, conductance, head)
-    parameters = UnconfinedAquiferParameters(nc, config, inds, top, bottom, area)
+function UnconfinedAquifer(dataset, config, indices, top, bottom, area, conductance, head)
+    parameters = UnconfinedAquiferParameters(dataset, config, indices, top, bottom, area)
 
     volume = @. (min(top, head) - bottom) * area * parameters.specific_yield
     variables = UnconfinedAquiferVariables{Float}(head, conductance, volume)
@@ -368,16 +379,16 @@ end
     index::Vector{Int} | "-"
 end
 
-function ConstantHead(nc, config, inds)
+function ConstantHead(dataset, config, indices)
     constanthead = ncread(
-        nc,
+        dataset,
         config,
         "lateral.subsurface.constant_head";
-        sel = inds,
+        sel = indices,
         type = Float,
         fill = mv,
     )
-    n = length(inds)
+    n = length(indices)
     index_constanthead = filter(i -> !isequal(constanthead[i], mv), 1:n)
     head = constanthead[index_constanthead]
     variables = ConstantHeadVariables{Float}(head)
diff --git a/src/groundwater/boundary_conditions.jl b/src/groundwater/boundary_conditions.jl
index 5e3a45740..52927733c 100644
--- a/src/groundwater/boundary_conditions.jl
+++ b/src/groundwater/boundary_conditions.jl
@@ -33,26 +33,32 @@ end
     index::Vector{Int} | "-"
 end
 
-function River(nc, config, inds, index)
+function River(dataset, config, indices, index)
     infiltration_conductance = ncread(
-        nc,
+        dataset,
         config,
         "lateral.subsurface.infiltration_conductance";
-        sel = inds,
+        sel = indices,
         type = Float,
     )
     exfiltration_conductance = ncread(
-        nc,
+        dataset,
         config,
         "lateral.subsurface.exfiltration_conductance";
-        sel = inds,
+        sel = indices,
+        type = Float,
+    )
+    bottom = ncread(
+        dataset,
+        config,
+        "lateral.subsurface.river_bottom";
+        sel = indices,
         type = Float,
     )
-    bottom = ncread(nc, config, "lateral.subsurface.river_bottom"; sel = inds, type = Float)
 
     parameters =
         RiverParameters{Float}(infiltration_conductance, exfiltration_conductance, bottom)
-    n = length(inds)
+    n = length(indices)
     variables = RiverVariables(n)
     river = River(parameters, variables, index)
     return river
@@ -90,20 +96,20 @@ end
     index::Vector{Int} | "-"
 end
 
-function Drainage(nc, config, inds, index)
+function Drainage(dataset, config, indices, index)
     drain_elevation = ncread(
-        nc,
+        dataset,
         config,
         "lateral.subsurface.drain_elevation";
-        sel = inds,
+        sel = indices,
         type = Float,
         fill = mv,
     )
     drain_conductance = ncread(
-        nc,
+        dataset,
         config,
         "lateral.subsurface.drain_conductance";
-        sel = inds,
+        sel = indices,
         type = Float,
         fill = mv,
     )
diff --git a/src/horizontal_process.jl b/src/horizontal_process.jl
index fe9c57252..e51c5f9af 100644
--- a/src/horizontal_process.jl
+++ b/src/horizontal_process.jl
@@ -1,17 +1,17 @@
 "Convert a gridded drainage direction to a directed graph"
-function flowgraph(ldd::AbstractVector, inds::AbstractVector, pcr_dir::AbstractVector)
+function flowgraph(ldd::AbstractVector, indices::AbstractVector, pcr_dir::AbstractVector)
     # prepare a directed graph to be filled
-    n = length(inds)
+    n = length(indices)
     graph = DiGraph(n)
 
     # loop over ldd, adding the edge to the downstream node
-    for (from_node, from_index) in enumerate(inds)
+    for (from_node, from_index) in enumerate(indices)
         ldd_val = ldd[from_node]
         # skip pits to prevent cycles
         ldd_val == 5 && continue
         to_index = from_index + pcr_dir[ldd_val]
         # find the node id of the downstream cell
-        to_node = searchsortedfirst(inds, to_index)
+        to_node = searchsortedfirst(indices, to_index)
         add_edge!(graph, from_node, to_node)
     end
     if is_cyclic(graph)
diff --git a/src/parameters.jl b/src/parameters.jl
index f2d2d582a..91dc3e662 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -19,47 +19,47 @@
 end
 
 "Initialize (shared) vegetation parameters"
-function VegetationParameters(nc, config, inds)
-    n = length(inds)
+function VegetationParameters(dataset, config, indices)
+    n = length(indices)
     rootingdepth = ncread(
-        nc,
+        dataset,
         config,
         "vertical.vegetation_parameter_set.rootingdepth";
-        sel = inds,
+        sel = indices,
         defaults = 750.0,
         type = Float,
     )
     kc = ncread(
-        nc,
+        dataset,
         config,
         "vertical.vegetation_parameter_set.kc";
-        sel = inds,
+        sel = indices,
         defaults = 1.0,
         type = Float,
     )
     if haskey(config.input.vertical.vegetation_parameter_set, "leaf_area_index")
         storage_specific_leaf = ncread(
-            nc,
+            dataset,
             config,
             "vertical.vegetation_parameter_set.storage_specific_leaf";
             optional = false,
-            sel = inds,
+            sel = indices,
             type = Float,
         )
         storage_wood = ncread(
-            nc,
+            dataset,
             config,
             "vertical.vegetation_parameter_set.storage_wood";
             optional = false,
-            sel = inds,
+            sel = indices,
             type = Float,
         )
         kext = ncread(
-            nc,
+            dataset,
             config,
             "vertical.vegetation_parameter_set.kext";
             optional = false,
-            sel = inds,
+            sel = indices,
             type = Float,
         )
         vegetation_parameter_set = VegetationParameters(;
@@ -74,18 +74,18 @@ function VegetationParameters(nc, config, inds)
         )
     else
         canopygapfraction = ncread(
-            nc,
+            dataset,
             config,
             "vertical.vegetation_parameter_set.canopygapfraction";
-            sel = inds,
+            sel = indices,
             defaults = 0.1,
             type = Float,
         )
         cmax = ncread(
-            nc,
+            dataset,
             config,
             "vertical.vegetation_parameter_set.cmax";
-            sel = inds,
+            sel = indices,
             defaults = 1.0,
             type = Float,
         )
diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index 7fb339030..a2ecb9be8 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -255,12 +255,12 @@ end
     hq::Vector{Union{HQ, Missing}}                  # data for rating curve    
 end
 
-function LakeParameters(config, nc, inds_riv, nriv, pits, dt)
+function LakeParameters(config, dataset, inds_riv, nriv, pits, dt)
     # read only lake data if lakes true
     # allow lakes only in river cells
     # note that these locations are only the lake outlet pixels
     lakelocs_2d = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.lake.locs";
         optional = false,
@@ -271,7 +271,7 @@ function LakeParameters(config, nc, inds_riv, nriv, pits, dt)
 
     # this holds the same ids as lakelocs, but covers the entire lake
     lakecoverage_2d = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.lake.areas";
         optional = false,
@@ -304,7 +304,7 @@ function LakeParameters(config, nc, inds_riv, nriv, pits, dt)
     end
 
     lakearea = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.lake.area";
         optional = false,
@@ -313,7 +313,7 @@ function LakeParameters(config, nc, inds_riv, nriv, pits, dt)
         fill = 0,
     )
     lake_b = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.lake.b";
         optional = false,
@@ -322,7 +322,7 @@ function LakeParameters(config, nc, inds_riv, nriv, pits, dt)
         fill = 0,
     )
     lake_e = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.lake.e";
         optional = false,
@@ -331,7 +331,7 @@ function LakeParameters(config, nc, inds_riv, nriv, pits, dt)
         fill = 0,
     )
     lake_threshold = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.lake.threshold";
         optional = false,
@@ -340,7 +340,7 @@ function LakeParameters(config, nc, inds_riv, nriv, pits, dt)
         fill = 0,
     )
     linked_lakelocs = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.lake.linkedlakelocs";
         sel = inds_lake,
@@ -349,7 +349,7 @@ function LakeParameters(config, nc, inds_riv, nriv, pits, dt)
         fill = 0,
     )
     lake_storfunc = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.lake.storfunc";
         optional = false,
@@ -358,7 +358,7 @@ function LakeParameters(config, nc, inds_riv, nriv, pits, dt)
         fill = 0,
     )
     lake_outflowfunc = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.lake.outflowfunc";
         optional = false,
@@ -367,7 +367,7 @@ function LakeParameters(config, nc, inds_riv, nriv, pits, dt)
         fill = 0,
     )
     lake_waterlevel = ncread(
-        nc,
+        dataset,
         config,
         "lateral.river.lake.waterlevel";
         optional = false,
diff --git a/src/sbm.jl b/src/sbm.jl
index ffed27f18..0e9f286c7 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -12,29 +12,29 @@
 end
 
 "Initialize land hydrology model with SBM soil model"
-function LandHydrologySBM(nc, config, riverfrac, inds)
+function LandHydrologySBM(dataset, config, riverfrac, indices)
     dt = Second(config.timestepsecs)
-    n = length(inds)
+    n = length(indices)
 
     atmospheric_forcing = AtmosphericForcing(n)
-    vegetation_parameter_set = VegetationParameters(nc, config, inds)
+    vegetation_parameter_set = VegetationParameters(dataset, config, indices)
     if dt >= Hour(23)
         interception_model =
-            GashInterceptionModel(nc, config, inds, vegetation_parameter_set)
+            GashInterceptionModel(dataset, config, indices, vegetation_parameter_set)
     else
         interception_model = RutterInterceptionModel(vegetation_parameter_set, n)
     end
 
     modelsnow = get(config.model, "snow", false)::Bool
     if modelsnow
-        snow_model = SnowHbvModel(nc, config, inds, dt)
+        snow_model = SnowHbvModel(dataset, config, indices, dt)
     else
         snow_model = NoSnowModel{Float}()
     end
     modelglacier = get(config.model, "glacier", false)::Bool
     if modelsnow && modelglacier
         glacier_bc = SnowStateBC{Float}(; snow_storage = snow_model.variables.snow_storage)
-        glacier_model = GlacierHbvModel(nc, config, inds, dt, glacier_bc)
+        glacier_model = GlacierHbvModel(dataset, config, indices, dt, glacier_bc)
     elseif modelsnow == false && modelglacier == true
         @warn string(
             "Glacier processes can be modelled when snow modelling is enabled. To include ",
@@ -44,9 +44,9 @@ function LandHydrologySBM(nc, config, riverfrac, inds)
     else
         glacier_model = NoGlacierModel{Float}()
     end
-    runoff_model = OpenWaterRunoff(nc, config, inds, riverfrac)
+    runoff_model = OpenWaterRunoff(dataset, config, indices, riverfrac)
 
-    soil_model = SbmSoilModel(nc, config, vegetation_parameter_set, inds, dt)
+    soil_model = SbmSoilModel(dataset, config, vegetation_parameter_set, indices, dt)
     @. vegetation_parameter_set.rootingdepth = min(
         soil_model.parameters.soilthickness * 0.99,
         vegetation_parameter_set.rootingdepth,
@@ -54,8 +54,9 @@ function LandHydrologySBM(nc, config, riverfrac, inds)
 
     do_water_demand = haskey(config.model, "water_demand")
     allocation =
-        do_water_demand ? AllocationLand(nc, config, inds) : NoAllocationLand{Float}()
-    demand = do_water_demand ? Demand(nc, config, inds, dt) : NoDemand{Float}()
+        do_water_demand ? AllocationLand(dataset, config, indices) :
+        NoAllocationLand{Float}()
+    demand = do_water_demand ? Demand(dataset, config, indices, dt) : NoDemand{Float}()
 
     args = (demand, allocation)
     land_hydrology_model = LandHydrologySBM{Float, typeof.(args)...}(;
diff --git a/src/snow/snow.jl b/src/snow/snow.jl
index a60521b13..c67c54ea5 100644
--- a/src/snow/snow.jl
+++ b/src/snow/snow.jl
@@ -68,45 +68,45 @@ end
 struct NoSnowModel{T} <: AbstractSnowModel{T} end
 
 "Initialize snow HBV model parameters"
-function SnowHbvParameters(nc, config, inds, dt)
+function SnowHbvParameters(dataset, config, indices, dt)
     cfmax =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.snow.parameters.cfmax";
-            sel = inds,
+            sel = indices,
             defaults = 3.75653,
             type = Float,
         ) .* (dt / basetimestep)
     tt = ncread(
-        nc,
+        dataset,
         config,
         "vertical.snow.parameters.tt";
-        sel = inds,
+        sel = indices,
         defaults = 0.0,
         type = Float,
     )
     tti = ncread(
-        nc,
+        dataset,
         config,
         "vertical.snow.parameters.tti";
-        sel = inds,
+        sel = indices,
         defaults = 1.0,
         type = Float,
     )
     ttm = ncread(
-        nc,
+        dataset,
         config,
         "vertical.snow.parameters.ttm";
-        sel = inds,
+        sel = indices,
         defaults = 0.0,
         type = Float,
     )
     whc = ncread(
-        nc,
+        dataset,
         config,
         "vertical.snow.parameters.whc";
-        sel = inds,
+        sel = indices,
         defaults = 0.1,
         type = Float,
     )
@@ -116,9 +116,9 @@ function SnowHbvParameters(nc, config, inds, dt)
 end
 
 "Initialize snow HBV model"
-function SnowHbvModel(nc, config, inds, dt)
-    n = length(inds)
-    params = SnowHbvParameters(nc, config, inds, dt)
+function SnowHbvModel(dataset, config, indices, dt)
+    n = length(indices)
+    params = SnowHbvParameters(dataset, config, indices, dt)
     vars = SnowVariables(Float, n)
     bc = SnowBC(Float, n)
     model = SnowHbvModel(; boundary_conditions = bc, parameters = params, variables = vars)
diff --git a/src/soil/soil.jl b/src/soil/soil.jl
index 157869c2c..bdd649e9d 100644
--- a/src/soil/soil.jl
+++ b/src/soil/soil.jl
@@ -200,18 +200,28 @@ end
 end
 
 "Initialize SBM soil model hydraulic conductivity depth profile"
-function sbm_kv_profiles(nc, config, inds, kv_0, f, maxlayers, nlayers, sumlayers, dt)
+function sbm_kv_profiles(
+    dataset,
+    config,
+    indices,
+    kv_0,
+    f,
+    maxlayers,
+    nlayers,
+    sumlayers,
+    dt,
+)
     kv_profile_type = get(config.input.vertical, "ksat_profile", "exponential")::String
-    n = length(inds)
+    n = length(indices)
     if kv_profile_type == "exponential"
         kv_profile = KvExponential(kv_0, f)
     elseif kv_profile_type == "exponential_constant"
         z_exp = ncread(
-            nc,
+            dataset,
             config,
             "vertical.soil.parameters.z_exp";
             optional = false,
-            sel = inds,
+            sel = indices,
             type = Float,
         )
         exp_profile = KvExponential(kv_0, f)
@@ -219,10 +229,10 @@ function sbm_kv_profiles(nc, config, inds, kv_0, f, maxlayers, nlayers, sumlayer
     elseif kv_profile_type == "layered" || kv_profile_type == "layered_exponential"
         kv =
             ncread(
-                nc,
+                dataset,
                 config,
                 "vertical.soil.parameters.kv";
-                sel = inds,
+                sel = indices,
                 defaults = 1000.0,
                 type = Float,
                 dimname = :layer,
@@ -236,11 +246,11 @@ function sbm_kv_profiles(nc, config, inds, kv_0, f, maxlayers, nlayers, sumlayer
             kv_profile = KvLayered(svectorscopy(kv, Val{maxlayers}()))
         else
             z_layered = ncread(
-                nc,
+                dataset,
                 config,
                 "vertical.soil.parameters.z_layered";
                 optional = false,
-                sel = inds,
+                sel = indices,
                 type = Float,
             )
             nlayers_kv = fill(0, n)
@@ -331,7 +341,7 @@ end
 end
 
 "Initialize SBM soil model parameters"
-function SbmSoilParameters(nc, config, vegetation_parameter_set, inds, dt)
+function SbmSoilParameters(dataset, config, vegetation_parameter_set, indices, dt)
     config_thicknesslayers = get(config.model, "thicknesslayers", Float[])
     if length(config_thicknesslayers) > 0
         thicknesslayers = SVector(Tuple(push!(Float.(config_thicknesslayers), mv)))
@@ -344,152 +354,152 @@ function SbmSoilParameters(nc, config, vegetation_parameter_set, inds, dt)
     end
     w_soil =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.soil.parameters.w_soil";
-            sel = inds,
+            sel = indices,
             defaults = 0.1125,
             type = Float,
         ) .* (dt / basetimestep)
     cf_soil = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.cf_soil";
-        sel = inds,
+        sel = indices,
         defaults = 0.038,
         type = Float,
     )
     # soil parameters
     theta_s = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.theta_s";
-        sel = inds,
+        sel = indices,
         defaults = 0.6,
         type = Float,
     )
     theta_r = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.theta_r";
-        sel = inds,
+        sel = indices,
         defaults = 0.01,
         type = Float,
     )
     kv_0 =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.soil.parameters.kv_0";
-            sel = inds,
+            sel = indices,
             defaults = 3000.0,
             type = Float,
         ) .* (dt / basetimestep)
     f = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.f";
-        sel = inds,
+        sel = indices,
         defaults = 0.001,
         type = Float,
     )
     hb = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.hb";
-        sel = inds,
+        sel = indices,
         defaults = -10.0,
         type = Float,
     )
     h1 = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.h1";
-        sel = inds,
+        sel = indices,
         defaults = 0.0,
         type = Float,
     )
     h2 = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.h2";
-        sel = inds,
+        sel = indices,
         defaults = -100.0,
         type = Float,
     )
     h3_high = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.h3_high";
-        sel = inds,
+        sel = indices,
         defaults = -400.0,
         type = Float,
     )
     h3_low = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.h3_low";
-        sel = inds,
+        sel = indices,
         defaults = -1000.0,
         type = Float,
     )
     h4 = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.h4";
-        sel = inds,
+        sel = indices,
         defaults = -15849.0,
         type = Float,
     )
     alpha_h1 = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.alpha_h1";
-        sel = inds,
+        sel = indices,
         defaults = 1.0,
         type = Float,
     )
     soilthickness = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.soilthickness";
-        sel = inds,
+        sel = indices,
         defaults = 2000.0,
         type = Float,
     )
     infiltcappath =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.soil.parameters.infiltcappath";
-            sel = inds,
+            sel = indices,
             defaults = 10.0,
             type = Float,
         ) .* (dt / basetimestep)
     infiltcapsoil =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.soil.parameters.infiltcapsoil";
-            sel = inds,
+            sel = indices,
             defaults = 100.0,
             type = Float,
         ) .* (dt / basetimestep)
     maxleakage =
         ncread(
-            nc,
+            dataset,
             config,
             "vertical.soil.parameters.maxleakage";
-            sel = inds,
+            sel = indices,
             defaults = 0.0,
             type = Float,
         ) .* (dt / basetimestep)
 
     c = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.c";
-        sel = inds,
+        sel = indices,
         defaults = 10.0,
         type = Float,
         dimname = :layer,
@@ -500,10 +510,10 @@ function SbmSoilParameters(nc, config, vegetation_parameter_set, inds, dt)
         error("$parname needs a layer dimension of size $maxlayers, but is $size1")
     end
     kvfrac = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.kvfrac";
-        sel = inds,
+        sel = indices,
         defaults = 1.0,
         type = Float,
         dimname = :layer,
@@ -515,36 +525,36 @@ function SbmSoilParameters(nc, config, vegetation_parameter_set, inds, dt)
     end
     # fraction compacted area
     pathfrac = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.pathfrac";
-        sel = inds,
+        sel = indices,
         defaults = 0.01,
         type = Float,
     )
 
     # vegetation parameters
     rootdistpar = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.rootdistpar";
-        sel = inds,
+        sel = indices,
         defaults = -500.0,
         type = Float,
     )
     cap_hmax = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.cap_hmax";
-        sel = inds,
+        sel = indices,
         defaults = 2000.0,
         type = Float,
     )
     cap_n = ncread(
-        nc,
+        dataset,
         config,
         "vertical.soil.parameters.cap_n";
-        sel = inds,
+        sel = indices,
         defaults = 2.0,
         type = Float,
     )
@@ -554,20 +564,20 @@ function SbmSoilParameters(nc, config, vegetation_parameter_set, inds, dt)
     nlayers = number_of_active_layers.(act_thickl)
 
     if length(config_thicknesslayers) > 0
-        # root fraction read from nc file, in case of multiple soil layers and TOML file
+        # root fraction read from dataset file, in case of multiple soil layers and TOML file
         # includes "vertical.rootfraction"
         if haskey(config.input.vertical.soil.parameters, "rootfraction")
             rootfraction = ncread(
-                nc,
+                dataset,
                 config,
                 "vertical.soil.parameters.rootfraction";
-                sel = inds,
+                sel = indices,
                 optional = false,
                 type = Float,
                 dimname = :layer,
             )
         else
-            n = length(inds)
+            n = length(indices)
             (; rootingdepth) = vegetation_parameter_set
             # default root fraction in case of multiple soil layers
             rootfraction = zeros(Float, maxlayers, n)
@@ -590,12 +600,21 @@ function SbmSoilParameters(nc, config, vegetation_parameter_set, inds, dt)
         rootfraction = ones(Float, maxlayers, n)
     end
 
-    kv_profile =
-        sbm_kv_profiles(nc, config, inds, kv_0, f, maxlayers, nlayers, sumlayers, dt)
+    kv_profile = sbm_kv_profiles(
+        dataset,
+        config,
+        indices,
+        kv_0,
+        f,
+        maxlayers,
+        nlayers,
+        sumlayers,
+        dt,
+    )
 
     soilwatercapacity = @. soilthickness * (theta_s - theta_r)
 
-    n = length(inds)
+    n = length(indices)
     sbm_params = SbmSoilParameters(;
         maxlayers,
         nlayers,
@@ -639,9 +658,9 @@ end
 end
 
 "Initialize SBM soil model"
-function SbmSoilModel(nc, config, vegetation_parameter_set, inds, dt)
-    n = length(inds)
-    params = SbmSoilParameters(nc, config, vegetation_parameter_set, inds, dt)
+function SbmSoilModel(dataset, config, vegetation_parameter_set, indices, dt)
+    n = length(indices)
+    params = SbmSoilParameters(dataset, config, vegetation_parameter_set, indices, dt)
     vars = SbmSoilVariables(n, params)
     bc = SbmSoilBC(n)
     model = SbmSoilModel(; boundary_conditions = bc, parameters = params, variables = vars)
diff --git a/src/surfacewater/runoff.jl b/src/surfacewater/runoff.jl
index c8dc5cb36..b53e55e5b 100644
--- a/src/surfacewater/runoff.jl
+++ b/src/surfacewater/runoff.jl
@@ -34,13 +34,13 @@ end
 end
 
 "Initialize open water runoff parameters"
-function OpenWaterRunoffParameters(nc, config, inds, riverfrac)
+function OpenWaterRunoffParameters(dataset, config, indices, riverfrac)
     # fraction open water
     waterfrac = ncread(
-        nc,
+        dataset,
         config,
         "vertical.runoff.parameters.waterfrac";
-        sel = inds,
+        sel = indices,
         defaults = 0.0,
         type = Float,
     )
@@ -73,11 +73,11 @@ end
 end
 
 "Initialize open water runoff model"
-function OpenWaterRunoff(nc, config, inds, riverfrac)
+function OpenWaterRunoff(dataset, config, indices, riverfrac)
     n = length(riverfrac)
     vars = OpenWaterRunoffVariables(Float, n)
     bc = OpenWaterRunoffBC(Float, n)
-    params = OpenWaterRunoffParameters(nc, config, inds, riverfrac)
+    params = OpenWaterRunoffParameters(dataset, config, indices, riverfrac)
     model =
         OpenWaterRunoff(; boundary_conditions = bc, parameters = params, variables = vars)
     return model
diff --git a/src/vegetation/canopy.jl b/src/vegetation/canopy.jl
index 98762ea50..2ab44e7ce 100644
--- a/src/vegetation/canopy.jl
+++ b/src/vegetation/canopy.jl
@@ -39,16 +39,16 @@ end
 end
 
 "Initialize Gash interception model"
-function GashInterceptionModel(nc, config, inds, vegetation_parameter_set)
+function GashInterceptionModel(dataset, config, indices, vegetation_parameter_set)
     e_r = ncread(
-        nc,
+        dataset,
         config,
         "vertical.interception.parameters.e_r";
-        sel = inds,
+        sel = indices,
         defaults = 0.1,
         type = Float,
     )
-    n = length(inds)
+    n = length(indices)
     params =
         GashParameters(; e_r = e_r, vegetation_parameter_set = vegetation_parameter_set)
     vars = InterceptionVariables(Float, n)