Review comments part I

src/integrated/land.jl

@@ -481,10 +481,20 @@ function soil_boundary_fluxes!(
 ) where {FT}
     bc = soil.boundary_conditions.top
     p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, soil, Y, p, t)
-    influx = @. p.drivers.P_liq +
-       p.snow.water_runoff * p.snow.snow_cover_fraction +
-       p.excess_water_flux +
-       (1 - p.snow.snow_cover_fraction) * p.soil.turbulent_fluxes.vapor_flux_liq
+    # influx = maximum possible rate of infiltration given precip, snowmelt, evaporation/condensation
+    # but if this exceeds infiltration capacity of the soil, runoff will
+    # be generated.
+    # Use top_bc.water as temporary storage to avoid allocation
+    influx = p.soil.top_bc.water
+    @. influx =
+        p.drivers.P_liq +
+        p.snow.water_runoff * p.snow.snow_cover_fraction +
+        p.excess_water_flux +
+        (1 - p.snow.snow_cover_fraction) *
+        p.soil.turbulent_fluxes.vapor_flux_liq
+    # The update_runoff! function computes how much actually infiltrates
+    # given influx and our runoff model bc.runoff, and updates
+    # p.soil.infiltration in place
     Soil.Runoff.update_runoff!(p, bc.runoff, influx, Y, t, soil)
     @. p.soil.top_bc.water = p.soil.infiltration
     @. p.soil.top_bc.heat =

src/integrated/soil_canopy_model.jl

@@ -370,14 +370,16 @@ function soil_boundary_fluxes!(
 ) where {FT}
     bc = soil.boundary_conditions.top
     p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, soil, Y, p, t)
-    Soil.Runoff.update_runoff!(
-        p,
-        bc.runoff,
-        p.drivers.P_liq .+ p.soil.turbulent_fluxes.vapor_flux_liq,
-        Y,
-        t,
-        soil,
-    )
+    # influx = maximum possible rate of infiltration given precip, snowmelt, evaporation/condensation
+    # but if this exceeds infiltration capacity of the soil, runoff will
+    # be generated.
+    # Use top_bc.water as temporary storage to avoid allocation
+    influx = p.soil.top_bc.water
+    @. influx = p.drivers.P_liq + p.soil.turbulent_fluxes.vapor_flux_liq
+    # The update_runoff! function computes how much actually infiltrates
+    # given influx and our runoff model bc.runoff, and updates
+    # p.soil.infiltration in place
+    Soil.Runoff.update_runoff!(p, bc.runoff, influx, Y, t, soil)
     @. p.soil.top_bc.water = p.soil.infiltration
     @. p.soil.top_bc.heat =
         -p.soil.R_n + p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf

src/integrated/soil_snow_model.jl

@@ -387,10 +387,20 @@ function soil_boundary_fluxes!(
     bc = soil.boundary_conditions.top
     p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, soil, Y, p, t)
     p.soil.R_n .= net_radiation(bc.radiation, soil, Y, p, t)
-    influx = @. p.drivers.P_liq +
-       p.snow.water_runoff * p.snow.snow_cover_fraction +
-       p.excess_water_flux +
-       (1 - p.snow.snow_cover_fraction) * p.soil.turbulent_fluxes.vapor_flux_liq
+    # influx = maximum possible rate of infiltration given precip, snowmelt, evaporation/condensation
+    # but if this exceeds infiltration capacity of the soil, runoff will
+    # be generated.
+    # Use top_bc.water as temporary storage to avoid allocation
+    influx = p.soil.top_bc.water
+    @. influx =
+        p.drivers.P_liq +
+        p.snow.water_runoff * p.snow.snow_cover_fraction +
+        p.excess_water_flux +
+        (1 - p.snow.snow_cover_fraction) *
+        p.soil.turbulent_fluxes.vapor_flux_liq
+    # The update_runoff! function computes how much actually infiltrates
+    # given influx and our runoff model bc.runoff, and updates
+    # p.soil.infiltration in place
     Soil.Runoff.update_runoff!(p, bc.runoff, influx, Y, t, soil)
     @. p.soil.top_bc.water = p.soil.infiltration
 

src/standalone/Snow/snow_parameterizations.jl

@@ -164,6 +164,12 @@ Computes the thermal conductivity, given the density
 of the snow, according to Equation
 5.33 from Bonan's textbook, which in turn is taken from
 Jordan (1991).
+
+We have adjusted the original equation to make the coefficients
+non-dimensional by multiplying by the first by x = ρ_ice/ρ_ice
+and the second by x², with ρ_ice in kg/m³.
+
+When ρ_snow = ρ_ice, we recover κ_snow = κ_ice.
 """
 function snow_thermal_conductivity(
     ρ_snow::FT,

src/standalone/Soil/boundary_conditions.jl

@@ -825,14 +825,16 @@ function soil_boundary_fluxes!(
 )
     p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, model, Y, p, t)
     p.soil.R_n .= net_radiation(bc.radiation, model, Y, p, t)
-    update_runoff!(
-        p,
-        bc.runoff,
-        p.drivers.P_liq .+ p.soil.turbulent_fluxes.vapor_flux_liq,
-        Y,
-        t,
-        model,
-    )
+    # influx = maximum possible rate of infiltration given precip, snowmelt, evaporation/condensation
+    # but if this exceeds infiltration capacity of the soil, runoff will
+    # be generated.
+    # Use top_bc.water as temporary storage to avoid allocation
+    influx = p.soil.top_bc.water
+    @. influx = p.drivers.P_liq + p.soil.turbulent_fluxes.vapor_flux_liq
+    # The update_runoff! function computes how much actually infiltrates
+    # given influx and our runoff model bc.runoff, and updates
+    # p.soil.infiltration in place
+    update_runoff!(p, bc.runoff, influx, Y, t, model)
     # We do not model the energy flux from infiltration.
     @. p.soil.top_bc.water = p.soil.infiltration
     @. p.soil.top_bc.heat =