Optimize freeW freeze curve and related functions

src/Physics/HeatConduction/HeatConduction.jl

@@ -11,6 +11,7 @@ using CryoGrid.Numerics: nonlineardiffusion!, harmonicmean!, harmonicmean, heavi
 using CryoGrid.Layers: Soil, porosity
 using CryoGrid.Utils
 
+using Base: @propagate_inbounds
 using DimensionalData
 using IfElse
 using Interpolations: Linear, Flat

src/Physics/HeatConduction/heat.jl

@@ -201,8 +201,14 @@ total water content (θw), and liquid water content (θl).
             I_c*(H/Lθ) + I_t
         end
     end
-    @inbounds @. state.θl = freezethaw(state.H, heat.L, state.θw)*state.θw
-    heatcapacity!(layer, heat, state) # update heat capacity, C
-    @inbounds @. state.T = enthalpyinv(state.H, state.C, heat.L, state.θw)
+    @inbounds for i in 1:length(state.H)
+        let θw = state.θw[i],
+            H = state.H[i],
+            L = heat.L;
+            state.θl[i] = θw*freezethaw(H, L, θw)
+            state.C[i] = heatcapacity(layer, heat, state, i) # update heat capacity, C
+            state.T[i] = enthalpyinv(H, state.C[i], L, θw)
+        end
+    end
     return nothing
 end

src/Physics/HeatConduction/heat_bc.jl

@@ -12,7 +12,7 @@ struct TemperatureGradient{E,F} <: BoundaryProcess
 end
 BoundaryStyle(::Type{<:TemperatureGradient}) = Dirichlet()
 BoundaryStyle(::Type{<:TemperatureGradient{<:Damping}}) = Neumann()
-@inline boundaryvalue(bc::TemperatureGradient, l1, p2, l2, s1, s2) where {F} = bc.T(s1.t)
+@inline boundaryvalue(bc::TemperatureGradient, l1, p2, l2, s1, s2) = bc.T(s1.t)
 
 @with_kw struct NFactor{W,S} <: BoundaryEffect
     winterfactor::W = Param(1.0, bounds=(0.0,1.0)) # applied when Tair <= 0

src/Physics/HeatConduction/soil/soilheat.jl

+# Thermal properties
 @inline function heatcapacity(soil::Soil, totalwater, liquidwater, mineral, organic)
     @unpack cw,co,cm,ca,ci = soil.hc
     let air = 1.0 - totalwater - mineral - organic,
@@ -14,13 +15,16 @@ end
         (liq*kw^0.5 + ice*ki^0.5 + mineral*km^0.5 + organic*ko^0.5 + air*ka^0.5)^2
     end
 end
-@inline function heatcapacity!(soil::Soil, ::Heat, state)
-    @inbounds @. state.C = heatcapacity(soil, state.θw, state.θl, state.θm, state.θo)
+@inline @propagate_inbounds heatcapacity(soil::Soil, ::Heat, state, i) = heatcapacity(soil, state.θw[i], state.θl[i], state.θm[i], state.θo[i])
+@inline @propagate_inbounds function heatcapacity!(soil::Soil, ::Heat, state)
+    @. state.C = heatcapacity(soil, state.θw, state.θl, state.θm, state.θo)
 end
-@inline function thermalconductivity!(soil::Soil, ::Heat, state)
-    @inbounds @. state.kc = thermalconductivity(soil, state.θw, state.θl, state.θm, state.θo)
+@inline @propagate_inbounds thermalconductivity(soil::Soil, ::Heat, state, i) = thermalconductivity(soil, state.θw[i], state.θl[i], state.θm[i], state.θo[i])
+@inline @propagate_inbounds function thermalconductivity!(soil::Soil, ::Heat, state)
+    @. state.kc = thermalconductivity(soil, state.θw, state.θl, state.θm, state.θo)
 end
 
+# SFCC
 include("sfcc.jl")
 
 """