"""
    spinup(setup::Tile, 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::Tile, tspan::NTuple{2,DateTime}, p, tol, layername; maxdepth=100u"m", maxiter=1000, saveat=3*3600.0, solver=CGEuler(), dt=60.0, solve_args...)
    u0, du0 = initialcondition!(tile, tspan, p)
    prob = CryoGridProblem(setup, u0, tspan, p)
    @info "Running initial solve ..."
    sol = solve(prob, solver, dt=dt, saveat=saveat, solve_args...)
    out = CryoGridOutput(sol)
    H = out.vars[layername].H
    initialstate = getstate(setup, prob.u0, similar(prob.u0), prob.tspan[1])
    grid = initialstate[layername].grid
    max_ind = argmin(abs.(grid*1u"m" .- maxdepth))
    dz = Δ(grid)[1:max_ind]
    H_sub = H[1:max_ind,:]
    # transpose and dot with cell size to integrate over grid
    tspan_mean = mean(H_sub' * dz)
    tspan_mean_prev = Inf
    Δμ = abs.(tspan_mean - tspan_mean_prev)
    itercount = 1
    while Δμ > tol && itercount <= maxiter
        @info "[Iteration $itercount] tspan mean: $tspan_mean, Δμ: $Δμ"
        tspan_mean_prev = tspan_mean
        prob = remake(prob, u0=sol.u[end])
        sol = solve(prob, solver, dt=dt, saveat=saveat, solve_args...)
        out = CryoGridOutput(sol)
        H = out.vars[layername].H
        H_sub = H[1:max_ind,:]
        tspan_mean = mean(H_sub' * dz)
        Δμ = abs.(tspan_mean - tspan_mean_prev)
        itercount += 1
    end
    if Δμ > tol
        @warn "Energy state did not converge! Stopping after $maxiter iterations; tspan mean: $tspan_mean, Δμ: $Δμ"
    else
        @info "Finished after $itercount iterations; tspan mean: $tspan_mean, Δμ: $Δμ"
    end
    return out
end
