diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl
index dc4d2da8df057789a3e340f0be33ab03579a4254..2cc4cbaafb12ce4dda6f7e54e43ab359e24f4f69 100644
--- a/examples/heat_freeW_lake_lite_implicit.jl
+++ b/examples/heat_freeW_lake_lite_implicit.jl
@@ -6,7 +6,12 @@ Plots.plotly()
 
 CryoGrid.debug(true)
 
-forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+raw_forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+Tair = raw_forcings.data.Tair
+Ptot = uconvert.(u"m/s", raw_forcings.data.Ptot)
+rainfall = Ptot.*(Tair .> 0u"°C")
+snowfall = Ptot.*(Tair .<= 0u"°C")
+forcings = rebuild(raw_forcings; Tair, rainfall, snowfall);
 soilprofile = SoilProfile(
     0.0u"m" => SimpleSoil(por=0.80,sat=0.9,org=0.75),
     0.1u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.25),
diff --git a/examples/heat_freeW_lite_implicit.jl b/examples/heat_freeW_lite_implicit.jl
index aa88e5574dd0859035b4ec239ea82f74d14e46ad..f818f2f8c8382bbe5ec0824c6f38de2b2c7922de 100644
--- a/examples/heat_freeW_lite_implicit.jl
+++ b/examples/heat_freeW_lite_implicit.jl
@@ -7,11 +7,17 @@
 using CryoGrid
 using CryoGrid.LiteImplicit
 
-# Load forcings and build stratigraphy like before.
-forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
-forcings = Base.rename(forcings, :Ptot => :precip)
-precip_values = forcings.precip.interpolant.coefs
-precip_values .*= 2
+# First we prepare the forcing data. This data includes only air temperature and
+# *total* precipitation, so we need to separate this into rainfall and snowfall
+# using air temperature.
+raw_forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+Tair = raw_forcings.data.Tair
+Ptot = uconvert.(u"m/s", raw_forcings.data.Ptot)
+rainfall = Ptot.*(Tair .> 0u"°C")
+snowfall = Ptot.*(Tair .<= 0u"°C")
+forcings = rebuild(raw_forcings; Tair, rainfall, snowfall);
+
+# Next we define the boundary conditions and stratigraphy.
 z_top = -2.0u"m"
 z_bot = 1000.0u"m"
 upperbc = WaterHeatBC(
@@ -29,7 +35,7 @@ soilprofile = SoilProfile(
 )
 heat = HeatBalance(heatop)
 water = WaterBalance()
-soil_layers = map(para -> Ground(para; heat, water), soilprofile)
+soil_layers = map(para -> Ground(para; heat, water), soilprofile);
 strat = @Stratigraphy(
     z_top => Top(upperbc),
     z_top => Snowpack(Snow.LiteGridded(); heat),
@@ -40,11 +46,11 @@ modelgrid = Grid(vcat(z_top:0.02u"m":-0.02u"m", CryoGrid.DefaultGrid_2cm))
 tile = Tile(strat, modelgrid, forcings, ssinit);
 
 # Since the solver can take daily timesteps, we can easily specify longer simulation time spans at minimal cost.
-# Here we specify a time span of 30 years.
-tspan = (DateTime(1990,10,1), DateTime(2020,10,1))
+# Here we specify a time span of 20 years.
+tspan = (DateTime(2000,10,1), DateTime(2020,10,1))
 u0, du0 = initialcondition!(tile, tspan);
-prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:T,:θw,:θwi,:kc,:ρsn))
-sol = @time solve(prob, LiteImplicitEuler())
+prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600, savevars=(:T,:θw))
+sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600)
 out = CryoGridOutput(sol)
 
 # Plot the results!
@@ -54,20 +60,6 @@ cg = Plots.cgrad(:copper,rev=true);
 Plots.plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", title="", leg=false, dpi=150)
 Plots.plot(out.snowpack.dsn)
 
-integrator = init(prob, LiteImplicitEuler())
-for i in integrator
-    state = getstate(integrator)
-    last_above_ground_idx = findlast(<(0), cells(state.grid))
-    if state.snowpack.dsn[1] > 0.0 && state.snowpack.swe[last_above_ground_idx] <= 0 && state.snowpack.snowfall[1] <= 0
-        println("stopping")
-        break
-    end
-end
-
-CryoGrid.debug(true)
-
-step!(integrator)
-
 # CryoGridLite can also be embedded into integrators from OrdinaryDiffEq.jl via the `NLCGLite` nonlinear solver interface.
 # Note that these sovers generally will not be faster (in execution time) but may be more stable in some cases. Adaptive timestepping can be employed by
 # removing the `adaptive=false` argument.
diff --git a/examples/heat_sfcc_richardseq_samoylov.jl b/examples/heat_sfcc_richardseq_samoylov.jl
index 877d07f80adc593d20e8f6247ad07a6fe0fc1214..a2e57485e8bc70c946abbbf3865b7c3aa0ca4cf2 100644
--- a/examples/heat_sfcc_richardseq_samoylov.jl
+++ b/examples/heat_sfcc_richardseq_samoylov.jl
@@ -5,6 +5,8 @@
 # pressure gradients.
 
 using CryoGrid
+
+
 forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
 tspan = (DateTime(2011,1,1),DateTime(2011,12,31))
 T0 = forcings.Tair(tspan[1])
diff --git a/src/IO/InputOutput.jl b/src/IO/InputOutput.jl
index b055d861aab31c6a0201e7e8ed4275193ee83493..aeb237e906dce23a1a8d2b3e89f4b94aed77b87c 100644
--- a/src/IO/InputOutput.jl
+++ b/src/IO/InputOutput.jl
@@ -67,7 +67,7 @@ export Input, InputProvider, InputFunctionProvider
 export inputs
 include("input.jl")
 
-export ForcingFormat, ForcingFormatJSON, ForcingFormatNCD
+export ForcingFormat, ForcingFormatJSON, ForcingFormatNCD, Interpolated1D
 export loadforcings
 include("forcings/forcings.jl")
 
diff --git a/src/IO/forcings/providers.jl b/src/IO/forcings/providers.jl
index dfbaf81bb8169a32e7bd08c55a6b7185f17728ac..0d3cde62d4cb57846aac157f945a4f92f9da1a8d 100644
--- a/src/IO/forcings/providers.jl
+++ b/src/IO/forcings/providers.jl
@@ -41,6 +41,11 @@ function Base.show(io::IO, mimetype::MIME"text/plain", data::Interpolated1D)
     show(io, mimetype, getfield(data, :data))
 end
 
+function DimensionalData.rebuild(data::Interpolated1D, default_interp_mode=Linear(); stripunits=true, interp_modes=(;), datasets...)
+    newdata = DimStack((; datasets...))
+    return Interpolated1D(newdata, default_interp_mode; stripunits, interp_modes...)
+end
+
 function inputmap(data::Interpolated1D)
     return map(_wrap_interpolant, getfield(data, :interpolants))
 end
diff --git a/src/Physics/Heat/heat_water.jl b/src/Physics/Heat/heat_water.jl
index 18a403e5cab4a767e1fb9d32293db568a791e043..abdc2cfc772c5bc5096e4337436b25aedfe60694 100644
--- a/src/Physics/Heat/heat_water.jl
+++ b/src/Physics/Heat/heat_water.jl
@@ -18,7 +18,6 @@ advectiveflux(jw, T₁, T₂, cw, L) = jw*(cw*T₁*(jw > zero(jw)) + cw*T₂*(jw
 Adds advective energy fluxes for all internal grid cell faces.
 """
 function water_energy_advection!(jH, jw, T, cw::Real, L::Real)
-    water, heat = ps
     @inbounds for i in 2:length(jw)-1
         let jw = jw[i],
             T₁ = T[i-1],
@@ -86,18 +85,23 @@ function CryoGrid.interact!(
     end
 end
 
-# Flux calculation
 function CryoGrid.computeprognostic!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state)
     water, heat = ps
-    cw = heatcapacitywater(sub, state)
-    L = heat.prop.L
     CryoGrid.computeprognostic!(sub, water, state)
     if heat.advection
+        cw = heatcapacitywater(sub, state)
+        L = heat.prop.L
         water_energy_advection!(state.jH, state.jw, state.T, cw, L)
     end
     CryoGrid.computeprognostic!(sub, heat, state)
 end
 
+function CryoGrid.computeprognostic!(sub::SubSurface, ps::Coupled(WaterBalance{NoFlow}, HeatBalance), state)
+    water, heat = ps
+    CryoGrid.computeprognostic!(sub, water, state)
+    CryoGrid.computeprognostic!(sub, heat, state)
+end
+
 function heatcapacitywater(sub::SubSurface, state)
     if hasproperty(state, :ch_w)
         # TODO: this is a temporary solution that works so long as the
diff --git a/src/Tiles/stratigraphy.jl b/src/Tiles/stratigraphy.jl
index e7fab3649dbdfbb69288c3cd727bdae78a642c72..2d2d65e35f28ffd996755793cd59f4caa0e5064f 100644
--- a/src/Tiles/stratigraphy.jl
+++ b/src/Tiles/stratigraphy.jl
@@ -93,7 +93,9 @@ Base.NamedTuple(strat::Stratigraphy) = layers(strat)
 function Base.show(io::IO, ::MIME"text/plain", strat::Stratigraphy)
     print(io, "Stratigraphy:\n")
     for (k,v,b) in zip(keys(strat), values(strat), boundaries(strat))
-        print(repeat(" ", Base.indent_width), "$b: $k :: $(nameof(typeof(v)))\n")
+        procs = isa(processes(v), Tuple) ? processes(v) : (processes(v),)
+        procnames = map(p -> string(nameof(typeof(p))), procs)
+        print(repeat(" ", Base.indent_width), "$(b*u"m"): $k :: $(nameof(typeof(v))) with processes $(join(procnames, ", "))\n")
     end
 end
 
diff --git a/src/Tiles/tile.jl b/src/Tiles/tile.jl
index 3b57f17d884aac90853f5506faaa0803f17d91b4..96f2cd4fd919e446ba24974978562b4faf0fcb01 100644
--- a/src/Tiles/tile.jl
+++ b/src/Tiles/tile.jl
@@ -34,7 +34,13 @@ Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:inits}}) = true
 Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:events}}) = true
 # all other fields are non-flattenable by default
 Flatten.flattenable(::Type{<:Tile}, ::Type{Val{name}}) where name = false
-Base.show(io::IO, ::MIME"text/plain", tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,iip}) where {TStrat,TGrid,TStates,TInits,TEvents,iip} = print(io, "Tile (iip=$iip) with layers $(keys(tile.strat)), $TGrid, $TStrat")
+
+function Base.show(io::IO, mime::MIME"text/plain", tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip}) where {TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip}
+    println(io, "Tile (iip=$iip) with layers $(keys(tile.strat))")
+    print("  "); show(io, mime, tile.grid); println();
+    print("  "); show(io, mime, tile.strat);
+    println("  metadata: $(tile.metadata)")
+end
 
 """
     Tile(