Merge with changes coming from main branch

Author: daverumph

src/prognostic_equations/implicit/manual_sparse_jacobian.jl

@@ -619,39 +619,44 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
         FT = eltype(params)
         (; vertical_diffusion, smagorinsky_lilly) = p.atmos
         (; ᶜp) = p.precomputed
+        ᶜK_u = p.scratch.ᶜtemp_scalar_4
+        ᶜK_h = p.scratch.ᶜtemp_scalar_5
         if vertical_diffusion isa DecayWithHeightDiffusion
-            ᶜK_h = ᶜcompute_eddy_diffusivity_coefficient(Y.c.ρ, vertical_diffusion)
-            ᶜK_u = ᶜK_h
+            ᶜK_h .= ᶜcompute_eddy_diffusivity_coefficient(Y.c.ρ, vertical_diffusion)
+            ᶜK_u .= ᶜK_h
         elseif vertical_diffusion isa VerticalDiffusion
-            ᶜK_h = ᶜcompute_eddy_diffusivity_coefficient(Y.c.uₕ, ᶜp, vertical_diffusion)
-            ᶜK_u = ᶜK_h
+            ᶜK_h .= ᶜcompute_eddy_diffusivity_coefficient(Y.c.uₕ, ᶜp, vertical_diffusion)
+            ᶜK_u .= ᶜK_h
         elseif is_smagorinsky_vertical(smagorinsky_lilly)
             set_smagorinsky_lilly_precomputed_quantities!(Y, p, smagorinsky_lilly)
-            ᶜK_u = p.precomputed.ᶜνₜ_v
-            ᶜK_h = p.precomputed.ᶜD_v
+            ᶜK_u .= p.precomputed.ᶜνₜ_v
+            ᶜK_h .= p.precomputed.ᶜD_v
         elseif turbconv_model isa AbstractEDMF
             (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed
             ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
             ᶜmixing_length_field = p.scratch.ᶜtemp_scalar_3
             ᶜmixing_length_field .= ᶜmixing_length(Y, p)
-            ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field))
+            @. ᶜK_u = eddy_viscosity(turbconv_params, ᶜtke⁰, ᶜmixing_length_field)
+
             ᶜprandtl_nvec = @. lazy(
                 turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm),
             )
-            ᶜK_h = @. lazy(eddy_diffusivity(ᶜK_u, ᶜprandtl_nvec))
+            @. ᶜK_h = eddy_diffusivity(ᶜK_u, ᶜprandtl_nvec)
         end
 
+        ᶠρK_h = p.scratch.ᶠtemp_scalar
+        @. ᶠρK_h = ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_h)
         @. ᶜdiffusion_h_matrix =
-            ᶜadvdivᵥ_matrix() ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_h)) ⋅
-            ᶠgradᵥ_matrix()
+            ᶜadvdivᵥ_matrix() ⋅ DiagonalMatrixRow(ᶠρK_h) ⋅ ᶠgradᵥ_matrix()
         if (
             MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke)) ||
             !isnothing(p.atmos.turbconv_model) ||
             !disable_momentum_vertical_diffusion(p.atmos.vertical_diffusion)
         )
+            ᶠρK_u = p.scratch.ᶠtemp_scalar
+            @. ᶠρK_u = ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_u)
             @. ᶜdiffusion_u_matrix =
-                ᶜadvdivᵥ_matrix() ⋅
-                DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_u)) ⋅ ᶠgradᵥ_matrix()
+                ᶜadvdivᵥ_matrix() ⋅ DiagonalMatrixRow(ᶠρK_u) ⋅ ᶠgradᵥ_matrix()
         end
 
         ∂ᶜρe_tot_err_∂ᶜρ = matrix[@name(c.ρe_tot), @name(c.ρ)]