diff --git a/experiments/standalone/Vegetation/timestep_test.jl b/experiments/standalone/Vegetation/timestep_test.jl
index 782af58d66..6e68d9f362 100644
--- a/experiments/standalone/Vegetation/timestep_test.jl
+++ b/experiments/standalone/Vegetation/timestep_test.jl
@@ -150,71 +150,72 @@ canopy = ClimaLand.Canopy.CanopyModel{FT}(;
     radiation = radiation,
 );
 
-
-Y, p, coords = ClimaLand.initialize(canopy)
 exp_tendency! = make_exp_tendency(canopy)
 # imp_tendency! = make_imp_tendency(canopy)
 # jacobian! = make_jacobian(canopy);
 # jac_kwargs =
 #     (; jac_prototype = ClimaLand.ImplicitEquationJacobian(Y), Wfact = jacobian!);
 
-
-ψ_leaf_0 = FT(-2e5 / 9800)
-ψ_stem_0 = FT(-1e5 / 9800)
-
-S_l_ini =
-    inverse_water_retention_curve.(
-        retention_model,
-        [ψ_stem_0, ψ_leaf_0],
-        ν,
-        S_s,
-    )
-
-Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini[1])
-Y.canopy.hydraulics.ϑ_l.:2 .= augmented_liquid_fraction.(ν, S_l_ini[2])
-
-
 seconds_per_day = 3600 * 24.0
 t0 = 150seconds_per_day
 N_days = 0.5
 tf = t0 + N_days * seconds_per_day
-evaluate!(Y.canopy.energy.T, atmos.T, t0)
-set_initial_cache! = make_set_initial_cache(canopy)
-set_initial_cache!(p, Y, t0);
-
-saveat = Array(t0:(3 * 3600):tf)
-updateat = Array(t0:(3600 * 3):tf)
-drivers = ClimaLand.get_drivers(canopy)
-updatefunc = ClimaLand.make_update_drivers(drivers)
-cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
-
-# timestepper = CTS.ARS111();
-# ode_algo = CTS.IMEXAlgorithm(
-#     timestepper,
-#     CTS.NewtonsMethod(
-#         max_iters = 6,
-#         update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
-#     ),
-# );
-timestepper = CTS.RK4();
-ode_algo = CTS.ExplicitAlgorithm(timestepper)
-
-prob = SciMLBase.ODEProblem(
-    CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!),
-    Y,
-    (t0, tf),
-    p,
-);
 
 savedir = joinpath(pkgdir(ClimaLand), "experiments/standalone/Vegetation");
 
 ref_dt = 0.001
+
+############ Initial setup for ref solution
+# Y, p, coords = ClimaLand.initialize(canopy)
+
+# ψ_leaf_0 = FT(-2e5 / 9800)
+# ψ_stem_0 = FT(-1e5 / 9800)
+
+# S_l_ini =
+#     inverse_water_retention_curve.(
+#         retention_model,
+#         [ψ_stem_0, ψ_leaf_0],
+#         ν,
+#         S_s,
+#     )
+
+# Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini[1])
+# Y.canopy.hydraulics.ϑ_l.:2 .= augmented_liquid_fraction.(ν, S_l_ini[2])
+
+# evaluate!(Y.canopy.energy.T, atmos.T, t0)
+# set_initial_cache! = make_set_initial_cache(canopy)
+# set_initial_cache!(p, Y, t0);
+
+# saveat = Array(t0:(3 * 3600):tf)
+# updateat = Array(t0:(3600 * 3):tf)
+# drivers = ClimaLand.get_drivers(canopy)
+# updatefunc = ClimaLand.make_update_drivers(drivers)
+# cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
+
+# # timestepper = CTS.ARS111();
+# # ode_algo = CTS.IMEXAlgorithm(
+# #     timestepper,
+# #     CTS.NewtonsMethod(
+# #         max_iters = 6,
+# #         update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
+# #     ),
+# # );
+# timestepper = CTS.RK4();
+# ode_algo = CTS.ExplicitAlgorithm(timestepper)
+
+# prob = SciMLBase.ODEProblem(
+#     CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!),
+#     Y,
+#     (t0, tf),
+#     p,
+# );
 # ref_sol =
 #     SciMLBase.solve(prob, ode_algo; dt = ref_dt, callback = cb, saveat = saveat);
 # ref_T = [parent(ref_sol.u[k].canopy.energy.T)[1] for k in 1:length(ref_sol.t)]
 # open(joinpath(savedir, "exp_T_dt$(ref_dt)_$(N_days)days.txt"), "w") do io
 #     writedlm(io, ref_T, ',')
 # end;
+###############
 
 ref_T = readdlm(joinpath(savedir, "exp_T_dt$(ref_dt)_$(N_days)days.txt"))
 
@@ -223,14 +224,39 @@ sols = []
 for dt in dts
     @info dt
 
+    # Set initial values of Y
     Y, p, coords = ClimaLand.initialize(canopy)
+    ψ_leaf_0 = FT(-2e5 / 9800)
+    ψ_stem_0 = FT(-1e5 / 9800)
+
+    S_l_ini =
+        inverse_water_retention_curve.(
+            retention_model,
+            [ψ_stem_0, ψ_leaf_0],
+            ν,
+            S_s,
+        )
+
+    Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini[1])
+    Y.canopy.hydraulics.ϑ_l.:2 .= augmented_liquid_fraction.(ν, S_l_ini[2])
+
+    # Set initial values of atmos driver and cache
     evaluate!(Y.canopy.energy.T, atmos.T, t0)
-    set_initial_cache!(p, Y, t0);
+    set_initial_cache!(p, Y, t0)
+
     saveat = Array(t0:(3 * 3600):tf)
     updateat = Array(t0:(3600 * 3):tf)
     updatefunc = ClimaLand.make_update_drivers(drivers)
     cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
 
+    # Construct problem using newly-initialized values
+    prob = SciMLBase.ODEProblem(
+        CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!),
+        Y,
+        (t0, tf),
+        p,
+    )
+
     @time sol =
         SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat)
     T = [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)]