Repository: Deltares/Wflow.jl
Pull request #122 (2021-04-15 12:41:12)

support Float32 structs
This makes the code and test changes needed to be able to run with Float32 number types in the structs.

Right now we hardcode `const Float = Float64`, but if you change this to `Float32` the tests also pass. To make that work I needed to increase the test tolerance a bit, most often by comparing to a Float32, which causes the approx test to use a larger relative tolerance. This can be merged as is, since it already touches quite some code. But after this we want to be able to set this from the configuration file, as well as updating more code, especially Float64 literals, to prevent functions converting to Float64 in between.


Patch:
diff --git a/src/Wflow.jl b/src/Wflow.jl
index 3a553c86d..ff9eb68bb 100644
--- a/src/Wflow.jl
+++ b/src/Wflow.jl
@@ -18,6 +18,7 @@ using CFTime
 @metadata get_units "mm"
 
 const BMI = BasicModelInterface
+const Float = Float64
 
 mutable struct Clock{T}
     time::T
diff --git a/src/hbv_model.jl b/src/hbv_model.jl
index f0ffa3575..c454ce71a 100644
--- a/src/hbv_model.jl
+++ b/src/hbv_model.jl
@@ -43,35 +43,35 @@ function initialize_hbv_model(config::Config)
             param(config, "input.vertical.cfmax", nothing);
             sel = inds,
             defaults = 3.75653,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     tt = ncread(
         nc,
         param(config, "input.vertical.tt", nothing);
         sel = inds,
         defaults = -1.41934,
-        type = Float64,
+        type = Float,
     )
     tti = ncread(
         nc,
         param(config, "input.vertical.tti", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
     )
     ttm = ncread(
         nc,
         param(config, "input.vertical.ttm", nothing);
         sel = inds,
         defaults = -1.41934,
-        type = Float64,
+        type = Float,
     )
     whc = ncread(
         nc,
         param(config, "input.vertical.whc", nothing);
         sel = inds,
         defaults = 0.1,
-        type = Float64,
+        type = Float,
     )
     # glacier parameters
     g_tt = ncread(
@@ -79,7 +79,7 @@ function initialize_hbv_model(config::Config)
         param(config, "input.vertical.g_tt", nothing);
         sel = inds,
         defaults = 0.0,
-        type = Float64,
+        type = Float,
         fill = 0.0,
     )
     g_cfmax =
@@ -88,7 +88,7 @@ function initialize_hbv_model(config::Config)
             param(config, "input.vertical.g_cfmax", nothing);
             sel = inds,
             defaults = 3.0,
-            type = Float64,
+            type = Float,
             fill = 0.0,
         ) .* (Δt / basetimestep)
     g_sifrac = ncread(
@@ -96,7 +96,7 @@ function initialize_hbv_model(config::Config)
         param(config, "input.vertical.g_sifrac", nothing);
         sel = inds,
         defaults = 0.001,
-        type = Float64,
+        type = Float,
         fill = 0.0,
     )
     glacierfrac = ncread(
@@ -104,7 +104,7 @@ function initialize_hbv_model(config::Config)
         param(config, "input.vertical.glacierfrac", nothing);
         sel = inds,
         defaults = 0.0,
-        type = Float64,
+        type = Float,
         fill = 0.0,
     )
     glacierstore = ncread(
@@ -112,7 +112,7 @@ function initialize_hbv_model(config::Config)
         param(config, "input.vertical.glacierstore", nothing);
         sel = inds,
         defaults = 5500.0,
-        type = Float64,
+        type = Float,
         fill = 0.0,
     )
     fc = ncread(
@@ -120,21 +120,21 @@ function initialize_hbv_model(config::Config)
         param(config, "input.vertical.fc", nothing);
         sel = inds,
         defaults = 260.0,
-        type = Float64,
+        type = Float,
     )
     betaseepage = ncread(
         nc,
         param(config, "input.vertical.betaseepage", nothing);
         sel = inds,
         defaults = 1.8,
-        type = Float64,
+        type = Float,
     )
     lp = ncread(
         nc,
         param(config, "input.vertical.lp", nothing);
         sel = inds,
         defaults = 0.53,
-        type = Float64,
+        type = Float,
     )
     k4 =
         ncread(
@@ -142,7 +142,7 @@ function initialize_hbv_model(config::Config)
             param(config, "input.vertical.k4", nothing);
             sel = inds,
             defaults = 0.02307,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     kquickflow =
         ncread(
@@ -150,14 +150,14 @@ function initialize_hbv_model(config::Config)
             param(config, "input.vertical.kquickflow", nothing);
             sel = inds,
             defaults = 0.09880,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     suz = ncread(
         nc,
         param(config, "input.vertical.suz", nothing);
         sel = inds,
         defaults = 100.0,
-        type = Float64,
+        type = Float,
     )
     k0 =
         ncread(
@@ -165,7 +165,7 @@ function initialize_hbv_model(config::Config)
             param(config, "input.vertical.k0", nothing);
             sel = inds,
             defaults = 0.30,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     khq =
         ncread(
@@ -173,7 +173,7 @@ function initialize_hbv_model(config::Config)
             param(config, "input.vertical.khq", nothing);
             sel = inds,
             defaults = 0.09880,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     hq =
         ncread(
@@ -181,14 +181,14 @@ function initialize_hbv_model(config::Config)
             param(config, "input.vertical.hq", nothing);
             sel = inds,
             defaults = 3.27,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     alphanl = ncread(
         nc,
         param(config, "input.vertical.alphanl", nothing);
         sel = inds,
         defaults = 1.1,
-        type = Float64,
+        type = Float,
     )
     perc =
         ncread(
@@ -196,35 +196,35 @@ function initialize_hbv_model(config::Config)
             param(config, "input.vertical.perc", nothing);
             sel = inds,
             defaults = 0.4,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     cfr = ncread(
         nc,
         param(config, "input.vertical.cfr", nothing);
         sel = inds,
         defaults = 0.05,
-        type = Float64,
+        type = Float,
     )
     pcorr = ncread(
         nc,
         param(config, "input.vertical.pcorr", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
     )
     rfcf = ncread(
         nc,
         param(config, "input.vertical.rfcf", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
     )
     sfcf = ncread(
         nc,
         param(config, "input.vertical.sfcf", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
     )
     cflux =
         ncread(
@@ -232,35 +232,35 @@ function initialize_hbv_model(config::Config)
             param(config, "input.vertical.cflux", nothing);
             sel = inds,
             defaults = 2.0,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     icf = ncread(
         nc,
         param(config, "input.vertical.icf", nothing);
         sel = inds,
         defaults = 2.0,
-        type = Float64,
+        type = Float,
     )
     cevpf = ncread(
         nc,
         param(config, "input.vertical.cevpf", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
     )
     epf = ncread(
         nc,
         param(config, "input.vertical.epf", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
     )
     ecorr = ncread(
         nc,
         param(config, "input.vertical.ecorr", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
     )
 
     # read x, y coordinates and calculate cell length [m]
@@ -283,10 +283,10 @@ function initialize_hbv_model(config::Config)
     threshold = fc .* lp
 
     altitude =
-        ncread(nc, param(config, "input.vertical.altitude"); sel = inds, type = Float64)
+        ncread(nc, param(config, "input.vertical.altitude"); sel = inds, type = Float)
 
-    hbv = HBV{Float64}(
-        Δt = tosecond(Δt),
+    hbv = HBV{Float}(
+        Δt = Float(tosecond(Δt)),
         n = n,
         fc = fc,
         betaseepage = betaseepage,
@@ -322,9 +322,9 @@ function initialize_hbv_model(config::Config)
         glacierstore = glacierstore,
         glacierfrac = glacierfrac,
         # default (cold) states:
-        interceptionstorage = fill(0.0, n),
-        snow = fill(0.0, n),
-        snowwater = fill(0.0, n),
+        interceptionstorage = zeros(Float, n),
+        snow = zeros(Float, n),
+        snowwater = zeros(Float, n),
         soilmoisture = copy(fc),
         upperzonestorage = 0.2 .* fc,
         lowerzonestorage = 1.0 ./ (3.0 .* k4),
@@ -352,10 +352,10 @@ function initialize_hbv_model(config::Config)
     river_2d = ncread(nc, param(config, "input.river_location"); type = Bool, fill = false)
     river = river_2d[inds]
     riverwidth_2d =
-        ncread(nc, param(config, "input.lateral.river.width"); type = Float64, fill = 0)
+        ncread(nc, param(config, "input.lateral.river.width"); type = Float, fill = 0)
     riverwidth = riverwidth_2d[inds]
     riverlength_2d =
-        ncread(nc, param(config, "input.lateral.river.length"); type = Float64, fill = 0)
+        ncread(nc, param(config, "input.lateral.river.length"); type = Float, fill = 0)
     riverlength = riverlength_2d[inds]
 
     inds_riv, rev_inds_riv = active_indices(river_2d, 0)
@@ -385,7 +385,7 @@ function initialize_hbv_model(config::Config)
         ldd = set_pit_ldd(pits_2d, ldd, inds)
     end
 
-    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float64)
+    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float)
     clamp!(βₗ, 0.00001, Inf)
     dl = fill(mv, n)
     dw = fill(mv, n)
@@ -400,31 +400,31 @@ function initialize_hbv_model(config::Config)
         param(config, "input.lateral.land.n", nothing);
         sel = inds,
         defaults = 0.072,
-        type = Float64,
+        type = Float,
     )
 
-    alpha_pow = (2.0 / 3.0) * 0.6
-    β = 0.6
+    alpha_pow = Float((2.0 / 3.0) * 0.6)
+    β = Float(0.6)
     olf = SurfaceFlow(
         β = β,
         sl = βₗ,
         n = n_land,
         dl = dl,
-        q = fill(0.0, n),
-        qin = fill(0.0, n),
-        q_av = fill(0.0, n),
-        qlat = fill(0.0, n),
-        h = fill(0.0, n),
-        h_av = fill(0.0, n),
-        Δt = tosecond(Δt),
+        q = zeros(Float, n),
+        qin = zeros(Float, n),
+        q_av = zeros(Float, n),
+        qlat = zeros(Float, n),
+        h = zeros(Float, n),
+        h_av = zeros(Float, n),
+        Δt = Float(tosecond(Δt)),
         its = kw_land_tstep > 0 ? Int(cld(tosecond(Δt), kw_land_tstep)) : kw_land_tstep,
         width = dw,
         wb_pit = pits[inds],
         alpha_pow = alpha_pow,
         α = fill(mv, n),
-        eps = 1e-03,
-        cel = fill(0.0, n),
-        to_river = fill(0.0, n),
+        eps = Float(1e-3),
+        cel = zeros(Float, n),
+        to_river = zeros(Float, n),
         rivercells = fill(false, n),
         reservoir_index = fill(0, n),
         lake_index = fill(0, n),
@@ -439,7 +439,7 @@ function initialize_hbv_model(config::Config)
         nc,
         param(config, "input.lateral.river.slope");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     clamp!(riverslope, 0.00001, Inf)
     riverlength = riverlength_2d[inds_riv]
@@ -449,7 +449,7 @@ function initialize_hbv_model(config::Config)
         param(config, "input.lateral.river.n", nothing);
         sel = inds_riv,
         defaults = 0.036,
-        type = Float64,
+        type = Float,
     )
     ldd_riv = ldd_2d[inds_riv]
     if do_pits
@@ -466,22 +466,22 @@ function initialize_hbv_model(config::Config)
         sl = riverslope,
         n = n_river,
         dl = riverlength,
-        q = fill(0.0, nriv),
-        qin = fill(0.0, nriv),
-        q_av = fill(0.0, nriv),
-        qlat = fill(0.0, nriv),
-        h = fill(0.0, nriv),
-        h_av = fill(0.0, nriv),
-        Δt = tosecond(Δt),
+        q = zeros(Float, nriv),
+        qin = zeros(Float, nriv),
+        q_av = zeros(Float, nriv),
+        qlat = zeros(Float, nriv),
+        h = zeros(Float, nriv),
+        h_av = zeros(Float, nriv),
+        Δt = Float(tosecond(Δt)),
         its = kw_river_tstep > 0 ? ceil(Int(tosecond(Δt) / kw_river_tstep)) :
               kw_river_tstep,
         width = riverwidth,
         wb_pit = pits[inds_riv],
         alpha_pow = alpha_pow,
         α = fill(mv, nriv),
-        eps = 1e-03,
-        cel = fill(0.0, nriv),
-        to_river = fill(0.0, nriv),
+        eps = Float(1e-3),
+        cel = zeros(Float, nriv),
+        to_river = zeros(Float, nriv),
         reservoir_index = do_reservoirs ? resindex : fill(0, nriv),
         lake_index = do_lakes ? lakeindex : fill(0, nriv),
         reservoir = do_reservoirs ? reservoirs : nothing,
@@ -533,7 +533,7 @@ function initialize_hbv_model(config::Config)
     if reinit == false
         instate_path = joinpath(tomldir, config.state.path_input)
         state_ncnames = ncnames(config.state)
-        set_states(instate_path, model, state_ncnames; type = Float64)
+        set_states(instate_path, model, state_ncnames; type = Float)
     end
 
     # make sure the forcing is already loaded
diff --git a/src/io.jl b/src/io.jl
index 5c72a741b..5df23e3cd 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -781,7 +781,7 @@ function prepare_writer(
             time_units,
             maxlayers,
             sizeinmetres;
-            float_type = Float64,
+            float_type = Float,
         )
     else
         ds_outstate = nothing
@@ -1124,7 +1124,7 @@ end
 
 "Read a rating curve from CSV into a NamedTuple of vectors"
 function read_sh_csv(path)
-    data, header = readdlm(path, ',', Float64, header = true)
+    data, header = readdlm(path, ',', Float, header = true)
     names = vec(uppercase.(header))
     idx_h = findfirst(==("H"), names)
     idx_s = findfirst(==("S"), names)
@@ -1138,11 +1138,11 @@ end
 
 "Read a specific storage curve from CSV into a NamedTuple of vectors"
 function read_hq_csv(path)
-    data = readdlm(path, ',', Float64, skipstart = 1)
+    data = readdlm(path, ',', Float, skipstart = 1)
     # Q is a matrix with 365 columns, one for each day in the year
     return (H = data[:, 1], Q = data[:, 2:end])
 end
 
 # these represent the type of the rating curve and specific storage data
-const SH = NamedTuple{(:H, :S),Tuple{Vector{Float64},Vector{Float64}}}
-const HQ = NamedTuple{(:H, :Q),Tuple{Vector{Float64},Matrix{Float64}}}
+const SH = NamedTuple{(:H, :S),Tuple{Vector{Float},Vector{Float}}}
+const HQ = NamedTuple{(:H, :Q),Tuple{Vector{Float},Matrix{Float}}}
diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index 6171501c4..f83ad8756 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -70,42 +70,42 @@ function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits)
         nc,
         param(config, "input.lateral.river.reservoir.demand");
         sel = inds_res,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
     resmaxrelease = ncread(
         nc,
         param(config, "input.lateral.river.reservoir.maxrelease");
         sel = inds_res,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
     resmaxvolume = ncread(
         nc,
         param(config, "input.lateral.river.reservoir.maxvolume");
         sel = inds_res,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
     resarea = ncread(
         nc,
         param(config, "input.lateral.river.reservoir.area");
         sel = inds_res,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
     res_targetfullfrac = ncread(
         nc,
         param(config, "input.lateral.river.reservoir.targetfullfrac");
         sel = inds_res,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
     res_targetminfrac = ncread(
         nc,
         param(config, "input.lateral.river.reservoir.targetminfrac");
         sel = inds_res,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
 
@@ -114,7 +114,7 @@ function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits)
     pits[inds_res] .= true
 
     n = length(resarea)
-    reservoirs = SimpleReservoir{Float64}(
+    reservoirs = SimpleReservoir{Float}(
         demand = resdemand,
         maxrelease = resmaxrelease,
         maxvolume = resmaxvolume,
@@ -251,28 +251,28 @@ function initialize_natural_lake(config, path, nc, inds_riv, nriv, pits)
         nc,
         param(config, "input.lateral.river.lake.area");
         sel = inds_lake,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
     lake_b = ncread(
         nc,
         param(config, "input.lateral.river.lake.b");
         sel = inds_lake,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
     lake_e = ncread(
         nc,
         param(config, "input.lateral.river.lake.e");
         sel = inds_lake,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
     lake_threshold = ncread(
         nc,
         param(config, "input.lateral.river.lake.threshold");
         sel = inds_lake,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
     linked_lakelocs = ncread(
@@ -301,7 +301,7 @@ function initialize_natural_lake(config, path, nc, inds_riv, nriv, pits)
         nc,
         param(config, "input.lateral.river.lake.waterlevel");
         sel = inds_lake,
-        type = Float64,
+        type = Float,
         fill = 0,
     )
 
@@ -341,7 +341,7 @@ function initialize_natural_lake(config, path, nc, inds_riv, nriv, pits)
         end
     end
     n = length(lakearea)
-    lakes = NaturalLake{Float64}(
+    lakes = NaturalLake{Float}(
         lowerlake_ind = linked_lakelocs,
         area = lakearea,
         threshold = lake_threshold,
diff --git a/src/sbm.jl b/src/sbm.jl
index d6bbe6ff7..daed1925c 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -212,15 +212,15 @@ function initialize_canopy(nc, config, inds)
             nc,
             param(config, "input.vertical.specific_leaf");
             sel = inds,
-            type = Float64,
+            type = Float,
         )
         swood = ncread(
             nc,
             param(config, "input.vertical.storage_wood");
             sel = inds,
-            type = Float64,
+            type = Float,
         )
-        kext = ncread(nc, param(config, "input.vertical.kext"); sel = inds, type = Float64)
+        kext = ncread(nc, param(config, "input.vertical.kext"); sel = inds, type = Float)
         cmax = fill(mv, n)
         e_r = fill(mv, n)
         canopygapfraction = fill(mv, n)
@@ -234,21 +234,21 @@ function initialize_canopy(nc, config, inds)
             param(config, "input.vertical.cmax", nothing);
             sel = inds,
             defaults = 1.0,
-            type = Float64,
+            type = Float,
         )
         e_r = ncread(
             nc,
             param(config, "input.vertical.eoverr", nothing);
             sel = inds,
             defaults = 0.1,
-            type = Float64,
+            type = Float,
         )
         canopygapfraction = ncread(
             nc,
             param(config, "input.vertical.canopygapfraction", nothing);
             sel = inds,
             defaults = 0.1,
-            type = Float64,
+            type = Float,
         )
     end
     return cmax, e_r, canopygapfraction, sl, swood, kext
@@ -257,9 +257,9 @@ end
 function initialize_sbm(nc, config, riverfrac, inds)
 
     Δt = Second(config.timestepsecs)
-    config_thicknesslayers = get(config.model, "thicknesslayers", Float64[])
+    config_thicknesslayers = get(config.model, "thicknesslayers", Float[])
     if length(config_thicknesslayers) > 0
-        thicknesslayers = SVector(Tuple(push!(Float64.(config_thicknesslayers), mv)))
+        thicknesslayers = SVector(Tuple(push!(Float.(config_thicknesslayers), mv)))
         sumlayers = pushfirst(cumsum(thicknesslayers), 0.0)
         maxlayers = length(thicknesslayers) # max number of soil layers
     else
@@ -274,35 +274,35 @@ function initialize_sbm(nc, config, riverfrac, inds)
             param(config, "input.vertical.cfmax", nothing);
             sel = inds,
             defaults = 3.75653,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     tt = ncread(
         nc,
         param(config, "input.vertical.tt", nothing);
         sel = inds,
         defaults = 0.0,
-        type = Float64,
+        type = Float,
     )
     tti = ncread(
         nc,
         param(config, "input.vertical.tti", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
     )
     ttm = ncread(
         nc,
         param(config, "input.vertical.ttm", nothing);
         sel = inds,
         defaults = 0.0,
-        type = Float64,
+        type = Float,
     )
     whc = ncread(
         nc,
         param(config, "input.vertical.whc", nothing);
         sel = inds,
         defaults = 0.1,
-        type = Float64,
+        type = Float,
     )
     w_soil =
         ncread(
@@ -310,14 +310,14 @@ function initialize_sbm(nc, config, riverfrac, inds)
             param(config, "input.vertical.w_soil", nothing);
             sel = inds,
             defaults = 0.1125,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     cf_soil = ncread(
         nc,
         param(config, "input.vertical.cf_soil", nothing);
         sel = inds,
         defaults = 0.038,
-        type = Float64,
+        type = Float,
     )
     # glacier parameters
     g_tt = ncread(
@@ -325,7 +325,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         param(config, "input.vertical.g_tt", nothing);
         sel = inds,
         defaults = 0.0,
-        type = Float64,
+        type = Float,
         fill = 0.0,
     )
     g_cfmax =
@@ -334,7 +334,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             param(config, "input.vertical.g_cfmax", nothing);
             sel = inds,
             defaults = 3.0,
-            type = Float64,
+            type = Float,
             fill = 0.0,
         ) .* (Δt / basetimestep)
     g_sifrac = ncread(
@@ -342,7 +342,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         param(config, "input.vertical.g_sifrac", nothing);
         sel = inds,
         defaults = 0.001,
-        type = Float64,
+        type = Float,
         fill = 0.0,
     )
     glacierfrac = ncread(
@@ -350,7 +350,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         param(config, "input.vertical.glacierfrac", nothing);
         sel = inds,
         defaults = 0.0,
-        type = Float64,
+        type = Float,
         fill = 0.0,
     )
     glacierstore = ncread(
@@ -358,7 +358,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         param(config, "input.vertical.glacierstore", nothing);
         sel = inds,
         defaults = 5500.0,
-        type = Float64,
+        type = Float,
         fill = 0.0,
     )
     # soil parameters
@@ -367,14 +367,14 @@ function initialize_sbm(nc, config, riverfrac, inds)
         param(config, "input.vertical.θₛ", nothing);
         sel = inds,
         defaults = 0.6,
-        type = Float64,
+        type = Float,
     )
     θᵣ = ncread(
         nc,
         param(config, "input.vertical.θᵣ", nothing);
         sel = inds,
         defaults = 0.01,
-        type = Float64,
+        type = Float,
     )
     kv₀ =
         ncread(
@@ -382,28 +382,28 @@ function initialize_sbm(nc, config, riverfrac, inds)
             param(config, "input.vertical.kv₀", nothing);
             sel = inds,
             defaults = 3000.0,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     f = ncread(
         nc,
         param(config, "input.vertical.f", nothing);
         sel = inds,
         defaults = 0.001,
-        type = Float64,
+        type = Float,
     )
     hb = ncread(
         nc,
         param(config, "input.vertical.hb", nothing);
         sel = inds,
         defaults = 10.0,
-        type = Float64,
+        type = Float,
     )
     soilthickness = ncread(
         nc,
         param(config, "input.vertical.soilthickness", nothing);
         sel = inds,
         defaults = 2000.0,
-        type = Float64,
+        type = Float,
     )
     infiltcappath =
         ncread(
@@ -411,7 +411,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             param(config, "input.vertical.infiltcappath", nothing);
             sel = inds,
             defaults = 10.0,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     infiltcapsoil =
         ncread(
@@ -419,7 +419,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             param(config, "input.vertical.infiltcapsoil", nothing);
             sel = inds,
             defaults = 100.0,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
     maxleakage =
         ncread(
@@ -427,7 +427,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             param(config, "input.vertical.maxleakage", nothing);
             sel = inds,
             defaults = 0.0,
-            type = Float64,
+            type = Float,
         ) .* (Δt / basetimestep)
 
     c = ncread(
@@ -435,7 +435,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         param(config, "input.vertical.c", nothing);
         sel = inds,
         defaults = 10.0,
-        type = Float64,
+        type = Float,
         dimname = "layer",
     )
     if size(c, 1) != maxlayers
@@ -448,7 +448,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         param(config, "input.vertical.khfrac", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
         dimname = "layer",
     )
     if size(kvfrac, 1) != maxlayers
@@ -463,14 +463,14 @@ function initialize_sbm(nc, config, riverfrac, inds)
         param(config, "input.vertical.waterfrac", nothing);
         sel = inds,
         defaults = 0.0,
-        type = Float64,
+        type = Float,
     )
     pathfrac = ncread(
         nc,
         param(config, "input.vertical.pathfrac", nothing);
         sel = inds,
         defaults = 0.01,
-        type = Float64,
+        type = Float,
     )
 
     # vegetation parameters
@@ -479,28 +479,28 @@ function initialize_sbm(nc, config, riverfrac, inds)
         param(config, "input.vertical.rootingdepth", nothing);
         sel = inds,
         defaults = 750.0,
-        type = Float64,
+        type = Float,
     )
     rootdistpar = ncread(
         nc,
         param(config, "input.vertical.rootdistpar", nothing);
         sel = inds,
         defaults = -500.0,
-        type = Float64,
+        type = Float,
     )
     capscale = ncread(
         nc,
         param(config, "input.vertical.capscale", nothing);
         sel = inds,
         defaults = 100.0,
-        type = Float64,
+        type = Float,
     )
     et_reftopot = ncread(
         nc,
         param(config, "input.vertical.et_reftopot", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
     )
 
     cmax, e_r, canopygapfraction, sl, swood, kext = initialize_canopy(nc, config, inds)
@@ -508,8 +508,8 @@ function initialize_sbm(nc, config, riverfrac, inds)
     # these are filled in the loop below
     # TODO see if we can replace this approach
     nlayers = zeros(Int, n)
-    act_thickl = zeros(Float64, maxlayers, n)
-    s_layers = zeros(Float64, maxlayers + 1, n)
+    act_thickl = zeros(Float, maxlayers, n)
+    s_layers = zeros(Float, maxlayers + 1, n)
 
     for i = 1:n
         if length(config_thicknesslayers) > 0
@@ -538,7 +538,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
     vwc = fill(mv, maxlayers, n)
     vwc_perc = fill(mv, maxlayers, n)
 
-    sbm = SBM{Float64,maxlayers,maxlayers + 1}(
+    sbm = SBM{Float,maxlayers,maxlayers + 1}(
         Δt = tosecond(Δt),
         maxlayers = maxlayers,
         n = n,
@@ -556,7 +556,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         infiltcappath = infiltcappath,
         infiltcapsoil = infiltcapsoil,
         maxleakage = maxleakage,
-        waterfrac = max.(waterfrac .- riverfrac, 0.0),
+        waterfrac = max.(waterfrac .- riverfrac, Float(0.0)),
         pathfrac = pathfrac,
         rootingdepth = rootingdepth,
         rootdistpar = rootdistpar,
@@ -570,7 +570,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         satwaterdepth = satwaterdepth,
         zi = zi,
         soilwatercapacity = soilwatercapacity,
-        canopystorage = fill(0.0, n),
+        canopystorage = zeros(Float, n),
         cmax = cmax,
         canopygapfraction = canopygapfraction,
         e_r = e_r,
@@ -620,10 +620,10 @@ function initialize_sbm(nc, config, riverfrac, inds)
         whc = whc,
         w_soil = w_soil,
         cf_soil = cf_soil,
-        snow = fill(0.0, n),
-        snowwater = fill(0.0, n),
+        snow = zeros(Float, n),
+        snowwater = zeros(Float, n),
         rainfallplusmelt = fill(mv, n),
-        tsoil = fill(10.0, n),
+        tsoil = fill(Float(10.0), n),
         # glacier parameters
         g_tt = g_tt,
         g_sifrac = g_sifrac,
@@ -636,7 +636,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         kext = kext,
         leaf_area_index = fill(mv, n),
         waterlevel_land = fill(mv, n),
-        waterlevel_river = fill(0.0, n), #set to zero to account for cells outside river domain
+        waterlevel_river = zeros(Float, n), #set to zero to account for cells outside river domain
     )
 
     return sbm
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index 6859c5e0c..ce0a21e9a 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -53,14 +53,14 @@ function initialize_sbm_gwf_model(config::Config)
     river_2d = ncread(nc, param(config, "input.river_location"); type = Bool, fill = false)
     river = river_2d[inds]
     riverwidth_2d =
-        ncread(nc, param(config, "input.lateral.river.width"); type = Float64, fill = 0)
+        ncread(nc, param(config, "input.lateral.river.width"); type = Float, fill = 0)
     riverwidth = riverwidth_2d[inds]
     riverlength_2d =
-        ncread(nc, param(config, "input.lateral.river.length"); type = Float64, fill = 0)
+        ncread(nc, param(config, "input.lateral.river.length"); type = Float, fill = 0)
     riverlength = riverlength_2d[inds]
 
     altitude =
-        ncread(nc, param(config, "input.altitude"); sel = inds, type = Float64)
+        ncread(nc, param(config, "input.altitude"); sel = inds, type = Float)
     # read x, y coordinates and calculate cell length [m]
     y_nc = "y" in keys(nc.dim) ? ncread(nc, "y") : ncread(nc, "lat")
     x_nc = "x" in keys(nc.dim) ? ncread(nc, "x") : ncread(nc, "lon")
@@ -106,7 +106,7 @@ function initialize_sbm_gwf_model(config::Config)
     end
 
     # overland flow (kinematic wave)
-    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float64)
+    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float)
     clamp!(βₗ, 0.00001, Inf)
     ldd_2d = ncread(nc, param(config, "input.ldd"); allow_missing = true)
     ldd = ldd_2d[inds]
@@ -125,31 +125,31 @@ function initialize_sbm_gwf_model(config::Config)
         param(config, "input.lateral.land.n", nothing);
         sel = inds,
         defaults = 0.072,
-        type = Float64,
+        type = Float,
     )
 
-    alpha_pow = (2.0 / 3.0) * 0.6
-    β = 0.6
+    alpha_pow = Float((2.0 / 3.0) * 0.6)
+    β = Float(0.6)
     olf = SurfaceFlow(
         β = β,
         sl = βₗ,
         n = n_land,
         dl = dl,
-        q = fill(0.0, n),
-        qin = fill(0.0, n),
-        q_av = fill(0.0, n),
-        qlat = fill(0.0, n),
-        h = fill(0.0, n),
-        h_av = fill(0.0, n),
-        Δt = tosecond(Δt),
+        q = zeros(Float, n),
+        qin = zeros(Float, n),
+        q_av = zeros(Float, n),
+        qlat = zeros(Float, n),
+        h = zeros(Float, n),
+        h_av = zeros(Float, n),
+        Δt = Float(tosecond(Δt)),
         its = kw_land_tstep > 0 ? Int(cld(tosecond(Δt), kw_land_tstep)) : kw_land_tstep,
         width = sw,
         wb_pit = pits[inds],
         alpha_pow = alpha_pow,
         α = fill(mv, n),
-        eps = 1e-03,
-        cel = fill(0.0, n),
-        to_river = fill(0.0, n),
+        eps = Float(1e-03),
+        cel = zeros(Float, n),
+        to_river = zeros(Float, n),
         rivercells = fill(false, n),
         reservoir_index = fill(0, n),
         lake_index = fill(0, n),
@@ -165,7 +165,7 @@ function initialize_sbm_gwf_model(config::Config)
         nc,
         param(config, "input.lateral.river.slope");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     clamp!(riverslope, 0.00001, Inf)
     riverlength = riverlength_2d[inds_riv]
@@ -175,7 +175,7 @@ function initialize_sbm_gwf_model(config::Config)
         param(config, "input.lateral.river.n", nothing);
         sel = inds_riv,
         defaults = 0.036,
-        type = Float64,
+        type = Float,
     )
     ldd_riv = ldd_2d[inds_riv]
     graph_riv = flowgraph(ldd_riv, inds_riv, pcr_dir)
@@ -189,22 +189,22 @@ function initialize_sbm_gwf_model(config::Config)
         sl = riverslope,
         n = n_river,
         dl = riverlength,
-        q = fill(0.0, nriv),
-        q_av = fill(0.0, nriv),
-        qin = fill(0.0, nriv),
-        qlat = fill(0.0, nriv),
-        h = fill(0.0, nriv),
-        h_av = fill(0.0, nriv),
-        Δt = tosecond(Δt),
+        q = zeros(Float, nriv),
+        q_av = zeros(Float, nriv),
+        qin = zeros(Float, nriv),
+        qlat = zeros(Float, nriv),
+        h = zeros(Float, nriv),
+        h_av = zeros(Float, nriv),
+        Δt = Float(tosecond(Δt)),
         its = kw_river_tstep > 0 ? ceil(Int(tosecond(Δt) / kw_river_tstep)) :
               kw_river_tstep,
         width = riverwidth,
         wb_pit = pits[inds_riv],
         alpha_pow = alpha_pow,
         α = fill(mv, nriv),
-        eps = 1e-03,
-        cel = fill(0.0, nriv),
-        to_river = fill(0.0, nriv),
+        eps = Float(1e-03),
+        cel = zeros(Float, nriv),
+        to_river = zeros(Float, nriv),
         reservoir_index = do_reservoirs ? resindex : fill(0, nriv),
         lake_index = do_lakes ? lakeindex : fill(0, nriv),
         reservoir = do_reservoirs ? reservoirs : nothing,
@@ -218,7 +218,7 @@ function initialize_sbm_gwf_model(config::Config)
             nc,
             param(config, "input.lateral.subsurface.constant_head", nothing);
             sel = inds,
-            type = Float64,
+            type = Float,
             fill = mv,
         )
         index_constanthead = filter(i -> !isequal(constanthead[i], mv), 1:n)
@@ -229,17 +229,17 @@ function initialize_sbm_gwf_model(config::Config)
         nc,
         param(config, "input.lateral.subsurface.conductivity", nothing);
         sel = inds,
-        type = Float64,
+        type = Float,
     )
     specific_yield = ncread(
         nc,
         param(config, "input.lateral.subsurface.specific_yield", nothing);
         sel = inds,
-        type = Float64,
+        type = Float,
     )
 
     connectivity = Connectivity(inds, rev_inds, xl, yl)
-    initial_head = altitude .- 0.10 # cold state for groundwater head
+    initial_head = altitude .- Float(0.10) # cold state for groundwater head
     initial_head[index_river] = altitude[index_river]
 
     if do_constanthead
@@ -250,10 +250,10 @@ function initialize_sbm_gwf_model(config::Config)
         initial_head,
         conductivity,
         altitude,
-        altitude .- sbm.soilthickness ./ 1000.0,
+        altitude .- sbm.soilthickness ./ Float(1000.0),
         xl .* yl,
         specific_yield,
-        fill(0.0, connectivity.nconnection),  # conductance
+        zeros(Float, connectivity.nconnection),  # conductance
     )
 
     # reset zi and satwaterdepth with groundwater head from unconfined aquifer 
@@ -265,19 +265,19 @@ function initialize_sbm_gwf_model(config::Config)
         nc,
         param(config, "input.lateral.subsurface.infiltration_conductance", nothing);
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     exfiltration_conductance = ncread(
         nc,
         param(config, "input.lateral.subsurface.exfiltration_conductance", nothing);
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     river_bottom = ncread(
         nc,
         param(config, "input.lateral.subsurface.river_bottom", nothing);
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
 
     river_flux = fill(mv, nriv)
@@ -293,7 +293,7 @@ function initialize_sbm_gwf_model(config::Config)
 
     # recharge boundary of unconfined aquifer
     r = fill(mv, n)
-    recharge = Recharge(r, fill(0.0, n), collect(1:n))
+    recharge = Recharge(r, zeros(Float, n), collect(1:n))
 
     # drain boundary of unconfined aquifer (optional)
     if do_drains
@@ -309,14 +309,14 @@ function initialize_sbm_gwf_model(config::Config)
             nc,
             param(config, "input.lateral.subsurface.drain_elevation", nothing);
             sel = inds,
-            type = Float64,
+            type = Float,
             fill = mv,
         )
         drain_conductance = ncread(
             nc,
             param(config, "input.lateral.subsurface.drain_conductance", nothing);
             sel = inds,
-            type = Float64,
+            type = Float,
             fill = mv,
         )
         index_drain = filter(i -> !isequal(drain[i], 0), 1:n)
@@ -407,7 +407,7 @@ function initialize_sbm_gwf_model(config::Config)
     if reinit == false
         instate_path = joinpath(tomldir, config.state.path_input)
         state_ncnames = ncnames(config.state)
-        set_states(instate_path, model, state_ncnames, type = Float64)
+        set_states(instate_path, model, state_ncnames, type = Float)
     end
 
     # make sure the forcing is already loaded
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index 7bdb5daaa..3f187cd2f 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -43,10 +43,10 @@ function initialize_sbm_model(config::Config)
     river_2d = ncread(nc, param(config, "input.river_location"); type = Bool, fill = false)
     river = river_2d[inds]
     riverwidth_2d =
-        ncread(nc, param(config, "input.lateral.river.width"); type = Float64, fill = 0)
+        ncread(nc, param(config, "input.lateral.river.width"); type = Float, fill = 0)
     riverwidth = riverwidth_2d[inds]
     riverlength_2d =
-        ncread(nc, param(config, "input.lateral.river.length"); type = Float64, fill = 0)
+        ncread(nc, param(config, "input.lateral.river.length"); type = Float, fill = 0)
     riverlength = riverlength_2d[inds]
 
     # read x, y coordinates and calculate cell length [m]
@@ -106,9 +106,9 @@ function initialize_sbm_model(config::Config)
         param(config, "input.lateral.subsurface.ksathorfrac", nothing);
         sel = inds,
         defaults = 1.0,
-        type = Float64,
+        type = Float,
     )
-    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float64)
+    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float)
     clamp!(βₗ, 0.00001, Inf)
     kh₀ = khfrac .* sbm.kv₀
     dl = fill(mv, n)
@@ -121,7 +121,7 @@ function initialize_sbm_model(config::Config)
         sw[i] = river[i] ? max(dw[i] - riverwidth[i], 0.0) : dw[i]
     end
 
-    ssf = LateralSSF{Float64}(
+    ssf = LateralSSF{Float}(
         kh₀ = kh₀,
         f = sbm.f,
         zi = sbm.zi,
@@ -149,31 +149,31 @@ function initialize_sbm_model(config::Config)
         param(config, "input.lateral.land.n", nothing);
         sel = inds,
         defaults = 0.072,
-        type = Float64,
+        type = Float,
     )
 
-    alpha_pow = (2.0 / 3.0) * 0.6
-    β = 0.6
+    alpha_pow = Float((2.0 / 3.0) * 0.6)
+    β = Float(0.6)
     olf = SurfaceFlow(
         β = β,
         sl = βₗ,
         n = n_land,
         dl = dl,
-        q = fill(0.0, n),
-        qin = fill(0.0, n),
-        q_av = fill(0.0, n),
-        qlat = fill(0.0, n),
-        h = fill(0.0, n),
-        h_av = fill(0.0, n),
-        Δt = tosecond(Δt),
+        q = zeros(Float, n),
+        qin = zeros(Float, n),
+        q_av = zeros(Float, n),
+        qlat = zeros(Float, n),
+        h = zeros(Float, n),
+        h_av = zeros(Float, n),
+        Δt = Float(tosecond(Δt)),
         its = kw_land_tstep > 0 ? Int(cld(tosecond(Δt), kw_land_tstep)) : kw_land_tstep,
         width = sw,
         wb_pit = pits[inds],
         alpha_pow = alpha_pow,
         α = fill(mv, n),
-        eps = 1e-03,
-        cel = fill(0.0, n),
-        to_river = fill(0.0, n),
+        eps = Float(1e-3),
+        cel = zeros(Float, n),
+        to_river = zeros(Float, n),
         rivercells = fill(false, n),
         reservoir_index = fill(0, n),
         lake_index = fill(0, n),
@@ -189,7 +189,7 @@ function initialize_sbm_model(config::Config)
         nc,
         param(config, "input.lateral.river.slope");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     clamp!(riverslope, 0.00001, Inf)
     riverlength = riverlength_2d[inds_riv]
@@ -199,7 +199,7 @@ function initialize_sbm_model(config::Config)
         param(config, "input.lateral.river.n", nothing);
         sel = inds_riv,
         defaults = 0.036,
-        type = Float64,
+        type = Float,
     )
     ldd_riv = ldd_2d[inds_riv]
     if do_pits
@@ -216,22 +216,22 @@ function initialize_sbm_model(config::Config)
         sl = riverslope,
         n = n_river,
         dl = riverlength,
-        q = fill(0.0, nriv),
-        qin = fill(0.0, nriv),
-        q_av = fill(0.0, nriv),
-        qlat = fill(0.0, nriv),
-        h = fill(0.0, nriv),
-        h_av = fill(0.0, nriv),
-        Δt = tosecond(Δt),
+        q = zeros(Float, nriv),
+        qin = zeros(Float, nriv),
+        q_av = zeros(Float, nriv),
+        qlat = zeros(Float, nriv),
+        h = zeros(Float, nriv),
+        h_av = zeros(Float, nriv),
+        Δt = Float(tosecond(Δt)),
         its = kw_river_tstep > 0 ? ceil(Int(tosecond(Δt) / kw_river_tstep)) :
               kw_river_tstep,
         width = riverwidth,
         wb_pit = pits[inds_riv],
         alpha_pow = alpha_pow,
         α = fill(mv, nriv),
-        eps = 1e-03,
-        cel = fill(0.0, nriv),
-        to_river = fill(0.0, nriv),
+        eps = Float(1e-3),
+        cel = zeros(Float, nriv),
+        to_river = zeros(Float, nriv),
         reservoir_index = do_reservoirs ? resindex : fill(0, nriv),
         lake_index = do_lakes ? lakeindex : fill(0, nriv),
         reservoir = do_reservoirs ? reservoirs : nothing,
@@ -294,7 +294,7 @@ function initialize_sbm_model(config::Config)
     if reinit == false
         instate_path = joinpath(tomldir, config.state.path_input)
         state_ncnames = ncnames(config.state)
-        set_states(instate_path, model, state_ncnames; type = Float64)
+        set_states(instate_path, model, state_ncnames; type = Float)
         @unpack lateral, vertical = model
         # update zi for vertical sbm and α for river and overland flow
         zi = max.(0.0, vertical.soilthickness .- vertical.satwaterdepth ./ (vertical.θₛ .- vertical.θᵣ))
diff --git a/src/sediment.jl b/src/sediment.jl
index 939bb3273..726525b4b 100644
--- a/src/sediment.jl
+++ b/src/sediment.jl
@@ -114,62 +114,62 @@ function initialize_landsed(nc, config, river, riverfrac, xl, yl, inds)
     landtransportmethod = get(config.model, "landtransportmethod", "yalinpart")::String
 
     altitude =
-        ncread(nc, param(config, "input.vertical.altitude"); sel = inds, type = Float64)
+        ncread(nc, param(config, "input.vertical.altitude"); sel = inds, type = Float)
     canopyheight = ncread(
         nc,
         param(config, "input.vertical.canopyheight", nothing);
         sel = inds,
         defaults = 3.0,
-        type = Float64,
+        type = Float,
     )
     erosk = ncread(
         nc,
         param(config, "input.vertical.erosk", nothing);
         sel = inds,
         defaults = 0.6,
-        type = Float64,
+        type = Float,
     )
     erosspl = ncread(
         nc,
         param(config, "input.vertical.erosspl", nothing);
         sel = inds,
         defaults = 2.0,
-        type = Float64,
+        type = Float,
     )
     erosov = ncread(
         nc,
         param(config, "input.vertical.erosov", nothing);
         sel = inds,
         defaults = 0.9,
-        type = Float64,
+        type = Float,
     )
     pathfrac = ncread(
         nc,
         param(config, "input.vertical.pathfrac", nothing);
         sel = inds,
         defaults = 0.01,
-        type = Float64,
+        type = Float,
     )
     slope = ncread(
         nc,
         param(config, "input.vertical.slope", nothing);
         sel = inds,
         defaults = 0.01,
-        type = Float64,
+        type = Float,
     )
     usleC = ncread(
         nc,
         param(config, "input.vertical.usleC", nothing);
         sel = inds,
         defaults = 0.01,
-        type = Float64,
+        type = Float,
     )
     usleK = ncread(
         nc,
         param(config, "input.vertical.usleK", nothing);
         sel = inds,
         defaults = 0.1,
-        type = Float64,
+        type = Float,
     )
 
     cmax, _, canopygapfraction, sl, swood, kext = initialize_canopy(nc, config, inds)
@@ -191,56 +191,56 @@ function initialize_landsed(nc, config, river, riverfrac, xl, yl, inds)
         param(config, "input.vertical.dmclay", nothing);
         sel = inds,
         defaults = 2.0,
-        type = Float64,
+        type = Float,
     )
     dmsilt = ncread(
         nc,
         param(config, "input.vertical.dmsilt", nothing);
         sel = inds,
         defaults = 10.0,
-        type = Float64,
+        type = Float,
     )
     dmsand = ncread(
         nc,
         param(config, "input.vertical.dmsand", nothing);
         sel = inds,
         defaults = 200.0,
-        type = Float64,
+        type = Float,
     )
     dmsagg = ncread(
         nc,
         param(config, "input.vertical.dmsagg", nothing);
         sel = inds,
         defaults = 30.0,
-        type = Float64,
+        type = Float,
     )
     dmlagg = ncread(
         nc,
         param(config, "input.vertical.dmlagg", nothing);
         sel = inds,
         defaults = 500.0,
-        type = Float64,
+        type = Float,
     )
     pclay = ncread(
         nc,
         param(config, "input.vertical.pclay", nothing);
         sel = inds,
         defaults = 0.1,
-        type = Float64,
+        type = Float,
     )
     psilt = ncread(
         nc,
         param(config, "input.vertical.psilt", nothing);
         sel = inds,
         defaults = 0.1,
-        type = Float64,
+        type = Float,
     )
     rhos = ncread(
         nc,
         param(config, "input.vertical.rhosed", nothing);
         sel = inds,
         defaults = 2650.0,
-        type = Float64,
+        type = Float,
     )
 
     ### Initialize transport capacity variables ###
@@ -296,13 +296,13 @@ function initialize_landsed(nc, config, river, riverfrac, xl, yl, inds)
     end
 
     # Reservoir and lakes
-    wbcover = fill(0.0, n)
+    wbcover = zeros(Float, n)
     if do_reservoirs
         rescoverage_2d = ncread(
             nc,
             param(config, "input.vertical.resareas");
             sel = inds,
-            type = Float64,
+            type = Float,
             fill = 0.0,
         )
         wbcover = wbcover .+ rescoverage_2d
@@ -312,13 +312,13 @@ function initialize_landsed(nc, config, river, riverfrac, xl, yl, inds)
             nc,
             param(config, "input.vertical.lakeareas");
             sel = inds,
-            type = Float64,
+            type = Float,
             fill = 0.0,
         )
         wbcover = wbcover .+ lakecoverage_2d
     end
 
-    eros = LandSediment{Float64}(
+    eros = LandSediment{Float}(
         n = n,
         yl = yl,
         xl = xl,
@@ -391,7 +391,7 @@ function update_until_ols(eros::LandSediment, config, network)
     do_lai = haskey(config.input.vertical, "leaf_area_index")
     rainerosmethod = get(config.model, "rainerosmethod", "answers")::String
     Δt = Second(config.timestepsecs)
-    ts = Float64(Δt.value)
+    ts = Float(Δt.value)
 
     for i = 1:eros.n
 
@@ -487,7 +487,7 @@ function update_until_oltransport(ols::LandSediment, config, network)
     do_river = get(config.model, "runrivermodel", false)::Bool
     tcmethod = get(config.model, "landtransportmethod", "yalinpart")::String
     Δt = Second(config.timestepsecs)
-    ts = Float64(Δt.value)
+    ts = Float(Δt.value)
 
     for i = 1:ols.n
         sinslope = sin(atan(ols.slope[i]))
@@ -885,38 +885,38 @@ function initialize_riversed(nc, config, riverwidth, riverlength, inds_riv)
     # Reservoir / lakes
     do_reservoirs = get(config.model, "doreservoir", false)::Bool
     do_lakes = get(config.model, "dolake", false)::Bool
-    wbcover = fill(0.0, nriv)
-    wblocs = fill(0.0, nriv)
-    wbarea = fill(0.0, nriv)
-    wbtrap = fill(0.0, nriv)
+    wbcover = zeros(Float, nriv)
+    wblocs = zeros(Float, nriv)
+    wbarea = zeros(Float, nriv)
+    wbtrap = zeros(Float, nriv)
 
     if do_reservoirs
         reslocs = ncread(
             nc,
             param(config, "input.lateral.river.reslocs");
             sel = inds_riv,
-            type = Float64,
+            type = Float,
             fill = 0,
         )
         rescoverage_2d = ncread(
             nc,
             param(config, "input.lateral.river.resareas");
             sel = inds_riv,
-            type = Float64,
+            type = Float,
             fill = 0,
         )
         resarea = ncread(
             nc,
             param(config, "input.lateral.river.resarea");
             sel = inds_riv,
-            type = Float64,
+            type = Float,
             fill = 0.0,
         )
         restrapefficiency = ncread(
             nc,
             param(config, "input.lateral.river.restrapeff");
             sel = inds_riv,
-            type = Float64,
+            type = Float,
             defaults = 1.0,
             fill = 0.0,
         )
@@ -932,21 +932,21 @@ function initialize_riversed(nc, config, riverwidth, riverlength, inds_riv)
             nc,
             param(config, "input.lateral.river.lakelocs");
             sel = inds_riv,
-            type = Float64,
+            type = Float,
             fill = 0,
         )
         lakecoverage_2d = ncread(
             nc,
             param(config, "input.lateral.river.lakeareas");
             sel = inds_riv,
-            type = Float64,
+            type = Float,
             fill = 0,
         )
         lakearea = ncread(
             nc,
             param(config, "input.lateral.river.lakearea");
             sel = inds_riv,
-            type = Float64,
+            type = Float,
             fill = 0.0,
         )
 
@@ -959,7 +959,7 @@ function initialize_riversed(nc, config, riverwidth, riverlength, inds_riv)
         nc,
         param(config, "input.lateral.river.slope");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     clamp!(riverslope, 0.00001, Inf)
     rhos = ncread(
@@ -967,100 +967,100 @@ function initialize_riversed(nc, config, riverwidth, riverlength, inds_riv)
         param(config, "input.lateral.river.rhosed", nothing);
         sel = inds_riv,
         defaults = 2650.0,
-        type = Float64,
+        type = Float,
     )
     dmclay = ncread(
         nc,
         param(config, "input.lateral.river.dmclay", nothing);
         sel = inds_riv,
         defaults = 2.0,
-        type = Float64,
+        type = Float,
     )
     dmsilt = ncread(
         nc,
         param(config, "input.lateral.river.dmsilt", nothing);
         sel = inds_riv,
         defaults = 10.0,
-        type = Float64,
+        type = Float,
     )
     dmsand = ncread(
         nc,
         param(config, "input.lateral.river.dmsand", nothing);
         sel = inds_riv,
         defaults = 200.0,
-        type = Float64,
+        type = Float,
     )
     dmsagg = ncread(
         nc,
         param(config, "input.lateral.river.dmsagg", nothing);
         sel = inds_riv,
         defaults = 30.0,
-        type = Float64,
+        type = Float,
     )
     dmlagg = ncread(
         nc,
         param(config, "input.lateral.river.dmlagg", nothing);
         sel = inds_riv,
         defaults = 500.0,
-        type = Float64,
+        type = Float,
     )
     dmgrav = ncread(
         nc,
         param(config, "input.lateral.river.dmgrav", nothing);
         sel = inds_riv,
         defaults = 2000.0,
-        type = Float64,
+        type = Float,
     )
     fclayriv = ncread(
         nc,
         param(config, "input.lateral.river.fclayriv");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     fsiltriv = ncread(
         nc,
         param(config, "input.lateral.river.fsiltriv");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     fsandriv = ncread(
         nc,
         param(config, "input.lateral.river.fsandriv");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     fgravriv = ncread(
         nc,
         param(config, "input.lateral.river.fgravriv");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     d50riv =
-        ncread(nc, param(config, "input.lateral.river.d50"); sel = inds_riv, type = Float64)
+        ncread(nc, param(config, "input.lateral.river.d50"); sel = inds_riv, type = Float)
     d50engelund = ncread(
         nc,
         param(config, "input.lateral.river.d50engelund");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     cbagnold = ncread(
         nc,
         param(config, "input.lateral.river.cbagnold");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
     ebagnold = ncread(
         nc,
         param(config, "input.lateral.river.ebagnold");
         sel = inds_riv,
-        type = Float64,
+        type = Float,
     )
 
     # Initialisation of parameters for Kodatie transport capacity
-    ak = fill(0.0, nriv)
-    bk = fill(0.0, nriv)
-    ck = fill(0.0, nriv)
-    dk = fill(0.0, nriv)
+    ak = zeros(Float, nriv)
+    bk = zeros(Float, nriv)
+    ck = zeros(Float, nriv)
+    dk = zeros(Float, nriv)
     if tcmethodriv == "kodatie"
         for i = 1:nriv
             if d50riv[i] <= 0.05
@@ -1090,20 +1090,20 @@ function initialize_riversed(nc, config, riverwidth, riverlength, inds_riv)
     # Bed and Bank from Shields diagram, Da Silva & Yalin (2017)
     E_ = (2.65 - 1) * 9.81
     E = (E_ .* (d50riv .* 1e-3) .^ 3 ./ 1e-12) .^ 0.33
-    TCrbed = (
-        E_ .* d50riv .* (
-            0.13 .* E .^ (-0.392) .* exp.(-0.015 .* E .^ 2) .+
-            0.045 .* (1 .- exp.(-0.068 .* E))
+    TCrbed = @. Float(
+        E_ * d50riv * (
+            0.13 * E ^ (-0.392) * exp(-0.015 * E ^ 2) +
+            0.045 * (1 - exp(-0.068 * E))
         )
     )
     TCrbank = TCrbed
     # kd from Hanson & Simon 2001
-    kdbank = 0.2 .* TCrbank .^ (-0.5) .* 1e-6
-    kdbed = 0.2 .* TCrbed .^ (-0.5) .* 1e-6
+    kdbank = @. Float(0.2 * TCrbank ^ (-0.5) * 1e-6)
+    kdbed = @. Float(0.2 * TCrbed ^ (-0.5) * 1e-6)
 
     rs = RiverSediment(
         n = nriv,
-        Δt = Float64(Δt.value),
+        Δt = Float(Δt.value),
         # Parameters
         sl = riverslope,
         dl = riverlength,
@@ -1135,44 +1135,44 @@ function initialize_riversed(nc, config, riverwidth, riverlength, inds_riv)
         h_riv = fill(mv, nriv),
         q_riv = fill(mv, nriv),
         # Input from land
-        inlandclay = fill(0.0, nriv),
-        inlandsilt = fill(0.0, nriv),
-        inlandsand = fill(0.0, nriv),
-        inlandsagg = fill(0.0, nriv),
-        inlandlagg = fill(0.0, nriv),
-        inlandsed = fill(0.0, nriv),
+        inlandclay = zeros(Float, nriv),
+        inlandsilt = zeros(Float, nriv),
+        inlandsand = zeros(Float, nriv),
+        inlandsagg = zeros(Float, nriv),
+        inlandlagg = zeros(Float, nriv),
+        inlandsed = zeros(Float, nriv),
         # States
-        sedload = fill(0.0, nriv),
-        clayload = fill(0.0, nriv),
-        siltload = fill(0.0, nriv),
-        sandload = fill(0.0, nriv),
-        saggload = fill(0.0, nriv),
-        laggload = fill(0.0, nriv),
-        gravload = fill(0.0, nriv),
-        sedstore = fill(0.0, nriv),
-        claystore = fill(0.0, nriv),
-        siltstore = fill(0.0, nriv),
-        sandstore = fill(0.0, nriv),
-        saggstore = fill(0.0, nriv),
-        laggstore = fill(0.0, nriv),
-        gravstore = fill(0.0, nriv),
-        outsed = fill(0.0, nriv),
-        outclay = fill(0.0, nriv),
-        outsilt = fill(0.0, nriv),
-        outsand = fill(0.0, nriv),
-        outsagg = fill(0.0, nriv),
-        outlagg = fill(0.0, nriv),
-        outgrav = fill(0.0, nriv),
+        sedload = zeros(Float, nriv),
+        clayload = zeros(Float, nriv),
+        siltload = zeros(Float, nriv),
+        sandload = zeros(Float, nriv),
+        saggload = zeros(Float, nriv),
+        laggload = zeros(Float, nriv),
+        gravload = zeros(Float, nriv),
+        sedstore = zeros(Float, nriv),
+        claystore = zeros(Float, nriv),
+        siltstore = zeros(Float, nriv),
+        sandstore = zeros(Float, nriv),
+        saggstore = zeros(Float, nriv),
+        laggstore = zeros(Float, nriv),
+        gravstore = zeros(Float, nriv),
+        outsed = zeros(Float, nriv),
+        outclay = zeros(Float, nriv),
+        outsilt = zeros(Float, nriv),
+        outsand = zeros(Float, nriv),
+        outsagg = zeros(Float, nriv),
+        outlagg = zeros(Float, nriv),
+        outgrav = zeros(Float, nriv),
         # Outputs
-        Sedconc = fill(0.0, nriv),
-        SSconc = fill(0.0, nriv),
-        Bedconc = fill(0.0, nriv),
-        maxsed = fill(0.0, nriv),
-        erodsed = fill(0.0, nriv),
-        erodsedbank = fill(0.0, nriv),
-        erodsedbed = fill(0.0, nriv),
-        depsed = fill(0.0, nriv),
-        insed = fill(0.0, nriv),
+        Sedconc = zeros(Float, nriv),
+        SSconc = zeros(Float, nriv),
+        Bedconc = zeros(Float, nriv),
+        maxsed = zeros(Float, nriv),
+        erodsed = zeros(Float, nriv),
+        erodsedbank = zeros(Float, nriv),
+        erodsedbed = zeros(Float, nriv),
+        depsed = zeros(Float, nriv),
+        insed = zeros(Float, nriv),
         # Reservoir / lake
         wbcover = wbcover,
         wblocs = wblocs,
diff --git a/src/sediment_model.jl b/src/sediment_model.jl
index 74182db0b..f8e17e69c 100644
--- a/src/sediment_model.jl
+++ b/src/sediment_model.jl
@@ -38,10 +38,10 @@ function initialize_sediment_model(config::Config)
     river_2d = ncread(nc, param(config, "input.river_location"); type = Bool, fill = false)
     river = river_2d[inds]
     riverwidth_2d =
-        ncread(nc, param(config, "input.lateral.river.width"); type = Float64, fill = 0)
+        ncread(nc, param(config, "input.lateral.river.width"); type = Float, fill = 0)
     riverwidth = riverwidth_2d[inds]
     riverlength_2d =
-        ncread(nc, param(config, "input.lateral.river.length"); type = Float64, fill = 0)
+        ncread(nc, param(config, "input.lateral.river.length"); type = Float, fill = 0)
     riverlength = riverlength_2d[inds]
 
     inds_riv, rev_inds_riv = active_indices(river_2d, 0)
@@ -79,7 +79,7 @@ function initialize_sediment_model(config::Config)
 
     # # lateral part sediment in overland flow
     rivcell = float(river)
-    ols = OverlandFlowSediment{Float64}(
+    ols = OverlandFlowSediment{Float}(
         n = n,
         rivcell = rivcell,
         h_riv = fill(mv, n),
@@ -116,7 +116,7 @@ function initialize_sediment_model(config::Config)
     # River processes
 
     # the indices of the river cells in the land(+river) cell vector
-    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float64)
+    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float)
     clamp!(βₗ, 0.00001, Inf)
     riverlength = riverlength_2d[inds_riv]
     riverwidth = riverwidth_2d[inds_riv]
@@ -168,7 +168,7 @@ function initialize_sediment_model(config::Config)
     if reinit == false
         instate_path = joinpath(tomldir, config.state.path_input)
         state_ncnames = ncnames(config.state)
-        set_states(instate_path, model, state_ncnames; type = Float64)
+        set_states(instate_path, model, state_ncnames; type = Float)
     end
 
     # make sure the forcing is already loaded
diff --git a/src/utils.jl b/src/utils.jl
index a328abb44..20a639c77 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -12,7 +12,7 @@ const pcrdir = [
     CartesianIndex(-1, 1),  # 9
 ]
 
-const mv = NaN
+const mv = Float(NaN)
 
 # timestep that the parameter units are defined in
 const basetimestep = Second(Day(1))
@@ -77,7 +77,7 @@ function active_indices(network, key::Tuple)
     end
 end
 
-function lattometres(lat::Float64)
+function lattometres(lat::Real)
     m1 = 111132.92     # latitude calculation term 1
     m2 = -559.82       # latitude calculation term 2
     m3 = 1.175         # latitude calculation term 3
@@ -176,6 +176,9 @@ function ncread(
 
     if isnothing(var)
         @assert !isnothing(defaults)
+        if !isnothing(type)
+            defaults = convert(type, defaults)
+        end
         if isnothing(dimname)
             return Base.fill(defaults, length(sel))
         else
@@ -236,13 +239,13 @@ function ncread(
 end
 
 """
-    set_layerthickness(d::Float64, sl::SVector)
+    set_layerthickness(d::Real, sl::SVector)
 
 Calculate actual soil thickness of layers based on a reference depth (e.g. soil depth or water table depth) `d`,
 a SVector `sl` with cumulative soil depth starting at soil surface (0), and a SVector `tl` with actual thickness
 per soil layer.
 """
-function set_layerthickness(d::Float64, sl::SVector, tl::SVector)
+function set_layerthickness(d::Real, sl::SVector, tl::SVector)
 
     act_d = tl .* mv
     for i = 1:length(act_d)
diff --git a/src/vertical_process.jl b/src/vertical_process.jl
index bcba2b925..dc1a67ae9 100644
--- a/src/vertical_process.jl
+++ b/src/vertical_process.jl
@@ -369,7 +369,7 @@ function snowpack_hbv(
 
     # fraction of precipitation which falls as rain
     rainfrac = if iszero(tti)
-        Float64(temperature > tt)
+        Float(temperature > tt)
     else
         frac = (temperature - (tt - tti / 2.0)) / tti
         min(frac, 1.0)


Test patch:
diff --git a/test/bmi.jl b/test/bmi.jl
index e30e73e9b..89a0825ad 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -32,12 +32,12 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
         @test_throws ErrorException BMI.get_var_grid(model, "lateral.river.lake.volume")
         # Vector{Float64} printing on Julia 1.6+
         @test BMI.get_var_type(model, "lateral.river.reservoir.inflow") in
-              ["Array{Float64,1}", "Vector{Float64}"]
+              ["Array{$Float,1}", "Vector{$Float}"]
         @test BMI.get_var_type(model, "vertical.n") == "Int64"
         @test BMI.get_var_units(model, "vertical.θₛ") == "mm mm-1"
-        @test BMI.get_var_itemsize(model, "lateral.subsurface.ssf") == 8
+        @test BMI.get_var_itemsize(model, "lateral.subsurface.ssf") == sizeof(Float)
         @test BMI.get_var_nbytes(model, "vertical.n") == 8
-        @test BMI.get_var_nbytes(model, "lateral.river.q") == 45240
+        @test BMI.get_var_nbytes(model, "lateral.river.q") == 5655 * sizeof(Float)
         @test BMI.get_var_location(model, "lateral.river.q") == "node"
     end
 
diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl
index b4a4ba9b7..59280900d 100644
--- a/test/reservoir_lake.jl
+++ b/test/reservoir_lake.jl
@@ -62,9 +62,9 @@ sh = [
 ]
 @testset "linked lakes (HBV)" begin
     @test keys(sh[1]) == (:H, :S)
-    @test typeof(values(sh[1])) == Tuple{Vector{Float64},Vector{Float64}}
+    @test typeof(values(sh[1])) == Tuple{Vector{Float},Vector{Float}}
 
-    lake = Wflow.NaturalLake{Float64}(
+    lake = Wflow.NaturalLake{Float}(
         lowerlake_ind = [2, 0],
         area = [472461536.0, 60851088.0],
         threshold = [393.7, 0.0],
@@ -90,17 +90,17 @@ sh = [
 
     Wflow.update(lake, 1, 500.0, 15, 86400.0)
     Wflow.update(lake, 2, 500.0, 15, 86400.0)
-    @test lake.outflow ≈ [214.80170846121263, 236.83281600000214]
-    @test lake.totaloutflow ≈ [1.855886761104877e7, 2.0462355302400187e7]
-    @test lake.storage ≈ [1.2737435094769483e9, 2.6019755340159863e8]
-    @test lake.waterlevel ≈ [395.0912274997361, 395.2101079057371]
+    @test lake.outflow ≈ [214.80170846121263, 236.83281600000214] atol=1e-2
+    @test lake.totaloutflow ≈ [1.855886761104877e7, 2.0462355302400187e7] atol=1e3
+    @test lake.storage ≈ [1.2737435094769483e9, 2.6019755340159863e8] atol=1e4
+    @test lake.waterlevel ≈ [395.0912274997361, 395.2101079057371] atol=1e-2
     lake.totaloutflow .= 0.0
     lake.inflow .= 0.0
     Wflow.update(lake, 1, 500.0, 15, 86400.0)
     Wflow.update(lake, 2, 500.0, 15, 86400.0)
-    @test lake.outflow ≈ [0.0, 239.66710359986183]
-    @test lake.totaloutflow ≈ [-2.2446764487487033e7, 2.070723775102806e7]
-    @test lake.storage ≈ [1.3431699662524352e9, 2.6073035986708355e8]
-    @test lake.waterlevel ≈ [395.239782021054, 395.21771942667266]
+    @test lake.outflow ≈ [0.0, 239.66710359986183] atol=1e-2
+    @test lake.totaloutflow ≈ [-2.2446764487487033e7, 2.070723775102806e7] atol=1e3
+    @test lake.storage ≈ [1.3431699662524352e9, 2.6073035986708355e8] atol=1e4
+    @test lake.waterlevel ≈ [395.239782021054, 395.21771942667266] atol=1e-2
 
 end
diff --git a/test/run_hbv.jl b/test/run_hbv.jl
index 26735695c..1db8db33f 100644
--- a/test/run_hbv.jl
+++ b/test/run_hbv.jl
@@ -13,23 +13,23 @@ flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
     row = csv_first_row(model.writer.csv_path)
 
     @test row.time == DateTime("2000-01-01T00:00:00")
-    @test row.Q ≈ 517.9761666607884
-    @test row.temp_bycoord ≈ 2.965437173843384
-    @test row.temp_byindex ≈ 1.1716821193695068
-    @test row.Q_1 ≈ 501.6252231440609
-    @test row.perc_33 ≈ 2.308000087738037
-    @test row.perc_34 ≈ 1.8980000019073486
-    @test row.perc_35 ≈ 2.7100000381469727
-    @test row.perc_36 ≈ 3.818000078201294
-    @test row.perc_37 ≈ 2.1440000534057617
+    @test row.Q ≈ 517.9761666607884f0
+    @test row.temp_bycoord ≈ 2.965437173843384f0
+    @test row.temp_byindex ≈ 1.1716821193695068f0
+    @test row.Q_1 ≈ 501.6252231440609f0
+    @test row.perc_33 ≈ 2.308000087738037f0
+    @test row.perc_34 ≈ 1.8980000019073486f0
+    @test row.perc_35 ≈ 2.7100000381469727f0
+    @test row.perc_36 ≈ 3.818000078201294f0
+    @test row.perc_37 ≈ 2.1440000534057617f0
 end
 
 @testset "first timestep" begin
     hbv = model.vertical
     @test hbv.tt[1] ≈ 0.0
     @test model.clock.iteration == 2
-    @test hbv.soilmoisture[1] == 134.35299682617188
-    @test hbv.runoff[1] == 7.406898120121746
+    @test hbv.soilmoisture[1] ≈ 134.35299682617188f0
+    @test hbv.runoff[1] ≈ 7.406898120121746f0
     @test hbv.soilevap[1] == 0.0
     @test hbv.snow[1] == 0.0
 end
@@ -39,26 +39,26 @@ model = Wflow.update(model)
 
 @testset "second timestep" begin
     hbv = model.vertical
-    @test hbv.soilmoisture[1] == 134.35299682617188
-    @test hbv.runoff[1] == 4.353345933656929
+    @test hbv.soilmoisture[1] ≈ 134.35299682617188f0
+    @test hbv.runoff[1] ≈ 4.3533463f0
     @test hbv.soilevap[1] == 0.0
     @test hbv.snow[1] == 0.0
 end
 
 @testset "overland flow" begin
     q = model.lateral.land.q_av
-    @test sum(q) ≈ 3278.9423039643552
-    @test q[500] ≈ 0.2348646912841589
-    @test q[1566] ≈ 0.030463805623037462
-    @test q[network.land.order[end]] ≈ 0.296794455575309
+    @test sum(q) ≈ 3278.9423039643552f0
+    @test q[500] ≈ 0.2348646912841589f0
+    @test q[1566] ≈ 0.030463805623037462f0
+    @test q[network.land.order[end]] ≈ 0.296794455575309f0
 end
 
 @testset "river flow" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 51943.24072462076
-    @test q[500] ≈ 5.693405508134953
-    @test q[88] ≈ 8.839306661545436
-    @test q[network.river.order[end]] ≈ 363.5647479049738
+    @test sum(q) ≈ 51943.24072462076f0
+    @test q[500] ≈ 5.693405508134953f0
+    @test q[88] ≈ 8.839306661545436f0
+    @test q[network.river.order[end]] ≈ 363.5647479049738f0
 end
 
 benchmark = @benchmark Wflow.update(model)
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index d61ad2311..81e201f59 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -14,73 +14,73 @@ flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
     row = csv_first_row(model.writer.csv_path)
 
     @test row.time == DateTime("2000-01-01T00:00:00")
-    @test row.Q ≈ 7.815587720999557
-    @test row.volume ≈ 4.364251626782245e7
-    @test row.temp_bycoord ≈ 2.3279826641082764
-    @test row.temp_byindex ≈ 2.3279826641082764
-    @test row.Q_6336050 ≈ 0.02388490515704142
-    @test row.Q_6336510 ≈ 0.012411069754330572
-    @test row.Q_6836100 ≈ 0.004848609260578462
-    @test row.Q_6336500 ≈ 0.011474801875419624
-    @test row.Q_6836190 ≈ 0.0005262529349807199
-    @test row.Q_6336800 ≈ 0.013467698526109831
-    @test row.Q_6336900 ≈ 0.003408280699908697
-    @test row.Q_6336930 ≈ 0.09773275295006544
-    @test row.Q_6336910 ≈ 0.002147610203630289
-    @test row.Q_6336920 ≈ 0.002649393436607859
-    @test row.Q_6136100 ≈ 0.0008708128761381345
-    @test row.Q_6136500 ≈ 0.000729148906480041
-    @test row.Q_6136520 ≈ 0.002155395279410574
-    @test row.Q_6136150 ≈ 0.0022298329229975987
-    @test row.Q_6136151 ≈ 0.0031045027476524524
-    @test row.Q_6136160 ≈ 3.3423894540713786
-    @test row.Q_6136200 ≈ 1.358270076420503
-    @test row.Q_6136201 ≈ 5.942330218662147
-    @test row.Q_6136202 ≈ 1.680997341533662
-    @test row.recharge_1 ≈ -0.027398093386017383
+    @test row.Q ≈ 7.815587720999557f0
+    @test row.volume ≈ 4.364251626782245f7
+    @test row.temp_bycoord ≈ 2.3279826641082764f0
+    @test row.temp_byindex ≈ 2.3279826641082764f0
+    @test row.Q_6336050 ≈ 0.02388490515704142f0
+    @test row.Q_6336510 ≈ 0.012411069754330572f0
+    @test row.Q_6836100 ≈ 0.004848609260578462f0
+    @test row.Q_6336500 ≈ 0.011474801875419624f0
+    @test row.Q_6836190 ≈ 0.0005262529349807199f0
+    @test row.Q_6336800 ≈ 0.013467698526109831f0
+    @test row.Q_6336900 ≈ 0.003408280699908697f0
+    @test row.Q_6336930 ≈ 0.09773275295006544f0
+    @test row.Q_6336910 ≈ 0.002147610203630289f0
+    @test row.Q_6336920 ≈ 0.002649393436607859f0
+    @test row.Q_6136100 ≈ 0.0008708128761381345f0
+    @test row.Q_6136500 ≈ 0.000729148906480041f0
+    @test row.Q_6136520 ≈ 0.002155395279410574f0
+    @test row.Q_6136150 ≈ 0.0022298329229975987f0
+    @test row.Q_6136151 ≈ 0.0031045027476524524f0
+    @test row.Q_6136160 ≈ 3.3423894540713786f0
+    @test row.Q_6136200 ≈ 1.358270076420503f0
+    @test row.Q_6136201 ≈ 5.942330218662147f0
+    @test row.Q_6136202 ≈ 1.680997341533662f0
+    @test row.recharge_1 ≈ -0.027398093386017383f0
 end
 
 @testset "NetCDF scalar output" begin
     ds = model.writer.dataset_scalar
     @test ds["time"][1] == DateTime("2000-01-01T00:00:00")
     @test ds["Q"][:] ≈ [
-        0.023884907,
-        0.012411069,
-        0.004848609,
-        0.011474802,
-        0.0005262529,
-        0.013467698,
-        0.003408281,
-        0.09773275,
-        0.0021476103,
-        0.0026493936,
-        0.0008708129,
-        0.0007291489,
-        0.0021553952,
-        0.002229833,
-        0.003104503,
-        3.3423893,
-        1.35827,
-        5.9423304,
-        1.6809973,
+        0.023884907f0,
+        0.012411069f0,
+        0.004848609f0,
+        0.011474802f0,
+        0.0005262529f0,
+        0.013467698f0,
+        0.003408281f0,
+        0.09773275f0,
+        0.0021476103f0,
+        0.0026493936f0,
+        0.0008708129f0,
+        0.0007291489f0,
+        0.0021553952f0,
+        0.002229833f0,
+        0.003104503f0,
+        3.3423893f0,
+        1.35827f0,
+        5.9423304f0,
+        1.6809973f0,
     ]
     @test ds["Q_gauges"].attrib["cf_role"] == "timeseries_id"
-    @test ds["temp_index"][:] ≈ [2.3279827]
-    @test ds["temp_coord"][:] ≈ [2.3279827]
+    @test ds["temp_index"][:] ≈ [2.3279827f0]
+    @test ds["temp_coord"][:] ≈ [2.3279827f0]
     @test keys(ds.dim) == ["time", "Q_gauges", "temp_bycoord", "temp_byindex"]
 end
 
 @testset "first timestep" begin
     sbm = model.vertical
 
-    @test sbm.tt[1] ≈ 1.2999999523162842
+    @test sbm.tt[1] ≈ 1.2999999523162842f0
 
     @test model.clock.iteration == 2
 
-    @test sbm.θₛ[1] ≈ 0.48343977332115173
+    @test sbm.θₛ[1] ≈ 0.48343977332115173f0
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] == 0.0
-    @test sbm.snow[1] ≈ 0.6029989752244306
+    @test sbm.snow[1] ≈ 0.6029989752244306f0
 end
 
 # run the second timestep
@@ -88,43 +88,43 @@ model = Wflow.update(model)
 
 @testset "second timestep" begin
     sbm = model.vertical
-    @test sbm.θₛ[1] ≈ 0.48343977332115173
+    @test sbm.θₛ[1] ≈ 0.48343977332115173f0
     @test sbm.runoff[1] == 0.0
-    @test sbm.soilevap[1] ≈ 0.005865651540305367
-    @test sbm.snow[1] ≈ 0.009696763863612956
+    @test sbm.soilevap[1] ≈ 0.005865651540305367f0
+    @test sbm.snow[1] ≈ 0.009696763863612956f0
 end
 
 @testset "subsurface flow" begin
     ssf = model.lateral.subsurface.ssf
-    @test sum(ssf) ≈ 6.368140761295825e16
-    @test ssf[network.land.order[1]] ≈ 3.0449782003445332e13
-    @test ssf[network.land.order[end-100]] ≈ 7.855716879739626e11
-    @test ssf[network.land.order[end]] ≈ 2.161246841709492e11
+    @test sum(ssf) ≈ 6.368140761295825f16
+    @test ssf[network.land.order[1]] ≈ 3.0449782003445332f13
+    @test ssf[network.land.order[end-100]] ≈ 7.855716879739626f11
+    @test ssf[network.land.order[end]] ≈ 2.161246841709492f11
 end
 
 @testset "overland flow" begin
     q = model.lateral.land.q_av
-    @test sum(q) ≈ 319.69538082148307
+    @test sum(q) ≈ 319.69538082148307f0
     @test q[26625] ≈ 0.0
     @test q[39308] ≈ 0.0
-    @test q[network.land.order[end]] ≈ 1.0978547520221912e-5
+    @test q[network.land.order[end]] ≈ 1.0978547520221912f-5
 end
 
 @testset "river flow" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 2807.884971209105
-    @test q[4061] ≈ 0.0016288040314320486
-    @test q[5617] ≈ 7.338169165884175
-    @test q[network.river.order[end]] ≈ 0.00610520650626283
+    @test sum(q) ≈ 2807.884971209105f0
+    @test q[4061] ≈ 0.0016288040314320486f0
+    @test q[5617] ≈ 7.338169165884175f0
+    @test q[network.river.order[end]] ≈ 0.00610520650626283f0
 end
 
 @testset "reservoir simple" begin
     res = model.lateral.river.reservoir
-    @test res.outflow[2] ≈ 0.2174998592483153
-    @test res.inflow[2] ≈ 50.170880189190626
-    @test res.volume[2] ≈ 2.776162917050312e7
-    @test res.precipitation[2] ≈ 0.1765228509902954
-    @test res.evaporation[2] ≈ 0.5372688174247742
+    @test res.outflow[2] ≈ 0.2174998592483153f0
+    @test res.inflow[2] ≈ 50.170880189190626f0
+    @test res.volume[2] ≈ 2.776162917050312f7
+    @test res.precipitation[2] ≈ 0.1765228509902954f0
+    @test res.evaporation[2] ≈ 0.5372688174247742f0
 end
 
 # set these variables for comparison in "changed dynamic parameters"
@@ -161,8 +161,8 @@ end
 # downstream from pit is at river index 1739, CartesianIndex(142, 85)
 @testset "river flow at and downstream of pit" begin
     q = model.lateral.river.q_av
-    @test q[1765] ≈ 8.133435566601927
-    @test q[1739] ≈ 0.008996045895155833
+    @test q[1765] ≈ 8.133435566601927f0
+    @test q[1739] ≈ 0.008996045895155833f0
 end
 
 # test changing forcing and cyclic LAI parameter
@@ -186,9 +186,9 @@ model = Wflow.update(model)
 @testset "changed dynamic parameters" begin
     res = model.lateral.river.reservoir
     vertical = model.vertical
-    @test vertical.precipitation[2] / precip[2] ≈ 2.0
-    @test (vertical.potential_evaporation[100] - 1.50) / evap[100] ≈ 3.0
-    @test vertical.leaf_area_index[100] / lai[100] ≈ 1.6
-    @test (res.evaporation[2] - 1.50) / res_evap[2] ≈ 3.0000012203408635
+    @test vertical.precipitation[2] / precip[2] ≈ 2.0f0
+    @test (vertical.potential_evaporation[100] - 1.50) / evap[100] ≈ 3.0f0
+    @test vertical.leaf_area_index[100] / lai[100] ≈ 1.6f0
+    @test (res.evaporation[2] - 1.50) / res_evap[2] ≈ 3.0000012203408635f0
 end
 Wflow.close_files(model, delete_output = false)
diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl
index 4596cbf99..53896624b 100644
--- a/test/run_sbm_gwf.jl
+++ b/test/run_sbm_gwf.jl
@@ -13,18 +13,18 @@ flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
     row = csv_first_row(model.writer.csv_path)
 
     @test row.time == DateTime("2000-06-01T00:00:00")
-    @test row.Q_av ≈ 0.016260519762890745
-    @test row.head ≈ 1.8416896245327632
+    @test row.Q_av ≈ 0.016260519762890745f0
+    @test row.head ≈ 1.8416896245327632f0
 end
 
 @testset "first timestep" begin
     sbm = model.vertical
 
     @test model.clock.iteration == 2
-    @test sbm.θₛ[1] == 0.44999998807907104
+    @test sbm.θₛ[1] ≈ 0.44999998807907104f0
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] == 0.0
-    @test sbm.transpiration[1] == 0.6117526566330049
+    @test sbm.transpiration[1] ≈ 0.6117526566330049f0
 end
 
 # run the second timestep
@@ -32,10 +32,10 @@ model = Wflow.update_sbm_gwf(model)
 
 @testset "second timestep" begin
     sbm = model.vertical
-    @test sbm.θₛ[1] == 0.44999998807907104
+    @test sbm.θₛ[1] ≈ 0.44999998807907104f0
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] == 0.0
-    @test sbm.transpiration[1] == 1.0122634204681036
+    @test sbm.transpiration[1] ≈ 1.0122634204681036f0
 end
 
 @testset "overland flow" begin
@@ -45,19 +45,19 @@ end
 
 @testset "river flow" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 0.03515257228971982
-    @test q[6] ≈ 0.007952670041402076
-    @test q[13] ≈ 0.0005980549969548235
-    @test q[network.river.order[end]] ≈ 0.008482648782941154
+    @test sum(q) ≈ 0.03515257228971982f0
+    @test q[6] ≈ 0.007952670041402076f0
+    @test q[13] ≈ 0.0005980549969548235f0
+    @test q[network.river.order[end]] ≈ 0.008482648782941154f0
 end
 
 @testset "groundwater" begin
     gw = model.lateral.subsurface
-    @test gw.river.stage[1] ≈ 1.21057153442702
-    @test gw.flow.aquifer.head[19] ≈ 1.7999999523162842
-    @test gw.river.flux[1] ≈ -50.751898489778426
+    @test gw.river.stage[1] ≈ 1.21057153442702f0
+    @test gw.flow.aquifer.head[19] ≈ 1.7999999523162842f0
+    @test gw.river.flux[1] ≈ -50.751898489778426f0
     @test gw.drain.flux[1] ≈ 0.0
-    @test gw.recharge.rate[19] ≈ -0.0014241196552847502
+    @test gw.recharge.rate[19] ≈ -0.0014241196552847502f0
 end
 
 Wflow.close_files(model)
diff --git a/test/run_sediment.jl b/test/run_sediment.jl
index b043c3ff5..eae48b0fc 100644
--- a/test/run_sediment.jl
+++ b/test/run_sediment.jl
@@ -11,13 +11,13 @@ model = Wflow.update(model)
 @testset "first timestep sediment model (vertical)" begin
     eros = model.vertical
 
-    @test eros.erosov[1] == 0.9
+    @test eros.erosov[1] ≈ 0.9f0
     @test model.clock.iteration == 2
-    @test mean(eros.leaf_area_index) ≈ 1.7120018886212223
-    @test eros.dmsand[1] == 200.0
-    @test eros.dmlagg[1] == 500.0
-    @test mean(eros.interception) ≈ 0.4767846753916875
-    @test mean(eros.soilloss) ≈ 0.00846987746234039
+    @test mean(eros.leaf_area_index) ≈ 1.7120018886212223f0
+    @test eros.dmsand[1] == 200.0f0
+    @test eros.dmlagg[1] == 500.0f0
+    @test mean(eros.interception) ≈ 0.4767846753916875f0
+    @test mean(eros.soilloss) ≈ 0.00846987746234039f0
 end
 
 # run the second timestep
@@ -26,29 +26,29 @@ model = Wflow.update(model)
 @testset "second timestep sediment model (vertical)" begin
     eros = model.vertical
 
-    @test mean(eros.soilloss) ≈ 0.07441259472436522
-    @test mean(eros.erosclay) ≈ 0.002198312957830507
-    @test mean(eros.erossand) ≈ 0.02515898751801652
-    @test mean(eros.eroslagg) ≈ 0.021613617766119256
+    @test mean(eros.soilloss) ≈ 0.07441259472436522f0
+    @test mean(eros.erosclay) ≈ 0.002198312957830507f0
+    @test mean(eros.erossand) ≈ 0.02515898751801652f0
+    @test mean(eros.eroslagg) ≈ 0.021613617766119256f0
     @test mean(eros.TCsed) == 0.0
-    @test mean(eros.TCsilt) ≈ 1.0988158364353527e6
-    @test mean(eros.TCsand) ≈ 1.0987090622888755e6
-    @test mean(eros.TCclay) ≈ 1.0992655197016734e6
-    @test mean(eros.TCsilt) ≈ 1.0988158364353527e6
+    @test mean(eros.TCsilt) ≈ 1.0988158364353527f6
+    @test mean(eros.TCsand) ≈ 1.0987090622888755f6
+    @test mean(eros.TCclay) ≈ 1.0992655197016734f6
+    @test mean(eros.TCsilt) ≈ 1.0988158364353527f6
 end
 
 @testset "second timestep sediment model (lateral)" begin
     lat = model.lateral
 
-    @test mean(lat.land.inlandsed) ≈ 0.0727443455767198
-    @test mean(lat.land.inlandclay) ≈ 0.002190037611542877
-    @test mean(lat.land.inlandsand) ≈ 0.024512058443074893
-    @test mean(lat.land.olclay) ≈ 0.006212441628146587
+    @test mean(lat.land.inlandsed) ≈ 0.0727443455767198f0
+    @test mean(lat.land.inlandclay) ≈ 0.002190037611542877f0
+    @test mean(lat.land.inlandsand) ≈ 0.024512058443074893f0
+    @test mean(lat.land.olclay) ≈ 0.006212441628146587f0
 
-    @test mean(lat.river.SSconc) ≈ 0.8454728220852514
-    @test mean(lat.river.inlandclay) ≈ 0.019393435838709512
-    @test lat.river.h_riv[network.river.order[end]] ≈ 0.006103649735450745
-    @test lat.river.outclay[1] ≈ 3.3447705344961993e-6
+    @test mean(lat.river.SSconc) ≈ 0.8454728220852514f0
+    @test mean(lat.river.inlandclay) ≈ 0.019393435838709512f0
+    @test lat.river.h_riv[network.river.order[end]] ≈ 0.006103649735450745f0
+    @test lat.river.outclay[1] ≈ 3.3447705344961993f-6
 end
 
 Wflow.close_files(model)
diff --git a/test/runtests.jl b/test/runtests.jl
index 10e5b97ba..046b9b927 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,6 +15,7 @@ using DelimitedFiles
 using Logging
 
 const BMI = BasicModelInterface
+const Float = Wflow.Float
 
 # ensure test data is present
 testdir = @__DIR__