diff --git a/src/boundary_layer.jl b/src/boundary_layer.jl
index 295fbc4..51c07bc 100644
--- a/src/boundary_layer.jl
+++ b/src/boundary_layer.jl
@@ -1,10 +1,11 @@
 function allocate_profile(heights)
     wind_speed = similar(heights, typeof(0.0u"m/s")) # output wind speeds
     height_array = similar(heights, typeof(0.0u"m"))
-    height_array[end:-1:begin] .= heights 
+    height_array[end:-1:begin] .= heights
     air_temperature = similar(heights, typeof(0.0u"K")) # output temperatures, need to do this otherwise get InexactError
     relative_humidity = similar(heights, Float64) # output relative humidities
-    return (; heights, height_array, air_temperature, wind_speed, relative_humidity)
+    obukhov_length_prev = Ref(-0.3u"m")  # warm-start across timesteps
+    return (; heights, height_array, air_temperature, wind_speed, relative_humidity, obukhov_length_prev)
 end
 
 """
@@ -91,7 +92,7 @@ function atmospheric_surface_profile!(buffers;
     (; roughness_height, karman_constant, dyer_constant, elevation) = micro_terrain
     (; atmospheric_pressure, reference_temperature, reference_wind_speed, reference_humidity, zenith_angle) = environment_instant
 
-    (; heights, height_array, air_temperature, wind_speed, relative_humidity) = buffers
+    (; heights, height_array, air_temperature, wind_speed, relative_humidity, obukhov_length_prev) = buffers
     N_heights = length(heights)
     if minimum(heights) < roughness_height
         throw(ArgumentError("The minimum height is not greater than the roughness height."))
@@ -138,10 +139,9 @@ function atmospheric_surface_profile!(buffers;
             air_temperature[i] = roughness_height_temp + (reference_temp - roughness_height_temp) * log(height_array[i] / z0 + 1.0) / log_z_ratio
         end
     else
-        obukhov_length = -0.3u"m" # initialise Obukhov length
-        # TODO just pass the environment_instant through here
-        Obukhov_out = calc_Obukhov_length(reference_temp, surface_temp, v_ref_height, z0, z, ρcpTκg, κ, log_z_ratio, ΔT, ρ_cp; max_iter=30, tol=1e-2)
+        Obukhov_out = calc_Obukhov_length(reference_temp, surface_temp, v_ref_height, z0, z, ρcpTκg, κ, log_z_ratio, ΔT, ρ_cp; max_iter=30, tol=1e-2, initial_obukhov_length=obukhov_length_prev[])
         obukhov_length = Obukhov_out.obukhov_length
+        obukhov_length_prev[] = obukhov_length
         roughness_height_temp = Obukhov_out.roughness_height_temperature
         convective_heat_flux = Obukhov_out.convective_heat_flux
         u_star = Obukhov_out.u_star
@@ -415,9 +415,9 @@ Iteratively solve for Monin-Obukhov length and convective heat flux.
 """
 @inline function calc_Obukhov_length(
     reference_temp, surface_temp, v_ref_height, z0, z, ρcpTκg, κ, log_z_ratio, ΔT, ρ_cp;
-    γ=16.0, max_iter=30, tol=1e-2
+    γ=16.0, max_iter=30, tol=1e-2, initial_obukhov_length=-0.3u"m"
 )
-    obukhov_length = -0.3u"m" # initial Monin-Obukhov length
+    obukhov_length = initial_obukhov_length
 
     # initialise with zeros
     convective_heat_flux = 0.0u"W/m^2"
diff --git a/src/simulation.jl b/src/simulation.jl
index ecd7294..e2b0f1f 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -442,6 +442,7 @@ function solve_soil!(output::MicroResult, mp::MicroProblem, solar_radiation_out;
                         T0 = temperature
                     end
                 end
+                init_soil_obukhov!(buffers, forcing, micro_terrain, heights, T0, i)
                 rain = hourly_rainfall ? mp.environment_hourly.rainfall[step] : environment_instant.rainfall
                 pool = clamp(pool + rain, 0.0u"kg/m^2", soil_moisture_model.maxpool)
                 if runmoist
@@ -456,9 +457,8 @@ function solve_soil!(output::MicroResult, mp::MicroProblem, solar_radiation_out;
                     output.surface_water[step] = pool
                     output.soil_temperature[step, :] .= T0
                     output.sky_temperature[step] = longwave_sky.sky_temperature
-                                    environment_instant = get_instant(environment_day, mp.environment_hourly, output, soil_moisture, step)
-
-                                    update_soil_properties!(output, buffers.soil_properties, soil_thermal_model;
+                    environment_instant = get_instant(environment_day, mp.environment_hourly, output, soil_moisture, step)
+                    update_soil_properties!(output, buffers.soil_properties, soil_thermal_model;
                         soil_temperature=T0, soil_moisture, atmospheric_pressure=environment_instant.atmospheric_pressure, step, vapour_pressure_equation
                     )
                 end
diff --git a/src/soil_balance.jl b/src/soil_balance.jl
index 74d9b03..f43124e 100644
--- a/src/soil_balance.jl
+++ b/src/soil_balance.jl
@@ -2,7 +2,8 @@ function allocate_soil_energy_balance(num_nodes::Int)
     layer_depths = fill(0.0u"cm", num_nodes + 1)
     heat_capacity = fill(1.0u"J/K/m^2", num_nodes)
     thermal_conductance = fill(1.0u"W/K/m^2", num_nodes)
-    return (; layer_depths, heat_capacity, thermal_conductance)
+    obukhov_length_prev = Ref(-0.3u"m")  # warm-start across hourly ODE solves
+    return (; layer_depths, heat_capacity, thermal_conductance, obukhov_length_prev)
 end
 
 # This is a 3-parameters OrdinaryDiffEq function
@@ -70,9 +71,7 @@ function soil_energy_balance(
         u_star = calc_u_star(; reference_wind_speed=wind_speed, log_z_ratio, κ=karman_constant)
         convective_heat_flux = calc_convection(; u_star, log_z_ratio, ΔT, ρ_cp, z0=roughness_height)
     else
-        # compute ρcpTκg (was a constant in original Fortran version)
-        ρcpTκg = 6.003e-8u"cal*minute^2/cm^4"
-        Obukhov_out = calc_Obukhov_length(air_temperature, surface_temperature, wind_speed, roughness_height, reference_height, ρcpTκg, karman_constant, log_z_ratio, ΔT, ρ_cp)
+        Obukhov_out = calc_soil_obukhov(air_temperature, surface_temperature, wind_speed, roughness_height, reference_height, karman_constant; initial_obukhov_length=buffers.soil_energy_balance.obukhov_length_prev[])
         convective_heat_flux = Obukhov_out.convective_heat_flux
     end
     heat_transfer_coefficient = max(abs(convective_heat_flux / (soil_temperature[1] - air_temperature)), 0.5u"W/m^2/K")
@@ -108,6 +107,26 @@ function soil_energy_balance(
 end
 
 
+function calc_soil_obukhov(air_temperature, surface_temperature, wind_speed, roughness_height, reference_height, karman_constant; initial_obukhov_length)
+    log_z_ratio = log(reference_height / roughness_height + 1)
+    ΔT = air_temperature - surface_temperature
+    ρ_cp = calc_ρ_cp((surface_temperature + air_temperature) / 2)
+    ρcpTκg = 6.003e-8u"cal*minute^2/cm^4"
+    return calc_Obukhov_length(air_temperature, surface_temperature, wind_speed, roughness_height, reference_height, ρcpTκg, karman_constant, log_z_ratio, ΔT, ρ_cp; initial_obukhov_length)
+end
+
+function init_soil_obukhov!(buffers, forcing, micro_terrain, heights, T0, i)
+    t_next = ((i - 1) * 60)u"minute"
+    (; air_temperature, wind_speed, zenith_angle) = interpolate_forcings(forcing, t_next)
+    surface_temperature = T0[1]
+    if air_temperature < surface_temperature && zenith_angle < 90u"°"
+        (; roughness_height, karman_constant) = micro_terrain
+        reference_height = last(heights)
+        Obukhov_out = calc_soil_obukhov(air_temperature, surface_temperature, wind_speed, roughness_height, reference_height, karman_constant; initial_obukhov_length=buffers.soil_energy_balance.obukhov_length_prev[])
+        buffers.soil_energy_balance.obukhov_length_prev[] = Obukhov_out.obukhov_length
+    end
+end
+
 function interpolate_forcings(f, t)
     t_m = ustrip(u"minute", t)
     return (;