Add flux solver with manually unrolled iterations

src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl

@@ -71,6 +71,9 @@ include("tabulated_albedo.jl")
 # Turbulent fluxes
 include("roughness_lengths.jl")
 include("interface_states.jl")
+
+abstract type AbstractTurbulentFluxFormulation{S} end
+
 include("compute_interface_state.jl")
 include("similarity_theory_turbulent_fluxes.jl")
 include("coefficient_based_turbulent_fluxes.jl")

src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl

@@ -7,7 +7,7 @@ A structure for computing turbulent fluxes using constant bulk transfer coeffici
 
 $(TYPEDFIELDS)
 """
-struct CoefficientBasedFluxes{CD, CH, CQ, S}
+struct CoefficientBasedFluxes{CD, CH, CQ, S} <: AbstractTurbulentFluxFormulation{S}
     "Coefficient for momentum transfer"
     drag_coefficient :: CD
     "Coefficient for sensible heat transfer"

src/EarthSystemModels/InterfaceComputations/compute_interface_state.jl

@@ -51,12 +51,49 @@ end
                                       interface_properties,
                                       atmosphere_properties,
                                       interior_properties)
+
         iteration += 1
     end
 
     return Ψₛⁿ
 end
 
+#####
+##### Solver with unrolled iterations
+#####
+
+struct TenUnrolledIterations end
+
+@inline function compute_interface_state(flux_formulation::AbstractTurbulentFluxFormulation{<:TenUnrolledIterations},
+                                         initial_interface_state,
+                                         atmosphere_state,
+                                         interior_state,
+                                         downwelling_radiation,
+                                         interface_properties,
+                                         atmosphere_properties,
+                                         interior_properties)
+
+    args = (downwelling_radiation,
+            interface_properties,
+            atmosphere_properties,
+            interior_properties)
+
+    Ψₐ = atmosphere_state
+    Ψᵢ = interior_state
+    Ψₛ⁰ = initial_interface_state
+    Ψₛ¹ = iterate_interface_state(flux_formulation, Ψₛ⁰, Ψₐ, Ψᵢ, args...)
+    Ψₛ² = iterate_interface_state(flux_formulation, Ψₛ¹, Ψₐ, Ψᵢ, args...)
+    Ψₛ³ = iterate_interface_state(flux_formulation, Ψₛ², Ψₐ, Ψᵢ, args...)
+    Ψₛ⁴ = iterate_interface_state(flux_formulation, Ψₛ³, Ψₐ, Ψᵢ, args...)
+    Ψₛ⁵ = iterate_interface_state(flux_formulation, Ψₛ⁴, Ψₐ, Ψᵢ, args...)
+    Ψₛ⁶ = iterate_interface_state(flux_formulation, Ψₛ⁵, Ψₐ, Ψᵢ, args...)
+    Ψₛ⁷ = iterate_interface_state(flux_formulation, Ψₛ⁶, Ψₐ, Ψᵢ, args...)
+    Ψₛ⁸ = iterate_interface_state(flux_formulation, Ψₛ⁷, Ψₐ, Ψᵢ, args...)
+    Ψₛ⁹ = iterate_interface_state(flux_formulation, Ψₛ⁸, Ψₐ, Ψᵢ, args...)
+
+    return iterate_interface_state(flux_formulation, Ψₛ⁹, Ψₐ, Ψᵢ, args...)
+end
+
 """
     iterate_interface_state(flux_formulation, Ψₛⁿ⁻¹, Ψₐ, Ψᵢ, Qᵣ, ℙₛ, ℙₐ, ℙᵢ)
 

src/EarthSystemModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl

@@ -13,7 +13,7 @@ import Thermodynamics.Parameters: Rv_over_Rd
 ##### Bulk turbulent fluxes based on similarity theory
 #####
 
-struct SimilarityTheoryFluxes{FT, UF, R, B, S}
+struct SimilarityTheoryFluxes{FT, UF, R, B, S} <: AbstractTurbulentFluxFormulation{S}
     von_karman_constant :: FT        # parameter
     turbulent_prandtl_number :: FT   # parameter
     gustiness_parameter :: FT        # bulk velocity parameter