src/Callbacks/courant_step.jl

 struct CFLHeatState{A}
     Δt::A
+    default_dt::Float64
 end
-function (fn::CryoGridCallbackFunction{<:CFLHeatState})(u,p,t)
+function (fn::CryoGridCallbackFunction{<:CFLHeatState{<:AbstractArray}})(u,p,t)
     let Δt = fn.state.Δt,
-        Δx = Δ(fn.setup.grid),
+        Δx = dustrip(Δ(fn.setup.grid)),
         Ceff = getvar(:Ceff, fn.setup, u), # apparent heat capacity
         kc = getvar(:kc, fn.setup, u); # thermal cond. at grid centers
-        @. Δt = Δx^2 * Ceff / kc
-        minimum(Δt)
+        @. Δt = Utils.adstrip(Δx^2 * Ceff / kc)
+        Δt_min = minimum(Δt)
+        IfElse.ifelse(isfinite(Δt_min) && Δt_min > 0, Δt_min, fn.state.default_dt)
     end
 end
-function CFLStepLimiter(setup::HeatOnlySetup; courant_number=1//2)
-    state = CFLHeatState(similar(Δ(setup.grid)))
+function (fn::CryoGridCallbackFunction{<:CFLHeatState{Nothing}})(u,p,t)
+    let Δx = dustrip(Δ(fn.setup.grid)),
+        Ceff = getvar(:Ceff, fn.setup, u), # apparent heat capacity
+        kc = getvar(:kc, fn.setup, u); # thermal cond. at grid centers
+        Δt = Utils.adstrip(Δx^2 * Ceff / kc)
+        Δt_min = minimum(Δt)
+        IfElse.ifelse(isfinite(Δt_min) && Δt_min > 0, Δt_min, fn.state.default_dt)
+    end
+end
+function CFLStepLimiter(setup::HeatOnlyLandModel; courant_number=1//2, default_dt=60.0, iip=true)
+    state = iip ? CFLHeatState(zero(dustrip(Δ(setup.grid))), default_dt) : CFLHeatState(nothing, default_dt)
     fn = CryoGridCallbackFunction(setup, state)
     StepsizeLimiter(fn; safety_factor=courant_number, max_step=true)
 end

src/Common/Forcings/Forcings.jl

-module Forcings
-
-using Dates
-using Interpolations
-using TimeSeries
-using Unitful
-
-export Forcing, TimeSeriesForcing
-
-"""
-      Forcing{T,N}
-"""
-abstract type Forcing{T,N} end
-
-include("timeseries.jl")
-
-end

src/Common/Numerics/variables.jl

-# Variable dimensions
-abstract type VarDim end
-struct OnGrid{S} <: VarDim
-    f::Function # G -> G' where G is the edge grid
-    OnGrid(::Type{S}, f::Function=x->x) where {S<:GridSpec} = new{S}(f)
-end
-struct Shape{S} <: VarDim Shape(dims::Int...) = new{dims}() end
-const Scalar = Shape()
-
-abstract type Var{name,T,S<:Union{<:VarDim,<:GridSpec}} end
-struct Prognostic{name,T,S} <: Var{name,T,S}
-    dim::S
-    Prognostic(name::Symbol, ::Type{T}, dims::Union{<:VarDim,<:GridSpec}) where {T} = new{name, T, typeof(dims)}(dims)
-end
-struct Algebraic{name,T,S} <: Var{name,T,S}
-    dim::S
-    # maybe a mass matrix init function?
-    Algebraic(name::Symbol, ::Type{T}, dims::Union{<:VarDim,<:GridSpec}) where {T} = new{name, T, typeof(dims)}(dims)
-end
-struct Diagnostic{name,T,S} <: Var{name,T,S}
-    dim::S
-    Diagnostic(name::Symbol, ::Type{T}, dims::Union{<:VarDim,<:GridSpec}) where {T} = new{name, T, typeof(dims)}(dims)
-end
-struct Parameter{name,T,S,D} <: Var{name,T,S}
-    dim::S
-    default_value::T
-    Parameter(name::Symbol, default_value::T, domain::Interval{L,R,I}=-Inf..Inf) where {L,R,I<:Real,T<:AbstractArray} =
-        new{name,T,typeof(Shape(size(default_value)...)),convert(Interval{L,R,Float64}, domain)}(Shape(size(default_value)...), default_value)
-    Parameter(name::Symbol, default_value::T, domain::Interval=-Inf..Inf) where {T<:Real} = Parameter(name, [default_value], domain)
-end
-==(var1::Var{N1,T1,D1},var2::Var{N2,T2,D2}) where {N1,N2,T1,T2,D1,D2} = (N1==N2) && (T1==T2) && (D1==D2)
-varname(::Var{name}) where {name} = name
-varname(::Type{<:Var{name}}) where {name} = name
-vartype(::Var{name,T}) where {name,T} = T
-vartype(::Type{<:Var{name,T}}) where {name,T} = T
-vardims(var::Var{name,T,S}) where {name,T,S} = var.dim
-isprognostic(::T) where {T<:Var} = T <: Prognostic
-isalgebraic(::T) where {T<:Var} = T <: Algebraic
-isdiagnostic(::T) where {T<:Var} = T <: Diagnostic
-isparameter(::T) where {T<:Var} = T <: Parameter
-domain(::Parameter{name,T,S,D}) where {name,T,S,D} = D
-export Var, Prognostic, Algebraic, Diagnostic, Parameter
-export VarDim, OnGrid, Shape, Scalar
-export varname, vartype, isprognostic, isalgebraic, isdiagnostic, isparameter, domain
-
-# parameter constraints
-function checkdomain(domain::Interval, f::Function, z)
-    @assert all(z .∈ 0..Inf) "value $z for parameter $name is outside of the given domain: $domain"
-    f.(z)
-end
-constrain(::Parameter{name,T,S,-Inf..Inf}, x) where {name,T,S} = x
-unconstrain(::Parameter{name,T,S,-Inf..Inf}, z) where {name,T,S} = z
-constrain(::Parameter{name,T,S,0.0..1.0}, x) where {name,T,S} = logistic.(x)
-unconstrain(::Parameter{name,T,S,0.0..1.0}, z) where {name,T,S} = checkdomain(0.0..1.0, logit, z)
-constrain(::Parameter{name,T,S,0..Inf}, x) where {name,T,S} = softplus.(x)
-unconstrain(::Parameter{name,T,S,0..Inf}, z) where {name,T,S} = checkdomain(0..Inf, softplusinv, z)
-constrain(::Parameter{name,T,S,1..Inf}, x) where {name,T,S} = (plusone ∘ softplus).(x)
-unconstrain(::Parameter{name,T,S,1..Inf}, z) where {name,T,S} = checkdomain(1..Inf, softplusinv ∘ minusone, z)
-
-export constrain, unconstrain

src/CryoGrid.jl

 module CryoGrid
 
 global CRYOGRID_DEBUG = haskey(ENV,"CG_DEBUG") && ENV["CG_DEBUG"] == "true"
+function debug(debug::Bool)
+    @warn "Enabling/disabling debug mode after the module has been loaded will not update existing types and may cause errors."
+    global CRYOGRID_DEBUG = debug
+end
 
+using Base: @propagate_inbounds
 using Reexport
 
+export Layer, SubSurface, Top, Bottom
+export Process, SubSurfaceProcess, BoundaryProcess, CompoundProcess, Coupled
+export BoundaryStyle, Dirichlet, Neumann
+export AbstractParameterization, Parameterization
+export variables, initialcondition!, diagnosticstep!, prognosticstep!, interact!, boundaryflux, boundaryvalue, observe
+
+# Common types and methods
+include("types.jl")
+include("methods.jl")
+
 # Submodules
-include("Common/Interface/Interface.jl")
-include("Common/Utils/Utils.jl")
-include("Common/Numerics/Numerics.jl")
-include("Common/Forcings/Forcings.jl")
-include("IO/InputOutput.jl")
+include("Utils/Utils.jl")
+include("Numerics/Numerics.jl")
 include("Layers/Layers.jl")
-include("Processes/Processes.jl")
-include("Setup/Setup.jl")
+include("Physics/Physics.jl")
+include("Land/Land.jl")
+include("IO/InputOutput.jl")
+include("Diagnostics/Diagnostics.jl")
+include("Drivers/Drivers.jl")
 include("Callbacks/Callbacks.jl")
 
-# Re-export submodules
-@reexport using .Interface
-@reexport using .Utils
-@reexport using .Numerics
-@reexport using .Forcings
-@reexport using .InputOutput: loadforcings, InputSpec, JsonSpec
+using .Utils
+using .Numerics
+# re-export initializer function from Numerics module
+export initializer
+# Re-exported submodules
+@reexport using .Utils: convert_tspan
+@reexport using .Numerics: Grid
 @reexport using .Layers
-@reexport using .Processes
-@reexport using .Setup
+@reexport using .Physics
+@reexport using .Land
+@reexport using .Drivers
+@reexport using .InputOutput
+@reexport using .Diagnostics
 @reexport using .Callbacks
 
 # Re-exported packages
 @reexport using Dates: Dates, DateTime
-@reexport using DimensionalData: DimArray, X, Y, Z, Dim, dims
+@reexport using DimensionalData: DimArray, Z, Dim, dims
 @reexport using IntervalSets
 @reexport using Unitful
 
-# Import parameters function into top-level scope;
-# We do not export it in order to avoid naming conflicts with other packages.
-import .Setup: parameters
-
-# Include Models submodule last to allow dependence on other submodules.
-include("Models/Models.jl")
-
-const CryoGridModels = Models
-export CryoGridModels
+# Include Presets submodule last to allow dependence on other submodules.
+include("Presets/Presets.jl")
 
 end # module

src/Diagnostics/Diagnostics.jl

+module Diagnostics
+
+using CryoGrid.Numerics
+using CryoGrid.Land
+using CryoGrid.Utils
+
+using DimensionalData
+using Statistics
+using TimeSeries
+using Unitful
+
+export zero_annual_amplitude, permafrosttable, permafrostbase, thawdepth, active_layer_thickness, mean_annual_ground_temperature
+
+include("spinup.jl")
+
+export spinup
+
+_check_arr_dims(T::AbstractDimArray) = @assert hasdim(T, Z) && hasdim(T, Ti) "Tay must have depth (Z) and time (Ti) dimensions"
+"""
+    zero_annual_amplitude(T::AbstractDimArray{<:TempQuantity}; threshold=0.5u"°C")
+
+Computes annual depth of zero amplitude (where `|max - min| < threshold`) and returns the
+result for each year. Assumes `T` to have dimensions `Ti` (time) and `Z` (depth) in any order.
+"""
+function zero_annual_amplitude(T::AbstractDimArray{<:TempQuantity}; threshold=0.5u"K")
+    _check_arr_dims(T)
+    Tarr = TimeArray(dims(T, Ti).val, permutedims(T, (Ti,Z)))
+    annual_min = collapse(Tarr, year, first, minimum)
+    annual_max = collapse(Tarr, year, first, maximum)
+    z_inds = argmax(abs.(values(annual_max) .- values(annual_min)) .< threshold, dims=2)
+    rebuild(T, dims(T, Z)[[i[2] for i in z_inds]], (Ti(timestamp(annual_max)),))[:,1]
+end
+"""
+    permafrosttable(T::AbstractDimArray{<:TempQuantity})
+
+Computes depth of permafrost table for all years, i.e. the closest depth to the surface at which
+the maximum annual temperature is strictly `< 0°C`. Assumes `T` to have dimensions `Ti` (time) and
+`Z` (depth) in any order.
+"""
+function permafrosttable(T::AbstractDimArray{<:TempQuantity})
+    _check_arr_dims(T)
+    Tarr = TimeArray(dims(T, Ti).val, permutedims(T, (Ti,Z)))
+    annual_max = collapse(Tarr, year, first, maximum)
+    # find argmax (first one/"true" bit) of annual_max < 0
+    z_inds = argmax(values(annual_max) .< 0u"°C", dims=2)
+    rebuild(T, dims(T, Z)[[i[2] for i in z_inds]], (Ti(timestamp(annual_max)),))[:,1]
+end
+"""
+    permafrostbase(T::AbstractDimArray{<:TempQuantity})
+
+Computes depth of permafrost base for all years, i.e. the closest depth to the "bottom" at which
+the maximum annual temperature is strictly `< 0°C`. Assumes `T` to have dimensions `Ti` (time) and
+`Z` (depth) in any order.
+"""
+function permafrostbase(T::AbstractDimArray{<:TempQuantity})
+    _check_arr_dims(T)
+    Tarr = TimeArray(dims(T, Ti).val, permutedims(T, (Ti,Z)))
+    annual_max = collapse(Tarr, year, first, maximum)
+    # find argmax (first one/"true" bit) of annual_max < 0
+    z_inds = size(Tarr,2) .- [i[2] for i in argmax(reverse(values(annual_max) .< 0u"°C", dims=2), dims=2)]
+    rebuild(T, dims(T, Z)[[i for i in z_inds]], (Ti(timestamp(annual_max)),))[:,1]
+end
+"""
+    thawdepth(T::AbstractDimArray{<:TempQuantity})
+
+Computes thaw depth (a.k.a freezing front) at all time steps. Assumes `T` to have dimensions `Ti` (time)
+and `Z` (depth) in any order.
+"""
+function thawdepth(T::AbstractDimArray{<:TempQuantity})
+    _check_arr_dims(T)
+    T = permutedims(T, (Z,Ti))
+    z_inds = argmax((T .<= 0u"°C").data, dims=1)
+    rebuild(T, dims(T, Z)[[i[1] for i in z_inds]], (dims(T,Ti),))[1,:]
+end
+"""
+    active_layer_thickness(T::AbstractDimArray{<:TempQuantity})
+
+Computes active layer thickness annually. The active layer thickness is defined here as the maximum thaw
+depth throughout the calendar year. Assumes `T` to have dimensions `Ti` (time) and `Z` (depth) in any order.
+"""
+function active_layer_thickness(T::AbstractDimArray{<:TempQuantity})
+    _check_arr_dims(T)
+    thaw_depth = thawdepth(T)
+    Tarr = TimeArray(dims(T, Ti).val, thaw_depth)
+    annual_max = collapse(Tarr, year, first, maximum)
+    rebuild(T, values(annual_max), (Ti(timestamp(annual_max)),))
+end
+"""
+    mean_annual_ground_temperature(T::AbstractDimArray; upper_limit=0u"m", lower_limit=10u"m")
+
+Computes mean annual ground temperature between `upper_limit` and `lower_limit`. Assumes `T` to have dimensions
+`Ti` (time) and `Z` (depth) in any order.
+"""
+function mean_annual_ground_temperature(T::AbstractDimArray{<:TempQuantity}; upper_limit=0u"m", lower_limit=10u"m")
+    _check_arr_dims(T)
+    T = permutedims(ustrip.(T), (Z,Ti))[Z(Between(lower_limit, upper_limit))]
+    return mean(T, dims=1)[1,:]
+end
+"""
+    integrate(X::AbstractDimArray, grid::Grid{Edges}; upper_limit=0u"m", lower_limit=10u"m")
+
+Integrates the quantity `X` over the given `grid`, which is assumed to be spatially alligned, i.e.
+`length(grid) == length(dims(X,Z)) + 1` and `cells(grid) .≈ dims(X,Z)` are necessary preconditions.
+"""
+function integrate(X::AbstractDimArray, grid::Grid{Edges,G,<:DistQuantity}; upper_limit=0u"m", lower_limit=10u"m") where G
+    _check_arr_dims(X)
+    X = permutedims(X, (Z,Ti))
+    @assert length(grid) == size(X,1) + 1
+    @assert all(cells(grid) .≈ dims(X,Z))
+    Δx = Δ(grid)
+    X_scaled = X.*Δx
+    return sum(X_scaled[Z(Between(lower_limit, upper_limit))], dims=1)[1,:].*area(grid)
+end
+
+end
\ No newline at end of file

src/Models/Models.jl → src/Diagnostics/spinup.jl

 """
-Pre-built CryoGrid configurations for rapid prototyping.
-"""
-module Models
-
-using CryoGrid
-using CryoGrid.InputOutput: Resource
-using Statistics
-
-include("presetgrids.jl")
-
-Forcings = (
-    Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044 = Resource("Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044", "json", "https://nextcloud.awi.de/s/F98s5WEo9xMPod7/download"),
-    Samoylov_ERA_MkL3_CCSM4_long_term = Resource("Samoylov_ERA_MkL3_CCSM4_long_term", "json", "https://nextcloud.awi.de/s/RSqBtp5sPwkCf45/download"),
-)
-
-Parameters = (
-    # Faroux et al. doi:10.1109/IGARSS.2007.4422971
-    EcoCLimMap_ULC_126_72 = Resource("EcoCLimMap_ULC_126_72", "json", "https://nextcloud.awi.de/s/nWiJr5pBoqFtw7p/download")
-)
-
-"""
-    SoilLayerConfig
-
-Helper type for representing site-specific soil layer configuration (e.g. soil and temperature profile).
-"""
-struct SoilLayerConfig{TSoilProfile,TTempProfile}
-    soilprofile::TSoilProfile
-    tempprofile::TTempProfile
-end
-
-const SamoylovDefault = SoilLayerConfig(
-    # soil profile: depth => (total water, liquid water, mineral organic, porosity)
-    SoilProfile(
-        0.0u"m" => SoilProperties(χ=0.0,ϕ=0.80,θ=1.0,ω=0.75), #(θw=0.80,θm=0.05,θo=0.15,ϕ=0.80),
-        0.1u"m" => SoilProperties(χ=0.0,ϕ=0.80,θ=1.0,ω=0.25), #(θw=0.80,θm=0.15,θo=0.05,ϕ=0.80),
-        0.4u"m" => SoilProperties(χ=0.30,ϕ=0.55,θ=1.0,ω=0.25), #(θw=0.80,θm=0.15,θo=0.05,ϕ=0.55),
-        3.0u"m" => SoilProperties(χ=0.0,ϕ=0.50,θ=1.0,ω=0.0), #(θw=0.50,θm=0.50,θo=0.0,ϕ=0.50),
-        10.0u"m" => SoilProperties(χ=0.0,ϕ=0.30,θ=1.0,ω=0.0), #(θw=0.30,θm=0.70,θo=0.0,ϕ=0.30),
-    ),
-    TempProfile(
-        0.0u"m" => -1.0u"°C",
-        2.0u"m" => -1.0u"°C",
-        5.0u"m" => -3.0u"°C",
-        10.0u"m" => -6.0u"°C",
-        25.0u"m" => -9.0u"°C",
-        100.0u"m" => -9.0u"°C",
-        1000.0u"m" => 10.2u"°C",
-    )
-)
-
-export SamoylovDefault
-
-"""
-    SoilHeat([heatvar=:H], upperbc::BoundaryProcess{Heat}, soilconfig::SoilLayerConfig; grid::Grid=DefaultGrid, freezecurve::F=FreeWater()) where {F<:FreezeCurve}
-
-Builds a simple one-layer soil/heat-conduction model with the given grid and configuration. Uses the "free water" freeze curve by default,
-but this can be changed via the `freezecurve` parameter. For example, to use the van Genuchten freeze curve, set `freezecurve=SFCC(VanGenuchten())`.
-"""
-function SoilHeat(heatvar, upperbc::BoundaryProcess{Heat}, soilconfig::SoilLayerConfig;
-    grid::Grid=DefaultGrid_5cm, freezecurve::F=FreeWater(), chunk_size=nothing) where {F<:FreezeCurve}
-    strat = Stratigraphy(
-        -2.0u"m" => Top(upperbc),
-        0.0u"m" => Ground(:soil, Soil(soilconfig.soilprofile), Heat{heatvar}(soilconfig.tempprofile, freezecurve=freezecurve)),
-        1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2"))
-    )
-    CryoGridSetup(strat,grid,chunk_size=chunk_size)
-end
-SoilHeat(upperbc::BoundaryProcess{Heat}, soilconfig::SoilLayerConfig; grid::Grid=DefaultGrid_2cm, freezecurve::F=FreeWater()) where {F<:FreezeCurve} =
-    SoilHeat(:H, upperbc, soilconfig; grid=grid, freezecurve=freezecurve)
-
-    """
-    spinup(setup::CryoGridSetup, tspan::NTuple{2,DateTime}, p, tol, layername; kwargs...)
+    spinup(setup::LandModel, tspan::NTuple{2,DateTime}, p, tol, layername; kwargs...)
 
 Implements a simple, iterative spin-up procedure.
 Runs the model specified by `setup` over `tspan` until the profile mean up to `maxdepth` over the whole time span changes only within the given tolerance `tol`.
 Returns the `ODESolution` generated by the final iteration.
 """
-function spinup(setup::CryoGridSetup, tspan::NTuple{2,DateTime}, p, tol, layername; maxdepth=100u"m", maxiter=1000, saveat=3*3600.0, solver=Euler(), dt=60.0, solve_args...)
+function spinup(setup::LandModel, tspan::NTuple{2,DateTime}, p, tol, layername; maxdepth=100u"m", maxiter=1000, saveat=3*3600.0, solver=Euler(), dt=60.0, solve_args...)
     prob = CryoGridProblem(setup, tspan, p)
     @info "Running initial solve ..."
     sol = solve(prob, solver, dt=dt, saveat=saveat, solve_args...)
     out = CryoGridOutput(sol)
     H = out.vars[layername].H
-    grid = setup.meta[layername].grids.H
+    initialstate = getstate(setup, prob.u0, similar(prob.u0), prob.tspan[1])
+    grid = initialstate[layername].grids.H
     max_ind = argmin(abs.(grid*1u"m" .- maxdepth))
     dz = Δ(grid)[1:max_ind]
     H_sub = H[1:max_ind,:]
@@ -110,5 +40,3 @@ function spinup(setup::CryoGridSetup, tspan::NTuple{2,DateTime}, p, tol, layerna
     end
     return out
 end
-
-end
\ No newline at end of file

src/Drivers/Drivers.jl

+"""
+Module containing types and functions for enabling time-integration of a CryoGrid land model.
+"""
+module Drivers
+
+using CryoGrid: Land, SubSurface, CompoundProcess, variables
+using CryoGrid.InputOutput
+using CryoGrid.Numerics
+using CryoGrid.Physics: Heat
+using CryoGrid.Utils
+
+using ComponentArrays
+using Dates
+using DimensionalData
+using Flatten
+using ModelParameters
+using LinearAlgebra
+using Reexport
+using Unitful
+
+import CryoGrid.Land: LandModel, Stratigraphy, StratComponent
+
+export DefaultJac, TridiagJac, JacobianStyle, HeatOnlyLandModel
+
+"""
+    JacobianStyle
+
+Trait for indicating Jacobian sparsity of a CryoGrid ODEProblem.
+"""
+abstract type JacobianStyle end
+struct DefaultJac <: JacobianStyle end
+struct TridiagJac <: JacobianStyle end
+"""
+    JacobianStyle(::Type{<:LandModel})
+
+Can be overriden/extended to specify Jacobian structure for specific `LandModel`s.
+"""
+JacobianStyle(::Type{<:LandModel}) = DefaultJac()
+# Auto-detect Jacobian sparsity for problems with one or more heat-only layers.
+# Note: This assumes that the processes/forcings on the boundary layers do not violate the tridiagonal structure!
+# Unfortunately, the Stratigraphy type signature is a bit nasty to work with :(
+const HeatOnlyLandModel = LandModel{<:Stratigraphy{N,<:Tuple{TTop,Vararg{<:Union{<:StratComponent{<:SubSurface, <:CompoundProcess{<:Tuple{<:Heat}}},TBot}}}}} where {N,TTop,TBot}
+JacobianStyle(::Type{<:HeatOnlyLandModel}) = TridiagJac()
+
+# DiffEq/SciML driver
+export CryoGridProblem
+include("diffeq.jl")
+
+end

src/Drivers/diffeq.jl

+using DiffEqBase
+using DiffEqCallbacks
+
+@reexport using OrdinaryDiffEq
+@reexport using DiffEqBase: solve, init, ODEProblem, SciMLBase
+
+"""
+Specialized problem type for CryoGrid `ODEProblem`s.
+"""
+struct CryoGridODEProblem end
+
+"""
+    CryoGridProblem(setup::LandModel, tspan::NTuple{2,Float64}, p=nothing;kwargs...)
+
+CryoGrid specialized constructor for ODEProblem that automatically generates the initial
+condition and necessary callbacks.
+"""
+function CryoGridProblem(
+    landmodel::LandModel,
+    u0::ComponentVector,
+    tspan::NTuple{2,Float64},
+    p=nothing;
+    saveat=3600.0,
+    savevars=(),
+    save_everystep=false,
+    callback=nothing,
+    kwargs...
+)
+    # workaround for bug in DiffEqCallbacks; see https://github.com/SciML/DifferentialEquations.jl/issues/326
+    # we have to manually expand single-number `saveat` (i.e. time interval for saving) to a step-range.
+    expandtstep(tstep::Number) = tspan[1]:tstep:tspan[end]
+    expandtstep(tstep::AbstractVector) = tstep
+    getsavestate(model::LandModel, u, du) = deepcopy(Land.getvars(model.state, u, du, savevars...))
+    savefunc(u, t, integrator) = getsavestate(integrator.f.f, Land.withaxes(u, integrator.f.f), get_du(integrator))
+    pmodel = Model(landmodel)
+    p = isnothing(p) ? dustrip.(collect(pmodel[:val])) : p
+    du0 = zero(u0)
+    # set up saving callback
+    stateproto = getsavestate(landmodel, u0, du0)
+    savevals = SavedValues(Float64, typeof(stateproto))
+    savingcallback = SavingCallback(savefunc, savevals; saveat=expandtstep(saveat), save_everystep=save_everystep)
+    callbacks = isnothing(callback) ? savingcallback : CallbackSet(savingcallback, callback)
+    # note that this implicitly discards any existing saved values in the model setup's state history
+    landmodel.hist.vals = savevals
+    # set up default mass matrix
+    M_diag = similar(landmodel.state.uproto)
+    M_idxmap = ComponentArrays.indexmap(getaxes(M_diag)[1])
+    allvars = Flatten.flatten(landmodel.state.vars, Flatten.flattenable, Var)
+    progvars = map(varname, filter(isprognostic, allvars))
+    algvars = map(varname, filter(isalgebraic, allvars))
+    for name in keys(M_idxmap)
+        M_diag_var = @view M_diag[name]
+        if name ∈ progvars
+            M_diag_var .= one(eltype(M_diag))
+        elseif name ∈ algvars
+            M_diag_var .= zero(eltype(M_diag))
+        end
+    end
+    # if no algebraic variables are present, use identity matrix
+    num_algebraic = length(M_diag) - sum(M_diag)
+    M = num_algebraic > 0 ? Diagonal(M_diag) : I
+	func = odefunction(landmodel, M, u0, p, tspan; kwargs...)
+	ODEProblem(func, u0, tspan, p, CryoGridODEProblem(); callback=callbacks, kwargs...)
+end
+"""
+    CryoGridProblem(setup::LandModel, tspan::NTuple{2,DateTime}, args...;kwargs...)
+"""
+CryoGridProblem(setup::LandModel, u0::ComponentVector, tspan::NTuple{2,DateTime}, args...;kwargs...) = CryoGridProblem(setup,u0,convert_tspan(tspan),args...;kwargs...)
+
+export CryoGridProblem
+
+"""
+    odefunction(setup::LandModel, M, u0, p, tspan; kwargs...)
+
+Constructs a SciML `ODEFunction` given the model setup, mass matrix M, initial state u0, parameters p, and tspan.
+Can (and should) be overridden by users to provide customized ODEFunction configurations for specific problem setups, e.g:
+```
+model = LandModel(strat,grid)
+function CryoGrid.Setup.odefunction(::DefaultJac, setup::typeof(model), M, u0, p, tspan)
+    ...
+    # make sure to return an instance of ODEFunction
+end
+...
+prob = CryoGridProblem(model, tspan, p)
+```
+
+`JacobianStyle` can also be extended to create custom traits which can then be applied to compatible `LandModel`s.
+"""
+odefunction(setup::TSetup, M, u0, p, tspan; kwargs...) where {TSetup<:LandModel} = odefunction(JacobianStyle(TSetup), setup, M, u0, p, tspan; kwargs...)
+odefunction(::DefaultJac, setup::TSetup, M, u0, p, tspan; kwargs...) where {TSetup<:LandModel} = ODEFunction(setup, mass_matrix=M; kwargs...)
+function odefunction(::TridiagJac, setup::LandModel, M, u0, p, tspan; kwargs...)
+    if :jac_prototype in keys(kwargs)
+        @warn "using user specified jac_prorotype instead of tridiagonal"
+        ODEFunction(setup, mass_matrix=M; kwargs...)
+    else
+        N = length(u0)
+        J = Tridiagonal(
+                similar(u0, eltype(p), N-1) |> Vector,
+                similar(u0, eltype(p), N) |> Vector,
+                similar(u0, eltype(p), N-1) |> Vector
+        )
+        ODEFunction(setup, mass_matrix=M, jac_prototype=J, kwargs...)
+    end
+end
+
+"""
+    getstate(layername::Symbol, integrator::SciMLBase.DEIntegrator)
+
+Builds the state named tuple for `layername` given an initialized integrator.
+"""
+Land.getstate(layername::Symbol, integrator::SciMLBase.DEIntegrator) = getstate(Val{layername}(), integrator)
+Land.getstate(::Val{layername}, integrator::SciMLBase.DEIntegrator) where {layername} = Land.getstate(integrator.f.f, integrator.u, get_du(integrator), integrator.t)
+"""
+    getvar(var::Symbol, integrator::SciMLBase.DEIntegrator)
+"""
+Land.getvar(var::Symbol, integrator::SciMLBase.DEIntegrator) = Land.getvar(Val{var}(), integrator.f.f, integrator.u)
+"""
+Constructs a `CryoGridOutput` from the given `ODESolution`.
+"""
+function InputOutput.CryoGridOutput(sol::TSol) where {TSol <: SciMLBase.AbstractODESolution}
+    # Helper functions for mapping variables to appropriate DimArrays by grid/shape.
+    withdims(::Var{name,T,<:OnGrid{Cells}}, arr, grid, ts) where {name,T} = DimArray(arr*oneunit(T), (Z(round.(typeof(1.0u"m"), cells(grid), digits=5)),Ti(ts)))
+    withdims(::Var{name,T,<:OnGrid{Edges}}, arr, grid, ts) where {name,T} = DimArray(arr*oneunit(T), (Z(round.(typeof(1.0u"m"), edges(grid), digits=5)),Ti(ts)))
+    withdims(::Var{name,T}, arr, zs, ts) where {name,T} = DimArray(arr*oneunit(T), (Ti(ts),))
+    model = sol.prob.f.f # LandModel
+    ts = model.hist.vals.t # use save callback time points
+    ts_datetime = Dates.epochms2datetime.(round.(ts*1000.0))
+    u_all = reduce(hcat, sol.(ts))
+    pax = ComponentArrays.indexmap(getaxes(model.state.uproto)[1])
+    savedstates = model.hist.vals.saveval
+    allvars = variables(model)
+    progvars = tuplejoin(filter(isprognostic, allvars), filter(isalgebraic, allvars))
+    diagvars = filter(isdiagnostic, allvars)
+    fluxvars = filter(isflux, allvars)
+    outputs = Dict{Symbol,Any}()
+    for var in progvars
+        name = varname(var)
+        outputs[name] = withdims(var, u_all[pax[name],:], model.grid, ts_datetime)
+    end
+    for var in filter(isongrid, tuplejoin(diagvars, fluxvars))
+        name = varname(var)
+        states = collect(skipmissing([name ∈ keys(state) ? state[name] : missing for state in savedstates]))
+        if length(states) == length(ts)
+            arr = reduce(hcat, states)
+            outputs[name] = withdims(var, arr, model.grid, ts_datetime)
+        end
+    end
+    # loop over remaining (user defined) log variables
+    # uservars = Dict()
+    # for varname in logvarnames
+    #     var_log = log[varname]
+    #     if eltype(var_log) <: AbstractVector
+    #         vardata = reduce(hcat, var_log)
+    #         uservars[varname] = DimArray(vardata, (Z(1:size(vardata,1)),Ti(ts)))
+    #     else
+    #         uservars[varname] = DimArray(vardata, (Ti(ts),))
+    #     end
+    # end
+    # if length(logvarnames) > 0
+    #     nt = NamedTuple{tuple(keys(uservars)...)}(tuple(values(uservars)...))
+    #     layerstates = merge(layerstates, (user=nt,))
+    # end
+    CryoGridOutput(ts_datetime, sol, (;outputs...))
+end
+
+"""
+Evaluates the continuous solution at time `t`.
+"""
+(out::CryoGridOutput{<:ODESolution})(t::Real) = withaxes(out.res(t), out.res.prob.f.f)
+(out::CryoGridOutput{<:ODESolution})(t::DateTime) = out(Dates.datetime2epochms(t)/1000.0)

src/IO/InputOutput.jl

 module InputOutput
 
+using CryoGrid.Numerics
+using CryoGrid.Land
+
 using DataStructures: DefaultDict
 using Dates
+using DimensionalData
 using Downloads
 using JSON3
 using Lazy: @>>
 using Unitful
 
+export loadforcings
+
 include("ioutils.jl")
 
+export CryoGridOutput
+include("output.jl")
+
 const INPUT_DIR = "input/"
-const FORCINGS_DIR = joinpath(INPUT_DIR, "forcings")
-const PARA_DIR = joinpath(INPUT_DIR, "para")
+const DEFAULT_FORCINGS_DIR = joinpath(INPUT_DIR, "forcings")
+const DEFAULT_PARA_DIR = joinpath(INPUT_DIR, "para")
 
 struct Resource
 	name::String
@@ -42,6 +51,10 @@ struct JsonSpec{Version} <: InputSpec end
 export InputSpec, JsonSpec
 
 """
+    loadforcings(filename::String, units...; spec::Type{T}=JsonSpec{1})
+    loadforcings(resource::Resource, units...; spec::Type{T}=JsonSpec{1}, outdir=DEFAULT_FORCINGS_DIR)
+    loadforcings(::Type{JsonSpec{1}}, filename::String, units::Pair{Symbol,<:Unitful.Units}...)
+
 Loads forcing data from the given file according to the format specified by `spec::InputSpec`. Default is JsonSpec{1}.
 Returns a NamedTuple of the form `(data=(...), timestamps=Array{DateTime,1})` where `data` is a NamedTuple matching
 the structure of the JSON file. Units can (and should) be supplied as additional pair arguments, e.g:
@@ -49,20 +62,29 @@ the structure of the JSON file. Units can (and should) be supplied as additional
 `loadforcings("example.json", :Tair=>u"°C", :Ptot=>u"mm")`
 """
 loadforcings(filename::String, units...; spec::Type{T}=JsonSpec{1}) where {T <: InputSpec} = loadforcings(T, filename, units...)
-loadforcings(resource::Resource, units...; spec::Type{T}=JsonSpec{1}) where {T <: InputSpec} = loadforcings(T, fetch(resource, FORCINGS_DIR), units...)
+loadforcings(resource::Resource, units...; spec::Type{T}=JsonSpec{1}, outdir=DEFAULT_FORCINGS_DIR) where {T <: InputSpec} = loadforcings(T, fetch(resource, outdir), units...)
 function loadforcings(::Type{JsonSpec{1}}, filename::String, units::Pair{Symbol,<:Unitful.Units}...)
-        dict = open(filename, "r") do file; JSON3.read(file) end
-        # convert JSON3 dict for data field to Julia dict
-        data = Dict(dict[:data]...)
-        # get timestamps and then remove from dict
-        ts = @>> data[:t_span] map(DateNumber) map(todatetime)
-        delete!(data, :t_span)
-        unitdict = Dict(units...)
-        vals_with_units = (haskey(unitdict, name) ? vals*unitdict[name] : vals for (name, vals) in data)
-        # construct new named tuple
-        (data = NamedTuple{Tuple(keys(data))}(tuple(vals_with_units...)), timestamps = Array(ts))
+    dict = open(filename, "r") do file; JSON3.read(file) end
+    # convert JSON3 dict for data field to Julia dict
+    data = Dict(dict[:data]...)
+    # get timestamps and then remove from dict
+    ts = @>> data[:t_span] map(DateNumber) map(todatetime)
+    delete!(data, :t_span)
+    unitdict = Dict(units...)
+    vals_with_units = (haskey(unitdict, name) ? vals*unitdict[name] : vals for (name, vals) in data)
+    # construct new named tuple
+    (data = NamedTuple{Tuple(keys(data))}(tuple(vals_with_units...)), timestamps = Array(ts))
+end
+function loadforcings(::Type{JsonSpec{2}}, filename::String, units::Pair{Symbol,<:Unitful.Units}...)
+    dict = open(filename, "r") do file; JSON3.read(file) end
+    # convert JSON3 dict for data field to Julia dict
+    data = Dict(dict[:data]...)
+    # get timestamps
+    ts = map(DateTime, dict[:timestamps])
+    unitdict = Dict(units...)
+    vals_with_units = (haskey(unitdict, name) ? vals*unitdict[name] : vals for (name, vals) in data)
+    # construct new named tuple
+    (data = NamedTuple{Tuple(keys(data))}(tuple(vals_with_units...)), timestamps = Array(ts))
 end
 
-export loadforcings
-
-end
\ No newline at end of file
+end

src/IO/output.jl

+"""
+    CryoGridOutput{TRes}
+
+Helper type that stores the raw output from a CryoGrid run along with `DimArray` views of all
+logged variables. `CryoGridOutput` overrides `Base.getproperty` to allow for direct dot-syntax
+access of state variables. For example, if your model has a grid variable named `T`, `out.T` returns a `DimArray`
+with indexed time and depth axes. For OrdinaryDiffEq.jl outputs, the `ODESolution` can be accessed via `out.res`,
+or for convenience, the continuous solution at time `t` can be computed via `out(t)` which is equivalent to
+`withaxes(out.res(t))`.
+"""
+struct CryoGridOutput{TRes}
+    ts::Vector{DateTime}
+    res::TRes
+    data::NamedTuple
+    CryoGridOutput(ts::Vector{DateTime}, res::TRes, data::NamedTuple) where TRes = new{TRes}(ts, res, data)
+end
+# Overrides from Base
+function Base.show(io::IO, out::CryoGridOutput)
+    countvars(x) = 1
+    countvars(nt::NamedTuple) = sum(map(countvars, nt))
+    describe(key, val) = "$key => $(typeof(val).name.wrapper) of $(eltype(val)) with dimensions $(size(val))"
+    describe(key, val::NamedTuple) = "$key => $(map(describe, keys(val), values(val)))"
+    data = out.data
+    nvars = countvars(data)
+    println(io, "CryoGridOutput with $(length(out.ts)) time steps and $(nvars != 1 ? "$nvars variables" : "1 variable"):")
+    strs = map(describe, keys(data), values(data))
+    for r in strs
+        println(io, "   $r")
+    end
+end
+function Base.getproperty(out::CryoGridOutput, sym::Symbol)
+    if sym in (:res,:ts,:data)
+        getfield(out, sym)
+    else
+        out.data[sym]
+    end
+end

src/Setup/Setup.jl → src/Land/Land.jl

-module Setup
+module Land
 
-import ForwardDiff
-import ReverseDiff
-import Zygote
+import CryoGrid: Layer, Top, Bottom, SubSurface, Process, SubSurfaceProcess, BoundaryProcess, CompoundProcess
+import CryoGrid: variables, initialcondition!, prognosticstep!, diagnosticstep!, interact!, observe
 
-import CryoGrid.Interface: Top, Bottom
-
-using CryoGrid.Interface
 using CryoGrid.Layers
 using CryoGrid.Numerics
-using CryoGrid.Processes
+using CryoGrid.Physics
 using CryoGrid.Utils
 
-using ArraysOfArrays
 using ComponentArrays
+using ConstructionBase
 using DataStructures: OrderedDict
 using Dates
 using DimensionalData
+using IfElse
 using IntervalSets
 using Lazy: @>>, @>, groupby
 using LinearAlgebra
+using Parameters
 using Unitful
-using RecursiveArrayTools
 using Reexport
+using Setfield
+
+import Flatten
 
-@reexport using DiffEqBase: solve, init, ODEProblem, SciMLBase
-@reexport using OrdinaryDiffEq
+@reexport using ModelParameters
 @reexport using SimulationLogs
 
+export Stratigraphy, @Stratigraphy
+export StratComponent, componentname, copmonenttypes, components, boundaries
+export top, bottom, subsurface
 include("stratigraphy.jl")
-include("setup.jl")
-include("problem.jl")
-include("output.jl")
+
+export LandModelState, LayerState
+include("state.jl")
+
+export LandModel, withaxes, getstate
+include("model.jl")
+
+export ParameterVector, LinearTrend, parameters
+include("params.jl")
 
 end
\ No newline at end of file

src/Land/model.jl

+mutable struct StateHistory
+    vals::Union{Missing,<:Any}
+    StateHistory() = new(missing)
+end
+
+"""
+    AbstractLandModel{iip}
+
+Base type for 1D land models. `iip` is a value of enum `InPlaceMode` that indicates
+whether the model operates on state variables in-place (overwriting arrays) or
+out-of-place (copying arrays).
+"""
+abstract type AbstractLandModel{iip} end
+"""
+    (model::AbstractLandModel{inplace})(du,u,p,t)
+    (model::AbstractLandModel{ooplace})(u,p,t)
+
+Invokes the corresponding `step` function to compute the time derivative du/dt.
+"""
+(model::AbstractLandModel{inplace})(du,u,p,t) = step!(model,du,u,p,t)
+(model::AbstractLandModel{ooplace})(u,p,t) = step(model,u,p,t)
+
+"""
+    step!(::T,du,u,p,t) where {T<:AbstractLandModel}
+
+In-place step function for model `T`. Computes du/dt and stores the result in `du`.
+"""
+step!(::T,du,u,p,t) where {T<:AbstractLandModel} = error("no implementation of in-place step! for $T")
+"""
+    step(::T,u,p,t) where {T<:AbstractLandModel}
+
+Out-of-place step function for model `T`. Computes and returns du/dt as vector with same size as `u`.
+"""
+step(::T,u,p,t) where {T<:AbstractLandModel} = error("no implementation of out-of-place step for $T")
+
+"""
+    LandModel{TStrat,TGrid,TStates,iip,obsv} <: AbstractLandModel{iip}
+
+Defines the full specification of a CryoGrid land model; i.e. stratigraphy, grid, and state variables.
+"""
+struct LandModel{TStrat,TGrid,TStates,iip,obsv} <: AbstractLandModel{iip}
+    strat::TStrat # stratigraphy
+    grid::TGrid # grid
+    state::TStates # state variables
+    hist::StateHistory # mutable "history" type for state tracking
+    function LandModel(
+        strat::TStrat,
+        grid::TGrid,
+        state::TStates,
+        hist::StateHistory=StateHistory(),
+        iip::InPlaceMode=inplace,
+        observe::Vector{Symbol}=Symbol[]) where
+        {TStrat<:Stratigraphy,TGrid<:Grid{Edges},TStates<:VarStates}
+        new{TStrat,TGrid,TStates,iip,tuple(observe...)}(strat,grid,state,hist)
+    end
+end
+ConstructionBase.constructorof(::Type{LandModel{TStrat,TGrid,TStates,iip,obsv}}) where {TStrat,TGrid,TStates,iip,obsv} =
+    (strat, grid, state, hist) -> LandModel(strat,grid,state,hist,iip,length(obsv) > 0 ? collect(obsv) : Symbol[])
+
+Base.show(io::IO, ::MIME"text/plain", model::LandModel{TStrat,TGrid,TStates,iip,obsv}) where {TStrat,TGrid,TStates,iip,obsv} = print(io, "LandModel ($iip) with layers $(map(componentname, components(model.strat))), observables=$obsv, $TGrid, $TStrat")
+
+"""
+Constructs a `LandModel` from the given stratigraphy and grid. `arrayproto` keyword arg should be an array instance
+(of any arbitrary length, including zero, contents are ignored) that will determine the array type used for all state vectors.
+"""
+function LandModel(
+    @nospecialize(strat::Stratigraphy),
+    @nospecialize(grid::Grid{Edges,<:Numerics.Geometry,<:DistQuantity});
+    arrayproto::Type{A}=Vector,
+    iip::InPlaceMode=inplace,
+    observe::Vector{Symbol}=Symbol[],
+    chunksize=nothing,
+) where {A<:AbstractArray}
+    vars = OrderedDict()
+    components = OrderedDict()
+    for (i,comp) in enumerate(strat)
+        name = componentname(comp)
+        # build layer
+        vars[name] = _collectvars(comp)
+        params = ModelParameters.params(comp)
+        if length(params) > 0
+            # create sub-model and add layer name to all parameters
+            submodel = Model(comp)
+            submodel[:layer] = repeat([name], length(params))
+            components[name] = parent(submodel)
+        else
+            components[name] = comp
+        end
+    end
+    # rebuild stratigraphy with updated parameters
+    strat = Stratigraphy(boundaries(strat), Tuple(values(components)))
+    para = params(strat)
+    # construct state variables
+    componentnames = [componentname(node) for node in strat]
+    ntvars = NamedTuple{Tuple(componentnames)}(Tuple(values(vars)))
+    npvars = (length(filter(isprognostic, layer)) + length(filter(isalgebraic, layer)) for layer in ntvars) |> sum
+    ndvars = (length(filter(isdiagnostic, layer)) for layer in ntvars) |> sum
+    @assert (npvars + ndvars) > 0 "No variable definitions found. Did you add a method definition for CryoGrid.variables(::L,::P) where {L<:Layer,P<:Process}?"
+    @assert npvars > 0 "At least one prognostic variable must be specified."
+    chunksize = isnothing(chunksize) ? length(para) : chunksize
+    states = VarStates(ntvars, Grid(dustrip(grid), grid.geometry), chunksize, arrayproto)
+    LandModel(strat,grid,states,StateHistory(),iip,observe)
+end
+LandModel(strat::Stratigraphy, grid::Grid{Cells}; kwargs...) = LandModel(strat, edges(grid); kwargs...)
+LandModel(strat::Stratigraphy, grid::Grid{Edges,<:Numerics.Geometry,T}; kwargs...) where {T} = error("grid must have values with units of length, e.g. try using `Grid((x)u\"m\")` where `x` are your grid points.")
+# mark only stratigraphy field as flattenable
+Flatten.flattenable(::Type{<:LandModel}, ::Type{Val{:strat}}) = true
+Flatten.flattenable(::Type{<:LandModel}, ::Type{Val{name}}) where name = false
+
+"""
+Generated step function (i.e. du/dt) for any arbitrary LandModel. Specialized code is generated and compiled
+on the fly via the @generated macro to ensure type stability. The generated code updates each layer in the stratigraphy
+in sequence, i.e for each layer 1 < i < N:
+
+diagnosticstep!(layer i, ...)
+interact!(layer i-1, ...)
+prognosticstep!(layer i, ...)
+
+Note for developers: All sections of code wrapped in quote..end blocks are generated. Code outside of quote blocks
+is only executed during compilation and will not appear in the compiled version.
+"""
+@generated function step!(model::LandModel{TStrat,TGrid,TStates,inplace,obsv}, _du,_u,_p,t) where {TStrat,TGrid,TStates,obsv}
+    nodetyps = componenttypes(TStrat)
+    N = length(nodetyps)
+    expr = Expr(:block)
+    # Declare variables
+    @>> quote
+    p = updateparams!(_p, model, _du, _u, t)
+    strat = Flatten.reconstruct(model.strat, p, ModelParameters.SELECT, ModelParameters.IGNORE)
+    _du .= zero(eltype(_du))
+    du = ComponentArray(_du, getaxes(model.state.uproto))
+    u = ComponentArray(_u, getaxes(model.state.uproto))
+    state = LandModelState(model.state, boundaries(strat), u, du, t, Val{inplace}())
+    end push!(expr.args)
+    # Initialize variables for all layers
+    for i in 1:N
+        n = componentname(nodetyps[i])
+        nstate = Symbol(n,:state)
+        nlayer = Symbol(n,:layer)
+        nprocess = Symbol(n,:process)
+        @>> quote
+        $nstate = state.$n
+        $nlayer = components(strat)[$i].layer
+        $nprocess = components(strat)[$i].process
+        end push!(expr.args)
+    end
+    # Diagnostic step
+    for i in 1:N
+        n = componentname(nodetyps[i])
+        nstate = Symbol(n,:state)
+        nlayer = Symbol(n,:layer)
+        nprocess = Symbol(n,:process)
+        @>> quote
+        diagnosticstep!($nlayer,$nprocess,$nstate)
+        end push!(expr.args)
+    end
+    # Interact
+    for i in 1:N-1
+        n1,n2 = componentname(nodetyps[i]), componentname(nodetyps[i+1])
+        n1state, n2state = Symbol(n1,:state), Symbol(n2,:state)
+        n1layer, n2layer = Symbol(n1,:layer), Symbol(n2,:layer)
+        n1process, n2process = Symbol(n1,:process), Symbol(n2,:process)
+        @>> quote
+        interact!($n1layer,$n1process,$n2layer,$n2process,$n1state,$n2state)
+        end push!(expr.args)
+    end
+    # Prognostic step
+    for i in 1:N
+        n = componentname(nodetyps[i])
+        nstate = Symbol(n,:state)
+        nlayer = Symbol(n,:layer)
+        nprocess = Symbol(n,:process)
+        @>> quote
+        prognosticstep!($nlayer,$nprocess,$nstate)
+        end push!(expr.args)
+    end
+    # Observables
+    for i in 1:N
+        n = componentname(nodetyps[i])
+        nstate = Symbol(n,:state)
+        nlayer = Symbol(n,:layer)
+        nprocess = Symbol(n,:process)
+        for name in obsv
+            nameval = Val{name}()
+            @>> quote
+                observe($nameval,$nlayer,$nprocess,$nstate)
+            end push!(expr.args)
+        end
+    end
+    # make sure compiled method returns no value
+    @>> :(return nothing) push!(expr.args)
+    # emit generated expression block
+    return expr
+end
+"""
+    initialcondition!(model::LandModel, tspan::NTuple{2,Float64}, p::AbstractVector, initializers::VarInit...)
+    initialcondition!(model::LandModel, tspan::NTuple{2,DateTime}, p::AbstractVector, initializers::VarInit...)
+
+Calls `initialcondition!` on all layers/processes and returns the fully constructed u0 and du0 states.
+"""
+initialcondition!(model::LandModel, tspan::NTuple{2,DateTime}, p::AbstractVector, args...) = initialcondition!(model, convert_tspan(tspan), p, args...)
+@generated function initialcondition!(model::LandModel{TStrat,TGrid,TStates,iip,obsv}, tspan::NTuple{2,Float64}, p::AbstractVector, initializers::Numerics.VarInit...) where {TStrat,TGrid,TStates,iip,obsv}
+    nodetyps = componenttypes(TStrat)
+    N = length(nodetyps)
+    expr = Expr(:block)
+    # Declare variables
+    @>> quote
+    du = zero(similar(model.state.uproto, eltype(p)))
+    u = zero(similar(model.state.uproto, eltype(p)))
+    strat = Flatten.reconstruct(model.strat, p, ModelParameters.SELECT, ModelParameters.IGNORE)
+    state = LandModelState(model.state, boundaries(strat), u, du, tspan[1], Val{iip}())
+    end push!(expr.args)
+    # Call initializers
+    for i in 1:N
+        for j in 1:length(initializers)
+            @>> quote
+            let layerstate = state[$i],
+                init = initializers[$j];
+                if haskey(layerstate.states, varname(init))
+                    initvar!(layerstate, strat, init)
+                end
+            end
+            end push!(expr.args)
+        end
+    end
+    # Iterate over layers
+    for i in 1:N-1
+        n1,n2 = componentname(nodetyps[i]), componentname(nodetyps[i+1])
+        # create variable names
+        n1,n2 = componentname(nodetyps[i]), componentname(nodetyps[i+1])
+        n1state, n2state = Symbol(n1,:state), Symbol(n2,:state)
+        n1layer, n2layer = Symbol(n1,:layer), Symbol(n2,:layer)
+        n1process, n2process = Symbol(n1,:process), Symbol(n2,:process)
+        # generated code for layer updates
+        @>> quote
+        $n1state = state.$n1
+        $n2state = state.$n2
+        $n1layer = components(strat)[$i].layer
+        $n1process = components(strat)[$i].process
+        $n2layer = components(strat)[$(i+1)].layer
+        $n2process = components(strat)[$(i+1)].process
+        end push!(expr.args)
+        if i == 1
+            # only invoke initialcondition! for layer i on first iteration to avoid duplicated calls
+            @>> quote
+            initialcondition!($n1layer,$n1process,$n1state)
+            end push!(expr.args)
+        end
+        # invoke initialcondition! for each layer, then for both (similar to interact!)
+        @>> quote
+        initialcondition!($n2layer,$n2process,$n2state)
+        initialcondition!($n1layer,$n1process,$n2layer,$n2process,$n1state,$n2state)
+        end push!(expr.args)
+    end
+    @>> quote
+    return u
+    end push!(expr.args)
+    return expr
+end
+"""
+    initvar!(state::LayerState, ::Stratigraphy, init::VarInit{varname}) where {varname}
+    initvar!(state::LayerState, ::Stratigraphy, init::InterpInit{varname})
+
+Calls the initializer for state variable `varname`.
+"""
+initvar!(state::LayerState, ::Stratigraphy, init::Numerics.VarInit{varname}) where {varname} = init!(state[varname], init)
+initvar!(state::LayerState, ::Stratigraphy, init::Numerics.InterpInit{varname}) where {varname} = init!(state[varname], init, state.grids[varname])
+"""
+    getvar(name::Symbol, model::LandModel, u)
+    getvar(::Val{name}, model::LandModel, u)
+
+Retrieves the (diagnostic or prognostic) grid variable from `model` given prognostic state `u`.
+If `name` is not a variable in the model, or if it is not a grid variable, `nothing` is returned.
+"""
+Numerics.getvar(name::Symbol, model::LandModel, u) = getvar(Val{name}(), model, u)
+Numerics.getvar(::Val{name}, model::LandModel, u) where name = getvar(Val{name}(), model.state, withaxes(u, model))
+"""
+    getstate(layername::Symbol, model::LandModel, u, du, t)
+    getstate(::Val{layername}, model::LandModel{TStrat,TGrid,<:VarStates{layernames},iip}, _u, _du, t)
+
+Constructs a `LayerState` representing the full state of `layername` given `model`, state vectors `u` and `du`, and the
+time step `t`.
+"""
+getstate(layername::Symbol, model::LandModel, u, du, t) = getstate(Val{layername}(), model, u, du, t)
+function getstate(::Val{layername}, model::LandModel{TStrat,TGrid,<:VarStates{layernames},iip}, _u, _du, t) where {layername,TStrat,TGrid,iip,layernames}
+    du = ComponentArray(_du, getaxes(model.state.uproto))
+    u = ComponentArray(_u, getaxes(model.state.uproto))
+    i = 1
+    for j in 1:length(model.strat)
+        if layernames[j] == layername
+            i = j
+            break
+        end
+    end
+    z = boundaryintervals(map(ustrip, stripparams(boundaries(model.strat))), ustrip(model.grid[end]))[i]
+    return LayerState(model.state, z, u, du, t, Val{layername}(), Val{iip}())
+end
+"""
+    variables(model::LandModel)
+
+Returns a tuple of all variables defined in the model.
+"""
+variables(model::LandModel) = Tuple(unique(Flatten.flatten(model.state.vars, Flatten.flattenable, Var)))
+"""
+    withaxes(u::AbstractArray, ::LandModel)
+
+Constructs a `ComponentArray` with labeled axes from the given state vector `u`. Assumes `u` to be of the same type/shape
+as `setup.uproto`.
+"""
+withaxes(u::AbstractArray, model::LandModel) = ComponentArray(u, getaxes(model.state.uproto))
+withaxes(u::ComponentArray, ::LandModel) = u
+"""
+Gets the 
+"""
+function getstate(model::LandModel{TStrat,TGrid,TStates,iip}, _u, _du, t) where {TStrat,TGrid,TStates,iip}
+    du = ComponentArray(_du, getaxes(model.state.uproto))
+    u = ComponentArray(_u, getaxes(model.state.uproto))
+    return LandModelState(model.strat, model.state, u, du, t, Val{iip}())
+end
+"""
+Collects and validates all declared variables (`Var`s) for the given strat component.
+"""
+function _collectvars(@nospecialize(comp::StratComponent))
+    layer, process = comp.layer, comp.process
+    all_vars = variables(layer, process)
+    @debug "Building layer $(componentname(comp)) with $(length(all_vars)) variables: $(all_vars)"
+    # check for (permissible) duplicates between variables, excluding parameters
+    groups = groupby(var -> varname(var), all_vars)
+    for (name,gvars) in filter(g -> length(g.second) > 1, groups)
+        # if any duplicate variable deifnitions do not match, raise an error
+        @assert reduce(==, gvars) "Found one or more conflicting definitions of $name in $gvars"
+    end
+    diag_vars = filter(isdiagnostic, all_vars)
+    prog_vars = filter(isprognostic, all_vars)
+    alg_vars = filter(isalgebraic, all_vars)
+    # check for duplicated algebraic/prognostic vars
+    prog_alg_duplicated = prog_vars ∩ alg_vars
+    @assert isempty(prog_alg_duplicated) "Variables $(prog_alg_duplicated) cannot be both prognostic and algebraic."
+    # check for re-definition of diagnostic variables as prognostic
+    prog_alg = prog_vars ∪ alg_vars
+    diag_prog = filter(v -> v ∈ prog_alg, diag_vars)
+    if !isempty(diag_prog)
+        @warn "Variables $(Tuple(map(varname,diag_prog))) declared as both prognostic/algebraic and diagnostic. In-place modifications outside of callbacks may degrade integration accuracy."
+    end
+    # check for conflicting definitions of differential vars
+    diff_varnames = map(v -> Symbol(:d, varname(v)), prog_alg)
+    @assert all((isempty(filter(v -> varname(v) == d, all_vars)) for d in diff_varnames)) "Variable names $(Tuple(diff_varnames)) are reserved for differentials."
+    # prognostic takes precedence, so we remove duplicated variables from the diagnostic variable set
+    diag_vars = filter(v -> v ∉ diag_prog, diag_vars)
+    # filter remaining duplicates
+    diag_vars = unique(diag_vars)
+    prog_vars = unique(prog_vars)
+    alg_vars = unique(alg_vars)
+    # convert back to tuples
+    diag_vars, prog_vars, alg_vars = Tuple(diag_vars), Tuple(prog_vars), Tuple(alg_vars)
+    return tuplejoin(diag_vars, prog_vars, alg_vars)
+end

src/Land/params.jl

+abstract type ParamTransform end
+
+struct ParamMapping{T,name,layer}
+    transform::T
+    ParamMapping(transform::T, name::Symbol, layer::Symbol) where {T<:ParamTransform} = new{T,name,layer}(transform)
+end
+
+struct ParameterVector{T,TV,P,M} <: DenseArray{T,1}
+    vals::TV # input/reparameterized param vector
+    params::P # parameters grouped by layer and name
+    mappings::M # mapping metadata
+    ParameterVector(vals::TV,params::P,mappings::ParamMapping...) where {T,TV<:ComponentVector{T},P<:NamedTuple} = new{T,TV,P,typeof(mappings)}(vals,params,mappings)
+end
+vals(rv::ParameterVector) = getfield(rv, :vals)
+mappings(rv::ParameterVector) = getfield(rv, :mappings)
+Base.axes(rv::ParameterVector) = axes(getfield(rv, :vals))
+Base.LinearIndices(rv::ParameterVector) = LinearIndices(getfield(rv, :vals))
+Base.IndexStyle(::Type{<:ParameterVector}) = Base.IndexLinear()
+Base.similar(rv::ParameterVector) = ParameterVector(similar(vals(rv)), similar(pout(rv)), mappings(rv))
+Base.similar(rv::ParameterVector, ::Type{T}) where T = ParameterVector(similar(vals(rv), T), similar(pout(rv), T), mappings(rv))
+Base.length(rv::ParameterVector) = length(getfield(rv, :vals))
+Base.size(rv::ParameterVector) = size(getfield(rv, :vals))
+Base.getproperty(rv::ParameterVector, sym::Symbol) = getproperty(getfield(rv, :vals), sym)
+Base.getindex(rv::ParameterVector, i) = getfield(rv, :vals)[i]
+Base.setproperty!(rv::ParameterVector, val, i) = setproperty!(getfield(rv, :vals), val, sym)
+Base.setindex!(rv::ParameterVector, val, i) = setindex!(getfield(rv, :vals), val, i)
+Base.show(io, rv::ParameterVector) = show(io, getfield(rv, :vals))
+ComponentArrays.ComponentArray(rv::ParameterVector) = getfield(rv, :vals)
+
+_paramval(p::Param) = ustrip(p.val) # extracts value from Param type and strips units
+function parameters(model::LandModel, transforms::Pair{Symbol,<:Pair{Symbol,<:ParamTransform}}...)
+    function getparam(p)
+        # currently, we assume only one variable of each name in each layer;
+        # this could be relaxed in the future but will need to be appropriately handled
+        @assert length(p) == 1 "Found more than one parameter with name $var in $layer; this is not currently supported."
+        return p[1]
+    end
+    m = Model(model)
+    nestedparams = mapflat(getparam, groupparams(m, :layer, :fieldname); maptype=NamedTuple)
+    mappedparams = nestedparams
+    mappings = ParamMapping[]
+    for (layer,(var,transform)) in transforms
+        @set! mappedparams[layer][var] = mapflat(getparam, groupparams(Model(transform), :fieldname); maptype=NamedTuple)
+        push!(mappings, ParamMapping(transform, var, layer))
+    end
+    mappedarr = ComponentArray(mapflat(_paramval, mappedparams))
+    return ParameterVector(mappedarr, mappedparams, mappings...)
+end
+@inline @generated function updateparams!(v::AbstractVector, model::LandModel, u, du, t)
+    quote
+        p = ModelParameters.update(ModelParameters.params(model), v)
+        return Utils.genmap(_paramval, p)
+    end
+end
+@inline @generated function updateparams!(rv::ParameterVector{T,TV,P,M}, model::LandModel, u, du, t) where {T,TV,P,M}
+    expr = quote
+        pvals = vals(rv)
+        pmodel = ModelParameters.update(getfield(rv, :params), pvals)
+    end
+    # apply parameter transforms
+    for i in 1:length(M.parameters)
+        push!(expr.args, :(pmodel = _updateparam(pmodel, mappings(rv)[$i], model, u, du, t)))
+    end
+    # flatten parameters and strip Param types
+    push!(expr.args, :(return Utils.genmap(_paramval, ModelParameters.params(pmodel))))
+    return expr
+end
+@inline @generated function _updateparam(pmodel, mapping::ParamMapping{T,name,layer}, model::LandModel, u, du, t) where {T,name,layer}
+    quote
+        state = getstate(Val($(QuoteNode(layer))), model, u, du, t)
+        p = pmodel.$layer.$name
+        # reconstruct transform with new parameter values
+        op = ModelParameters.update(mapping.transform, ModelParameters.params(p))
+        # apply transform and replace parameter in named tuple
+        @set! pmodel.$layer.$name = Param(transform(state, op))
+        return pmodel
+    end
+end
+
+# Transform implementations
+@with_kw struct LinearTrend{P} <: ParamTransform
+    slope::P = Param(0.0)
+    intercept::P = Param(0.0)
+    tstart::Float64 = 0.0
+    tstop::Float64 = Inf; @assert tstop > tstart
+    minval::Float64 = -Inf
+    maxval::Float64 = Inf
+end
+function transform(state, trend::LinearTrend)
+    let t = min(state.t - trend.tstart, trend.tstop),
+        β = trend.slope,
+        α = trend.intercept;
+        min(max(β*t + α, trend.minval), trend.maxval)
+    end
+end

src/Land/state.jl

+using CryoGrid.Numerics: Flux
+"""
+Enumeration for in-place vs out-of-place mode.
+"""
+@enum InPlaceMode inplace ooplace
+"""
+    LayerState{iip,TStates,TGrids,Tt,Tz,varnames}
+
+Represents the state of a single component (layer + processes) in the stratigraphy.
+"""
+struct LayerState{iip,TStates,TGrids,Tt,Tz,varnames}
+    grids::NamedTuple{varnames,TGrids}
+    states::NamedTuple{varnames,TStates}
+    t::Tt
+    z::NTuple{2,Tz}
+    LayerState(grids::NamedTuple{varnames,TG}, states::NamedTuple{varnames,TS}, t::Tt, z::NTuple{2,Tz}, ::Val{iip}=Val{inplace}()) where
+        {TG,TS,Tt,Tz,varnames,iip} = new{iip,TS,TG,Tt,Tz,varnames}(grids, states, t, z)
+end
+Base.getindex(state::LayerState, sym::Symbol) = getproperty(state, sym)
+function Base.getproperty(state::LayerState, sym::Symbol)
+    return if sym ∈ (:grids, :states, :t, :z)
+        getfield(state, sym)
+    else
+        getproperty(getfield(state, :states), sym)
+    end
+end
+@inline function LayerState(vs::VarStates, zs::NTuple{2,Tz}, u, du, t, ::Val{layername}, ::Val{iip}=Val{inplace}()) where {Tz,layername,iip}
+    z_inds = subgridinds(edges(vs.grid), zs[1]..zs[2])
+    return LayerState(
+        _makegrids(Val{layername}(), getproperty(vs.vars, layername), vs, z_inds),
+        _makestates(Val{layername}(), getproperty(vs.vars, layername), vs, z_inds, u, du, t),
+        t,
+        zs,
+        Val{iip}(),
+    )
+end
+
+"""
+    LandModelState{iip,TGrid,TStates,Tt,names}
+
+Represents the instantaneous state of a CryoGrid `LandModel`.
+"""
+struct LandModelState{iip,TGrid,TStates,Tt,names}
+    grid::TGrid
+    states::NamedTuple{names,TStates}
+    t::Tt
+    LandModelState(grid::TGrid, states::NamedTuple{names,TS}, t::Tt, ::Val{iip}=Val{inplace}()) where
+        {TGrid<:Numerics.AbstractDiscretization,TS<:Tuple{Vararg{<:LayerState}},Tt,names,iip} =
+            new{iip,TGrid,TS,Tt,names}(grid, states, t)
+end
+Base.getindex(state::LandModelState, sym::Symbol) = Base.getproperty(state, sym)
+Base.getindex(state::LandModelState, i::Int) = state.states[i]
+function Base.getproperty(state::LandModelState, sym::Symbol)
+    return if sym ∈ (:grid,:states,:t)
+        getfield(state, sym)
+    else
+        getproperty(getfield(state, :states), sym)
+    end
+end
+@inline @generated function LandModelState(vs::VarStates{names}, zs::NTuple, u=copy(vs.uproto), du=similar(vs.uproto), t=0.0, ::Val{iip}=Val{inplace}()) where {names,iip}
+    layerstates = (:(LayerState(vs, (ustrip(bounds[$i][1]), ustrip(bounds[$i][2])), u, du, t, Val{$(QuoteNode(names[i]))}(), Val{iip}())) for i in 1:length(names))
+    quote
+        bounds = boundaryintervals(zs, vs.grid[end])
+        return LandModelState(
+            vs.grid,
+            NamedTuple{tuple($(map(QuoteNode,names)...))}(tuple($(layerstates...))),
+            t,
+            Val{iip}(),
+        )
+    end
+end
+# internal method dispatches for type stable construction of state types
+@inline function _makegrid(::Var{name,T,<:OnGrid{Cells}}, vs::VarStates, z_inds) where {name,T}
+    let inds=infimum(z_inds)..supremum(z_inds)-1; # subtract one due to account for cell ofset
+        get!(vs.gridcache, inds) do
+            vs.grid[inds]
+        end
+    end
+end
+@inline _makegrid(::Var{name,T,<:OnGrid{Edges}}, vs::VarStates, z_inds) where {name,T} = get!(() -> vs.grid[z_inds], vs.gridcache, z_inds)
+@inline _makegrid(var::Var, vs::VarStates, z_inds) = 1:dimlength(vardims(var), vs.grid)
+@inline _makestate(::Val, ::Prognostic{name,T,<:OnGrid{Cells}}, vs::VarStates, z_inds, u, du, t) where {name,T} = view(view(u, Val{name}()), infimum(z_inds):supremum(z_inds)-1)
+@inline _makestate(::Val, ::Prognostic{name,T}, vs::VarStates, z_inds, u, du, t) where {name,T} = view(u, Val{name}())
+@inline _makestate(::Val, ::Algebraic{name,T,<:OnGrid{Cells}}, vs::VarStates, z_inds, u, du, t) where {name,T} = view(view(u, Val{name}()), infimum(z_inds):supremum(z_inds)-1)
+@inline _makestate(::Val, ::Algebraic{name,T}, vs::VarStates, z_inds, u, du, t) where {name,T} = view(u, Val{name}())
+@inline _makestate(::Val, ::Flux{dname,name,T,<:OnGrid{Cells}}, vs::VarStates, z_inds, u, du, t) where {dname,name,T} = view(view(du, Val{name}()), infimum(z_inds):supremum(z_inds)-1)
+@inline _makestate(::Val, ::Flux{dname,name,T}, vs::VarStates, z_inds, u, du, t) where {dname,name,T} = view(du, Val{name}())
+@inline _makestate(::Val, ::Diagnostic{name,T,<:OnGrid{Cells}}, vs::VarStates, z_inds, u, du, t) where {name,T} = view(retrieve(vs.griddiag[name], u, t), infimum(z_inds):supremum(z_inds)-1)
+@inline _makestate(::Val, ::Diagnostic{name,T,<:OnGrid{Edges}}, vs::VarStates, z_inds, u, du, t) where {name,T} = view(retrieve(vs.griddiag[name], u, t), infimum(z_inds):supremum(z_inds))
+@inline _makestate(::Val{layername}, ::Diagnostic{name,T}, vs::VarStates, z_inds, u, du, t) where {name,layername,T} = retrieve(vs.diag[layername][name], u, t)
+# these need to be @generated functions in order for the compiler to infer all of the types correctly
+@inline @generated function _makegrids(::Val{layername}, vars::NamedTuple{varnames}, vs::VarStates, z_inds::ClosedInterval) where {layername,varnames}
+    quote
+        NamedTuple{tuple($(map(QuoteNode, varnames)...))}(tuple($(map(n -> :(_makegrid(vs.vars.$layername.$n, vs, z_inds)), varnames)...)))
+    end
+end
+@inline @generated function _makestates(::Val{layername}, vars::NamedTuple{varnames}, vs::VarStates, z_inds::ClosedInterval, u, du, t) where {layername,varnames}
+    quote
+        NamedTuple{tuple($(map(QuoteNode, varnames)...))}(tuple($(map(n -> :(_makestate(Val{$(QuoteNode(layername))}(), vs.vars.$layername.$n, vs, z_inds, u, du, t)), varnames)...)))
+    end
+end

src/Land/stratigraphy.jl

+"""
+    StratComponent{TLayer,TProcess,name}
+
+Represents a single component (layer + processes) in the stratigraphy.
+"""
+struct StratComponent{TLayer,TProcess,name}
+    layer::TLayer
+    process::TProcess
+    StratComponent(name::Symbol, layer::TLayer, process::TProcess) where {TLayer<:Layer,TProcess<:CompoundProcess} =
+        new{TLayer,TProcess,name}(layer,process)
+end
+ConstructionBase.constructorof(::Type{StratComponent{TLayer,TProcess,name}}) where {TLayer,TProcess,name} = (layer,process) -> StratComponent(name, layer, process)
+"""
+Get the name of the given stratigraphy node.
+"""
+componentname(::StratComponent{L,P,name}) where {L,P,name} = name
+componentname(::Type{<:StratComponent{L,P,name}}) where {L,P,name} = name
+
+Base.show(io::IO, node::StratComponent{L,P,name}) where {L,P,name} = print(io, "$name($L,$P)")
+
+# Constructors for stratigraphy nodes
+top(bcs::BoundaryProcess...) = StratComponent(:top, Top(), CompoundProcess(bcs...))
+bottom(bcs::BoundaryProcess...) = StratComponent(:bottom, Bottom(), CompoundProcess(bcs...))
+subsurface(name::Symbol, layer::SubSurface, processes::SubSurfaceProcess...) = StratComponent(name, layer, CompoundProcess(processes...))
+
+"""
+    Stratigraphy{N,TComponents,Q}
+
+Defines a 1-dimensional stratigraphy by connecting a top and bottom layer to 1 or more subsurface layers.
+"""
+struct Stratigraphy{N,TComponents,Q}
+    boundaries::NTuple{N,Q}
+    components::TComponents
+    Stratigraphy(boundaries::NTuple{N,Q}, components::NTuple{N,StratComponent}) where {N,Q} = new{N,typeof(components),Q}(boundaries, components)
+    Stratigraphy(top::Pair{Q,<:StratComponent{Top}}, sub::Pair{Q,<:StratComponent{<:SubSurface}},
+        bot::Pair{Q,<:StratComponent{Bottom}}) where {Q<:DistQuantity} = Stratigraphy(top,(sub,),bot)
+    function Stratigraphy(
+        @nospecialize(top::Pair{Q,<:StratComponent{Top}}),
+        @nospecialize(sub::Tuple{Vararg{Pair{Q,<:StratComponent{<:SubSurface}}}}),
+        @nospecialize(bot::Pair{Q,<:StratComponent{Bottom}})
+    ) where {Q<:DistQuantity}
+        @assert length(sub) > 0 "At least one subsurface layer must be specified"
+        names = map(componentname, map(last, sub))
+        @assert length(unique(names)) == length(names) "All layer names in Stratigraphy must be unique"
+        boundaries = Tuple(map(pair -> Param(ustrip(first(pair)), units=unit(Q), layer=:strat), (top, sub..., bot)))
+        @assert issorted(boundaries, by=p -> p.val) "Stratigraphy boundary locations must be in strictly increasing order."
+        components = Tuple(map(last, (top, sub..., bot)))
+        new{length(components),typeof(components),eltype(boundaries)}(boundaries,components)
+    end
+end
+components(strat::Stratigraphy) = getfield(strat, :components)
+boundaries(strat::Stratigraphy) = getfield(strat, :boundaries)
+boundaryintervals(strat::Stratigraphy, z_bottom) = boundaryintervals(boundaries(strat), z_bottom)
+boundaryintervals(bounds::NTuple, z_bottom) = tuplejoin(map(tuple, bounds[1:end-1], bounds[2:end]), ((bounds[end], z_bottom),))
+componenttypes(::Type{<:Stratigraphy{N,TComponents}}) where {N,TComponents} = Tuple(TComponents.parameters)
+Base.getproperty(strat::Stratigraphy, sym::Symbol) = strat[Val{sym}()]
+# Array and iteration overrides
+Base.size(strat::Stratigraphy) = size(components(strat))
+Base.length(strat::Stratigraphy) = length(components(strat))
+Base.getindex(strat::Stratigraphy, i::Int) = components(strat)[i]
+Base.getindex(strat::Stratigraphy, sym::Symbol) = strat[Val{sym}]
+Base.getindex(strat::Stratigraphy, ::Val{sym}) where {sym} = strat[Val{sym}]
+@generated Base.getindex(strat::Stratigraphy{N,TC}, ::Type{Val{sym}}) where {N,TC,sym} = :(components(strat)[$(findfirst(T -> componentname(T) == sym, TC.parameters))])
+Base.iterate(strat::Stratigraphy) = (components(strat)[1],components(strat)[2:end])
+Base.iterate(strat::Stratigraphy, itrstate::Tuple) = (itrstate[1],itrstate[2:end])
+Base.iterate(strat::Stratigraphy, itrstate::Tuple{}) = nothing
+Base.show(io::IO, strat::Stratigraphy) = print(io, "Stratigraphy($(prod(("$b => $n, " for (n,b) in zip(components(strat),boundaries(strat)))))")
+
+"""
+Convenience macro for defining stratigraphies with multiple subsurface layers.
+"""
+macro Stratigraphy(args...)
+    @assert length(args) >= 3 "At least three stratigraphy nodes (top, subsurface, bottom) must be provided!"
+    if length(args) == 3
+        :(Stratigraphy($(esc(args[1])), $(esc(args[2])), $(esc(args[3]))))
+    elseif length(args) > 3
+        :(Stratigraphy($(esc(args[1])), tuple($(esc.(args[2:end-1])...)), $(esc(args[end]))))
+    end
+end

src/Layers/Layers.jl

 module Layers
 
-import CryoGrid.Interface: SubSurface, initialcondition!, variables
+import CryoGrid: SubSurface
+import CryoGrid: initialcondition!, variables
 
-using CryoGrid.Interface
 using CryoGrid.Numerics
 using CryoGrid.Utils
 
 using DimensionalData
 using IntervalSets
+using ModelParameters
 using Parameters
 using Unitful
 
+export Soil, SoilParameterization, SoilCharacteristicFractions, SoilProfile, SoilType, Sand, Silt, Clay, soilparameters
 include("soil.jl")
 
 end
\ No newline at end of file

src/Layers/soil.jl

 """
-Represents the general type of the soil. Sand, Silt, and Clay are provided by default.
+Represents the texture classification of the soil. Sand, Silt, and Clay are provided by default.
 """
-abstract type SoilType end
-struct Sand <: SoilType end
-struct Silt <: SoilType end
-struct Clay <: SoilType end
+abstract type SoilTexture end
+struct Sand <: SoilTexture end
+struct Silt <: SoilTexture end
+struct Clay <: SoilTexture end
 """
-Base type for soil parameterizations.
+    SoilParameterization
+
+Abstract base type for parameterizations of soil properties.
 """
 abstract type SoilParameterization end
-struct ByComposition <: SoilParameterization end
-struct ByXicePorSat <: SoilParameterization end
 """
+    SoilCharacteristicFractions{P1,P2,P3,P4} <: SoilParameterization
+
 Represents the composition of the soil in terms of fractions: excess ice, natural porosity, saturation, and organic/(mineral + organic).
 """
-@with_kw struct SoilProperties <: Params @deftype Float64
-    χ; @assert ϕ >= 0.0 && ϕ <= 1.0 # natural porosity
-    ϕ; @assert θ >= 0.0 && θ <= 1.0 # saturation
-    θ; @assert ω >= 0.0 && ω <= 1.0 # organic fraction of solid; mineral fraction is 1-ω
-    ω; @assert χ >= 0.0 && χ <= 1.0 # excess ice fraction
+struct SoilCharacteristicFractions{P1,P2,P3,P4} <: SoilParameterization
+    xic::P1 # excess ice fraction
+    por::P2 # natural porosity
+    sat::P3 # saturation
+    org::P4 # organic fraction of solid; mineral fraction is 1-org
+    SoilCharacteristicFractions(xic::P1, por::P2, sat::P3, org::P4) where {P1,P2,P3,P4} = new{P1,P2,P3,P4}(xic,por,sat,org)
+end
+function soilparameters(;xic=0.0, por=0.5, sat=1.0, org=0.5)
+    params = Tuple(Param(p, bounds=(0.0,1.0)) for p in [xic,por,sat,org])
+    SoilCharacteristicFractions(params...)
 end
+SoilProfile(pairs::Pair{<:DistQuantity,<:SoilParameterization}...) = Profile(pairs...)
+# Helper functions for obtaining soil compositions from characteristic fractions.
+soilcomp(::Val{:θp}, χ, ϕ, θ, ω) = (1-χ)*ϕ
+soilcomp(::Val{:θw}, χ, ϕ, θ, ω) = χ + (1-χ)*ϕ*θ
+soilcomp(::Val{:θm}, χ, ϕ, θ, ω) = (1-χ)*(1-ϕ)*(1-ω)
+soilcomp(::Val{:θo}, χ, ϕ, θ, ω) = (1-χ)*(1-ϕ)*ω
+θx(fracs::SoilCharacteristicFractions) = fracs.xic
+θp(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θp}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
+θw(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θw}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
+θm(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θm}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
+θo(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θo}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
 """
 Thermal conductivity constants.
 """
-@with_kw struct SoilTCParams <: Params @deftype Float"W/(m*K)"
+@with_kw struct SoilTCParams @deftype Float"W/(m*K)"
     kw = 0.57xu"W/(m*K)" #water [Hillel(1982)]
     ko = 0.25xu"W/(m*K)" #organic [Hillel(1982)]
     km = 3.8xu"W/(m*K)" #mineral [Hillel(1982)]
@@ -33,101 +51,20 @@ end
 """
 Heat capacity constants.
 """
-@with_kw struct SoilHCParams <: Params @deftype Float"J/(K*m^3)"
-    cw = 4.2*10^6xu"J/(K*m^3)" #[J/m^3K] heat capacity water
-    co = 2.5*10^6xu"J/(K*m^3)" #[J/m^3K]  heat capacity organic
-    cm = 2*10^6xu"J/(K*m^3)" #[J/m^3K]  heat capacity mineral
-    ca = 0.00125*10^6xu"J/(K*m^3)" #[J/m^3K]  heat capacity pore space
-    ci = 1.9*10^6xu"J/(K*m^3)" #[J/m^3K]  heat capacity ice
+@with_kw struct SoilHCParams @deftype Float"J/(K*m^3)"
+    cw = 4.2e6xu"J/(K*m^3)" #[J/m^3K] heat capacity water
+    co = 2.5e6xu"J/(K*m^3)" #[J/m^3K]  heat capacity organic
+    cm = 2e6xu"J/(K*m^3)" #[J/m^3K]  heat capacity mineral
+    ca = 0.00125e6xu"J/(K*m^3)" #[J/m^3K]  heat capacity pore space
+    ci = 1.9e6xu"J/(K*m^3)" #[J/m^3K]  heat capacity ice
 end
 """
-Parameter type for Soil layers, includes thermal conductivity and heat capacity
-constants as well as type/composition.
+Basic Soil layer.
 """
-@with_kw struct SoilParams{TType<:SoilType,S} <: Params
+@with_kw struct Soil{T,P<:SoilParameterization,S} <: SubSurface
+    texture::T = Sand()
+    para::P = soilparameters()
     tc::SoilTCParams = SoilTCParams()
     hc::SoilHCParams = SoilHCParams()
-    type::TType = Sand()
     sp::S = nothing # user-defined specialization
 end
-"""
-Basic Soil layer.
-"""
-struct Soil{TType,TPara,S} <: SubSurface
-    profile::DimArray
-    params::SoilParams{TType,S}
-    function Soil(profile::DimArray, ::TPara=Nonparametric(); kwargs...) where {P<:SoilParameterization,TPara<:Union{Nonparametric,Parametric{P}}}
-        params = SoilParams(;kwargs...)
-        new{typeof(params.type),TPara,typeof(params.sp)}(profile, params)
-    end
-end
-
-"""
-Alias/constructor for soil profile.
-"""
-function SoilProfile(vals::Pair{<:DistQuantity,SoilProperties}...)
-    points = [d => tuple(props...) for (d,props) in vals]
-    Profile(points...;names=fieldnames(SoilProperties))
-end
-
-export Soil, SoilProperties, SoilProfile, SoilParams, SoilType, Sand, Silt, Clay, ByComposition, ByXicePorSat
-
-# Helper functions for obtaining soil component fractions from soil properties.
-soilcomp(::Val{:θx}, χ, ϕ, θ, ω) = χ
-soilcomp(::Val{:θp}, χ, ϕ, θ, ω) = @. (1-χ)*ϕ*θ
-soilcomp(::Val{:θm}, χ, ϕ, θ, ω) = @. (1-χ)*(1-ϕ)*(1-ω)
-soilcomp(::Val{:θo}, χ, ϕ, θ, ω) = @. (1-χ)*(1-ϕ)*ω
-
-Base.show(io::IO, soil::Soil{T,P}) where {T,P} = print(io, "Soil{$T,$P}($(soil.params))")
-
-variables(::Soil{T,Nonparametric}) where T = (
-    Diagnostic(:θw, Float64, OnGrid(Cells)),
-    Diagnostic(:θp, Float64, OnGrid(Cells)),
-    Diagnostic(:θx, Float64, OnGrid(Cells)),
-    Diagnostic(:θl, Float64, OnGrid(Cells)),
-    Diagnostic(:θm, Float64, OnGrid(Cells)),
-    Diagnostic(:θo, Float64, OnGrid(Cells)),
-)
-
-variables(soil::Soil{T,Parametric{ByComposition}}) where T = (
-    Diagnostic(:θw, Float64, OnGrid(Cells)),
-    Diagnostic(:θp, Float64, OnGrid(Cells)),
-    Diagnostic(:θx, Float64, OnGrid(Cells)),
-    Diagnostic(:θl, Float64, OnGrid(Cells)),
-    Diagnostic(:θm, Float64, OnGrid(Cells)),
-    Diagnostic(:θo, Float64, OnGrid(Cells)),
-    Parameter(:θx, soilcomp(Val{:θx}(), soil.profile[Y(:χ)], soil.profile[Y(:ϕ)], soil.profile[Y(:θ)], soil.profile[Y(:ω)]), 0..1),
-    Parameter(:θp, soilcomp(Val{:θp}(), soil.profile[Y(:χ)], soil.profile[Y(:ϕ)], soil.profile[Y(:θ)], soil.profile[Y(:ω)]), 0..1),
-    Parameter(:θm, soilcomp(Val{:θm}(), soil.profile[Y(:χ)], soil.profile[Y(:ϕ)], soil.profile[Y(:θ)], soil.profile[Y(:ω)]), 0..1),
-    Parameter(:θo, soilcomp(Val{:θo}(), soil.profile[Y(:χ)], soil.profile[Y(:ϕ)], soil.profile[Y(:θ)], soil.profile[Y(:ω)]), 0..1),
-)
-
-variables(soil::Soil{T,Parametric{ByXicePorSat}}) where T = (
-    Diagnostic(:θw, Float64, OnGrid(Cells)),
-    Diagnostic(:θp, Float64, OnGrid(Cells)),
-    Diagnostic(:θx, Float64, OnGrid(Cells)),
-    Diagnostic(:θl, Float64, OnGrid(Cells)),
-    Diagnostic(:θm, Float64, OnGrid(Cells)),
-    Diagnostic(:θo, Float64, OnGrid(Cells)),
-    Parameter(:χ, soil.profile[var=:χ], 0..1),
-    Parameter(:ϕ, soil.profile[var=:ϕ], 0..1),
-    Parameter(:θ, soil.profile[var=:θ], 0..1),
-    Parameter(:ω, soil.profile[var=:ω], 0..1),
-)
-
-function initialcondition!(soil::Soil{T,P}, state) where {T,P}
-    # Helper functions for initializing soil composition state based on parameterization mode.
-    fromparams(::Val{var}, soil::Soil{T,Parametric{ByComposition}}, state) where {var,T} = state.params[var]
-    fromparams(::Val{var}, soil::Soil{T,Parametric{ByXicePorSat}}, state) where {var,T} = soilcomp(Val{var}(), state.params.χ, state.params.ϕ, state.params.θ, state.params.ω)
-    fromparams(::Val{var}, soil::Soil{T,Nonparametric}, state) where {var,T} = soilcomp(Val{var}(), soil.profile[var=:χ], soil.profile[var=:ϕ], soil.profile[var=:θ], soil.profile[var=:ω])
-    depths = length(size(soil.profile)) > 1 ? dims(soil.profile, :depth).val : [refdims(soil.profile)[1].val]
-    for var in [:θx,:θp,:θm,:θo]
-        arr = DimArray(similar(state[var], Union{Missing,eltype(state[var])}), (depth=state.grids[var]u"m",))
-        arr .= missing
-        arr_sub = @view arr[depth=Near(depths)]
-        arr_sub .= fromparams(Val{var}(), soil, state)
-        Utils.ffill!(arr)
-        state[var] .= skipmissing(arr)
-    end
-    @. state.θw = state.θx + state.θp
-end

src/Common/Numerics/Numerics.jl → src/Numerics/Numerics.jl

@@ -2,11 +2,16 @@ module Numerics
 
 import Base.==
 import ExprTools
+import ForwardDiff
+import PreallocationTools as Prealloc
 
 using CryoGrid.Utils
 
 using Base: @inbounds, @propagate_inbounds
-using DimensionalData: DimArray, dims
+using ConstructionBase
+using ComponentArrays
+using DimensionalData: AbstractDimArray, DimArray, Dim, At, dims, Z
+using Flatten
 using IfElse
 using Interpolations: Interpolations, Gridded, Linear, Flat, Line, interpolate, extrapolate
 using IntervalSets
@@ -14,13 +19,42 @@ using LinearAlgebra
 using LoopVectorization
 using RuntimeGeneratedFunctions
 using Unitful
+using StructTypes
 using Symbolics
 using SymbolicUtils
 
 RuntimeGeneratedFunctions.init(@__MODULE__)
 
+export GridSpec, Edges, Cells
+
+abstract type AbstractDiscretization{Q,N} <: DenseArray{Q,N} end
+
+abstract type GridSpec end
+struct Edges <: GridSpec end
+struct Cells <: GridSpec end
+
+abstract type Geometry end
+struct UnitVolume <: Geometry end
+
+export ∇
 include("math.jl")
+
+export Grid, cells, edges, indexmap, subgridinds, Δ, volume, area
 include("grid.jl")
+
+export Profile, profile2array, array2profile
+include("profile.jl")
+
+export Var, Prognostic, Algebraic, Diagnostic, VarDim, OnGrid, Shape, Scalar
+export varname, vartype, vardims, isprognostic, isalgebraic, isflux, isdiagnostic, isongrid, dimlength
 include("variables.jl")
 
+export VarStates, DiffCache, retrieve, getvar, getvars
+include("varstates.jl")
+
+include("discretize.jl")
+
+export initializer, init!
+include("init.jl")
+
 end

src/Numerics/discretize.jl

+"""
+    discretize([::Type{A}], ::T,  ::Var) where {T,N,D<:AbstractDiscretization{T,N},A<:AbstractArray{T,N}}
+
+Produces a discretization of the given variable based on `T` and array type `A`.
+"""
+discretize(::Type{A}, ::D, ::Var) where {Q,T,N,D<:AbstractDiscretization{Q,N},A<:AbstractArray{T,N}} = error("missing discretize implementation for $D")
+discretize(d::AbstractDiscretization{Q,N}, var::Var) where {Q,N} = discretize(Array{vartype(var),N}, d, var)
+# grid discretizations
+discretize(::Type{A}, grid::Grid, var::Var) where {A<:AbstractVector} = zero(similar(A{vartype(var)}, dimlength(var.dim, grid)))
+function discretize(::Type{A}, grid::Grid, pvars::Union{<:Prognostic,<:Algebraic}...) where {A<:AbstractVector}
+    # separate into grid and non-grid vars
+    gridvars = unique(filter(v -> isa(vardims(v), OnGrid), pvars))
+    pointvars = filter(v -> !isa(vardims(v), OnGrid), pvars)
+    # get lengths
+    gridvar_ns = map(v -> dimlength(vardims(v), grid), gridvars)
+    pointvar_ns = map(v -> dimlength(vardims(v), grid), pointvars)
+    Ng = length(gridvar_ns) > 0 ? sum(gridvar_ns) : 0
+    Np = length(pointvar_ns) > 0 ? sum(pointvar_ns) : 0
+    # build axis indices;
+    # non-grid prognostic variables get collected at the top of the vector, in the order provided
+    pointvar_ax = (;(varname(p) => i:(i+n) for (p,n,i) in zip(pointvars, pointvar_ns, cumsum(vcat([1],collect(pointvar_ns[1:end-1])))))...)
+    # grid variables get interlaced throughout the rest of the vector; i.e. for variable i, its grid points are:
+    # i:k:kn where k is the number of grid variables and n is the length of the grid.
+    gridvar_ax = (;(varname(p) => st:length(gridvars):Ng for (p,st) in zip(gridvars, (Np+1):(1+Np+length(gridvars))))...)
+    # allocate component array; assumes all variables have (and should!) have the same type
+    u = zero(similar(A{vartype(first(pvars))}, Ng+Np))
+    ComponentVector(u, (Axis(merge(pointvar_ax, gridvar_ax)),))
+end

src/Common/Numerics/grid.jl → src/Numerics/grid.jl

-abstract type GridSpec end
-struct Edges <: GridSpec end
-struct Cells <: GridSpec end
-
-export GridSpec, Edges, Cells
-
-abstract type Geometry end
-struct UnitVolume <: Geometry end
+const GridValues{Q,A} = NamedTuple{(:edges,:cells),NTuple{2,SubArray{Q,1,A,Tuple{StepRange{Int64,Int64}},true}}} where {Q,A<:AbstractArray}
 
 """
-    struct Grid{S,G,A,Q} <: AbstractVector{Q}
+    struct Grid{S,G,Q,A} <: AbstractDiscretization{Q,1}
 
 Represents the 1D spatial discretization on which time integration is performed. `S` is a `GridSpec`,
 either `Edges` or `Cells` (always edges upon initial construction). The grid representation can be
@@ -16,11 +9,12 @@ converted (allocation free) between grid edges and cells via the `cells` and `ed
 represents the geometry/volume on which the vertical 1D discretization is applied. `A` is the underlying
 array type, and `Q` is the numerical type (e.g. `Float64` or a `Unitful.Quantity`).
 """
-struct Grid{S,G,A,Q} <: AbstractVector{Q}
+struct Grid{S,G,Q,A} <: AbstractDiscretization{Q,1}
     geometry::G
-    values::NamedTuple{(:edges,:cells),NTuple{2,SubArray{Q,1,A,Tuple{StepRange{Int64,Int64}},true}}}
-    deltas::NamedTuple{(:edges,:cells),NTuple{2,SubArray{Q,1,A,Tuple{StepRange{Int64,Int64}},true}}}
-    function Grid(vals::TArray, geometry::G=UnitVolume()) where {G<:Geometry,Q<:Number,TArray<:AbstractArray{Q,1}}
+    values::GridValues{Q,A}
+    deltas::GridValues{Q,A}
+    Grid(::Type{S}, values::GridValues{Q,A}, deltas::GridValues{Q,A}, geom::G) where {S<:GridSpec,Q,A,G<:Geometry} = new{S,G,Q,A}(geom,values,deltas)
+    function Grid(vals::AbstractVector{Q}, geometry::G=UnitVolume()) where {G<:Geometry,Q<:Number}
         @assert issorted(vals) "grid values should be in ascending order"
         nedges = length(vals)
         ncells = nedges - 1 
@@ -31,43 +25,41 @@ struct Grid{S,G,A,Q} <: AbstractVector{Q}
         # midpoints should always be the same as distance between boundaries...
         Δedges = edges[2:end] .- edges[1:end-1]
         Δcells = cells[2:end] .- cells[1:end-1]
-        new{Edges,G,typeof(edges),Q}(geometry,
+        new{Edges,G,Q,typeof(edges)}(geometry,
             (edges=view(edges,1:1:nedges),cells=view(cells,1:1:ncells)),
             (edges=view(Δedges,1:1:nedges-1),cells=view(Δcells,1:1:ncells-1)),
         )
     end
-    function Grid(grid::Grid{Edges,G,A,Q}, interval::ClosedInterval{Int}) where {G,A,Q}
+    function Grid(grid::Grid{Edges,G,Q,A}, interval::ClosedInterval{Int}) where {G,Q,A}
         start, stop = interval.left, interval.right
         edges = @view grid.values.edges[start:stop]
         cells = @view grid.values.cells[start:stop-1]
         Δedges = @view grid.deltas.edges[start:stop-1]
         Δcells = @view grid.deltas.cells[start:stop-2]
-        new{Edges,G,A,Q}(grid.geometry, (edges=edges,cells=cells), (edges=Δedges,cells=Δcells))
+        new{Edges,G,Q,A}(grid.geometry, (edges=edges,cells=cells), (edges=Δedges,cells=Δcells))
     end
-    function Grid(grid::Grid{Cells,G,A,Q}, interval::ClosedInterval{Int}) where {G,A,Q}
+    function Grid(grid::Grid{Cells,G,Q,A}, interval::ClosedInterval{Int}) where {G,Q,A}
         start, stop = interval.left, interval.right
         edges = @view grid.values.edges[start:stop+1]
         cells = @view grid.values.cells[start:stop]
         Δedges = @view grid.deltas.edges[start:stop]
         Δcells = @view grid.deltas.cells[start:stop-1]
-        new{Cells,G,A,Q}(grid.geometry, (edges=edges,cells=cells), (edges=Δedges,cells=Δcells))
+        new{Cells,G,Q,A}(grid.geometry, (edges=edges,cells=cells), (edges=Δedges,cells=Δcells))
     end
-    Grid(::Type{Cells}, grid::Grid{Edges,G,A,Q}) where {G,A,Q} = new{Cells,G,A,Q}(grid.geometry,grid.values,grid.deltas)
-    Grid(::Type{Edges}, grid::Grid{Cells,G,A,Q}) where {G,A,Q} = new{Edges,G,A,Q}(grid.geometry,grid.values,grid.deltas)
+    Grid(::Type{Cells}, grid::Grid{Edges,G,Q,A}) where {G,Q,A} = new{Cells,G,Q,A}(grid.geometry,grid.values,grid.deltas)
+    Grid(::Type{Edges}, grid::Grid{Cells,G,Q,A}) where {G,Q,A} = new{Edges,G,Q,A}(grid.geometry,grid.values,grid.deltas)
 end
-
+ConstructionBase.constructorof(::Type{Grid{S,G,Q,A}}) where {S,G,Q,A} = (geom,values,deltas) -> Grid(S,values,deltas,geom)
 Base.show(io::IO, grid::Grid{S,G}) where {S,G} = print(io, "Grid{$S}($(grid[1])..$(grid[end])) of length $(length(grid)) with geometry $G")
 Base.show(io::IO, ::MIME{Symbol("text/plain")}, grid::Grid) = show(io, grid)
 
-function subgrid(grid::Grid{S,G,A,Q}, interval::Interval{L,R,Q}) where {S,G,A,Q,L,R}
-    l,r = interval.left,interval.right
+function subgridinds(grid::Grid{S,G,Q,A}, interval::Interval{L,R,Q}) where {S,G,Q,A,L,R}
+    @assert interval.left <= interval.right "Invalid interval: $interval"
     vals = values(grid)
     # Determine indices which lie in the given interval
-    l_ind = findfirst(x -> x ∈ interval, vals)
-    r_ind = findlast(x -> x ∈ interval, vals)
-    @assert !isnothing(l_ind) && !isnothing(r_ind) "No grid points in the given interval $interval"
-    # Map back to full grid indices
-    Grid(grid,l_ind..r_ind)
+    l_ind = searchsortedfirst(vals, interval.left)
+    r_ind = searchsortedlast(vals, interval.right)
+    return (L == :closed ? l_ind : l_ind + 1)..(R == :closed ? r_ind : r_ind - 1)
 end
 @inline Δ(grid::Grid{Edges}) = grid.deltas.edges
 @inline Δ(grid::Grid{Cells}) = grid.deltas.cells
@@ -82,47 +74,10 @@ end
 @inline Base.size(grid::Grid) = size(values(grid))
 @inline Base.length(grid::Grid) = length(values(grid))
 @propagate_inbounds Base.getindex(grid::Grid, i::Int) = values(grid)[i]
-@propagate_inbounds Base.getindex(grid::Grid{S,G,A,Q}, interval::Interval{L,R,Q}) where {S,G,A,Q,L,R} = subgrid(grid,interval)
+@propagate_inbounds Base.getindex(grid::Grid{S,G,Q,A}, interval::Interval{L,R,Q}) where {S,G,Q,A,L,R} = Grid(grid, subgridinds(grid,interval))
+@propagate_inbounds Base.getindex(grid::Grid, interval::Interval{L,R,Int}) where {L,R} = Grid(grid, interval)
 Base.setindex!(grid::Grid, args...) = error("setindex! is not allowed for Grid types")
 
-regrid(x::AbstractVector, xgrid::Grid, newgrid::Grid, interp=Linear(), bc=Line()) =
-    regrid!(similar(x,length(newgrid)), x, xgrid, newgrid, interp, bc)
-function regrid!(out::AbstractVector, x::AbstractVector, xgrid::Grid, newgrid::Grid, interp=Linear(), bc=Line())
-    let f = extrapolate(interpolate((xgrid,), x, Gridded(interp)), bc)
-        out .= f.(newgrid)
-    end
-end
-
-export Grid, cells, edges, indexmap, subgrid, Δ, regrid, regrid!
-
-"""
-    Profile(pairs...;names)
-
-Constructs a Profile from the given pairs Q => (x1,...,xn) where x1...xn are the values defined at Q.
-Column names for the resulting DimArray can be set via the names parameter which accepts an NTuple of symbols,
-where N must match the number of parameters given (i.e. n).
-"""
-function Profile(pairs::Pair{Q,NTuple{N,T}}...;names::Union{Nothing,NTuple{N,Symbol}}=nothing) where {T,N,Q<:DistQuantity}
-    depths, vals = zip(pairs...)
-    params = hcat(collect.(vals)...)'
-    names = isnothing(names) ? [Symbol(:x,:($i)) for i in 1:N] : collect(names)
-    DimArray(params, (depth=collect(depths), var=names))
-end
-
-"""
-    interpolateprofile(profile::Profile, state; interp=Linear())
-
-Interpolates the given profile to the corresponding variable grids. Assumes state to be indexable via the corresponding
-variable symbol and that the parameter names in state and profile match.
-"""
-function interpolateprofile!(profile::DimArray, state; interp=Linear())
-    let (depths,names) = dims(profile),
-        z = ustrip.(depths);
-        for p in names
-            f = extrapolate(interpolate((z,), profile[:,p], Gridded(interp)), Flat())
-            state[p] .= f.(state.grids[p]) .|> dustrip   # assume length(grid) == length(state.p)
-        end
-    end
-end
-
-export Profile, interpolateprofile!
+# unit volume
+@inline volume(grid::Grid{Cells,UnitVolume,Q}) where Q = Δ(edges(grid)).*oneunit(Q)^2
+@inline area(::Grid{Edges,UnitVolume,Q}) where Q = oneunit(Q)^2