CliMA · simone-silvestri · Feb 25, 2025 · Feb 25, 2025 · Feb 25, 2025 · Feb 25, 2025
diff --git a/docs/make.jl b/docs/make.jl
@@ -17,7 +17,7 @@ example_scripts = [
     "melting_in_spring.jl",
     "freezing_of_a_lake.jl",
     "ice_advected_by_anticyclone.jl",
-    # "ice_advected_on_coastline.jl",
+    "ice_advected_on_coastline.jl",
     "arctic_basin_seasonal_cycle.jl"
 ]
 
@@ -31,7 +31,7 @@ example_pages = [
     "Melting in Spring" => "literated/melting_in_spring.md",
     "Freezing of a Lake" => "literated/freezing_of_a_lake.md",
     "Ice advected by anticyclone" => "literated/ice_advected_by_anticyclone.md",
-    # "Ice advected on coastline" => "literated/ice_advected_on_coastline.md",
+    "Ice advected on coastline" => "literated/ice_advected_on_coastline.md",
     "Arctic basin seasonal cycle" => "literated/arctic_basin_seasonal_cycle.md"
 ]
 

diff --git a/examples/freezing_bucket.jl b/examples/freezing_bucket.jl
@@ -53,9 +53,9 @@ top_heat_boundary_condition = PrescribedTemperature(-10)
 # slab sea ice representation of ice_thermodynamics
 
 ice_thermodynamics = SlabSeaIceThermodynamics(grid;
-                                          internal_heat_flux,
-                                          phase_transitions,
-                                          top_heat_boundary_condition)
+                                              internal_heat_flux,
+                                              phase_transitions,
+                                              top_heat_boundary_condition)
 
 # We also prescribe a frazil ice heat flux that stops when the ice has reached a concentration of 1.
 # This heat flux represents the initial ice formation from a liquid bucket.

diff --git a/examples/ice_advected_by_anticyclone.jl b/examples/ice_advected_by_anticyclone.jl
@@ -25,7 +25,7 @@ L  = 512kilometers
 
 # 2 km domain
 grid = RectilinearGrid(arch;
-                       size = (128, 128), 
+                       size = (256, 256), 
                           x = (0, L), 
                           y = (0, L), 
                        halo = (7, 7),
@@ -35,19 +35,23 @@ grid = RectilinearGrid(arch;
 ##### Value boundary conditions for velocities
 #####
 
-u_bcs = FieldBoundaryConditions(north=ValueBoundaryCondition(0),
-                                south=ValueBoundaryCondition(0))
+u_bcs = FieldBoundaryConditions(north=OpenBoundaryCondition(0),
+                                south=OpenBoundaryCondition(0),
+                                west=OpenBoundaryCondition(0),
+                                east=OpenBoundaryCondition(0))
 
-v_bcs = FieldBoundaryConditions(west=ValueBoundaryCondition(0),
-                                east=ValueBoundaryCondition(0))
+v_bcs = FieldBoundaryConditions(north=OpenBoundaryCondition(0),
+                                south=OpenBoundaryCondition(0),
+                                west=OpenBoundaryCondition(0),
+                                east=OpenBoundaryCondition(0))
 
 #####
 ##### Ocean sea-ice stress
 #####
 
 # Constant ocean velocities corresponding to a cyclonic eddy
-Uₒ = XFaceField(grid)
-Vₒ = YFaceField(grid)
+Uₒ = Field{Face, Face, Nothing}(grid)
+Vₒ = Field{Face, Face, Nothing}(grid)
 
 set!(Uₒ, (x, y) -> 𝓋ₒ * (2y - L) / L)
 set!(Vₒ, (x, y) -> 𝓋ₒ * (L - 2x) / L)
@@ -59,8 +63,8 @@ fill_halo_regions!((Uₒ, Vₒ))
 #### Atmosphere - sea ice stress 
 ####
 
-Uₐ = XFaceField(grid)
-Vₐ = YFaceField(grid)
+Uₐ = Field{Face, Face, Nothing}(grid)
+Vₐ = Field{Face, Face, Nothing}(grid)
 
 # Atmospheric velocities corresponding to an anticyclonic eddy moving north-east
 @inline center(t) = 256kilometers + 51.2kilometers * t / 86400

diff --git a/examples/ice_advected_on_coastline.jl b/examples/ice_advected_on_coastline.jl
@@ -64,15 +64,10 @@ dynamics = SeaIceMomentumEquation(grid;
                                   rheology = ElastoViscoPlasticRheology(),
                                   solver = SplitExplicitSolver(substeps=150))
 
-u_bcs = FieldBoundaryConditions(top = nothing, bottom = nothing,
-                                north = ValueBoundaryCondition(0),
-                                south = ValueBoundaryCondition(0))
-
 #Define the model! 
 model = SeaIceModel(grid; 
                     advection = WENO(order=7),
                     dynamics = dynamics,
-                    boundary_conditions = (; u=u_bcs),
                     ice_thermodynamics = nothing)
 
 # We start with a concentration of ℵ = 1 everywhere

diff --git a/src/Rheologies/Rheologies.jl b/src/Rheologies/Rheologies.jl
@@ -1,25 +1,29 @@
 module Rheologies
 
 export ViscousRheology, ElastoViscoPlasticRheology
-export ∂ⱼ_σ₁ⱼ, ∂ⱼ_σ₂ⱼ, required_auxiliary_fields
+export ∂ⱼ_σ₁ⱼ, ∂ⱼ_σ₂ⱼ, rheology_auxiliary_fields
 
 using Oceananigans
 using Oceananigans.Operators
 using ClimaSeaIce: ice_mass
 
 # Nothing rheology
-required_auxiliary_fields(rheology, grid) = NamedTuple()
 initialize_rheology!(model, rheology) = nothing
-compute_stresses!(model, dynamics, rheology, Δt) = nothing
+compute_stresses!(model, dynamics, rheology, Δt, Ns) = nothing
 
 # Nothing rheology or viscous rheology
-@inline compute_substep_Δtᶠᶜᶜ(i, j, grid, Δt, rheology, substeps, fields) = Δt / substeps
-@inline compute_substep_Δtᶜᶠᶜ(i, j, grid, Δt, rheology, substeps, fields) = Δt / substeps
+@inline compute_substep_Δtᶠᶠᶜ(i, j, grid, Δt, rheology, substeps, fields) = Δt / substeps
 
 # Fallback
 @inline sum_of_forcing_u(i, j, k, grid, rheology, u_forcing, fields, Δt) = u_forcing(i, j, k, grid, fields)
 @inline sum_of_forcing_v(i, j, k, grid, rheology, v_forcing, fields, Δt) = v_forcing(i, j, k, grid, fields)
 
+# No needed extra tracers
+rheology_prognostic_tracers(rheology) = ()
+
+# No needed auxiliary fields
+rheology_auxiliary_fields(rheology, grid) = NamedTuple()
+
 @inline ice_stress_ux(i, j, k, grid, ::Nothing, args...) = zero(grid)
 @inline ice_stress_uy(i, j, k, grid, ::Nothing, args...) = zero(grid)
 @inline ice_stress_vx(i, j, k, grid, ::Nothing, args...) = zero(grid)

diff --git a/src/Rheologies/elasto_visco_plastic_rheology.jl b/src/Rheologies/elasto_visco_plastic_rheology.jl
@@ -105,19 +105,20 @@ function ElastoViscoPlasticRheology(FT::DataType = Float64;
                                       pressure_formulation)
 end
 
-function required_auxiliary_fields(r::ElastoViscoPlasticRheology, grid)
+function rheology_auxiliary_fields(r::ElastoViscoPlasticRheology, grid)
 
     # TODO: What about boundary conditions?
     σ₁₁ = Field{Center, Center, Nothing}(grid)
     σ₂₂ = Field{Center, Center, Nothing}(grid)
-    σ₁₂ = Field{Face, Face, Nothing}(grid)
+    σ₁₂ = Field{Center, Center, Nothing}(grid)
 
-    uⁿ = Field{Face,   Center, Nothing}(grid)
-    vⁿ = Field{Center, Face,   Nothing}(grid)
     P  = Field{Center, Center, Nothing}(grid)
-    α  = Field{Center, Center, Nothing}(grid) # Dynamic substeps a la Kimmritz et al (2016)
     ζ  = Field{Center, Center, Nothing}(grid)
     Δ  = Field{Center, Center, Nothing}(grid)
+    α  = Field{Center, Center, Nothing}(grid) # Dynamic substeps a la Kimmritz et al (2016)
+
+    uⁿ = Field{Face, Face, Nothing}(grid)
+    vⁿ = Field{Face, Face, Nothing}(grid)
 
     # An initial (safe) educated guess
     fill!(α, r.max_relaxation_parameter)
@@ -170,7 +171,7 @@ end
 @inline ice_strength(i, j, k, grid, P★, C, h, ℵ) = @inbounds P★ * h[i, j, k] * exp(- C * (1 - ℵ[i, j, k])) 
 
 # Specific compute stresses for the EVP rheology
-function compute_stresses!(model, dynamics, rheology::ElastoViscoPlasticRheology, Δt) 
+function compute_stresses!(model, dynamics, rheology::ElastoViscoPlasticRheology, Δt, Ns) 
 
     grid = model.grid
     arch = architecture(grid)
@@ -187,17 +188,25 @@ function compute_stresses!(model, dynamics, rheology::ElastoViscoPlasticRheology
 
     parameters = KernelParameters(-Hx+2:Nx+Hx-1, -Hy+2:Ny+Hy-1)
 
-    launch!(arch, grid, parameters, _compute_evp_viscosities!, fields, grid, rheology, u, v)
-    launch!(arch, grid, parameters, _compute_evp_stresses!, fields, grid, rheology, u, v, h, ℵ, ρᵢ, Δt)
+    launch!(arch, grid, parameters, _compute_evp_viscosities!, fields, grid, rheology, u, v, h, ℵ, ρᵢ, Δt)
 
     return nothing
 end
 
-@inline strain_rate_xx(i, j, k, grid, u, v) =  δxᶜᵃᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u) / Azᶜᶜᶜ(i, j, k, grid)
-@inline strain_rate_yy(i, j, k, grid, u, v) =  δyᵃᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v) / Azᶜᶜᶜ(i, j, k, grid)
-@inline strain_rate_xy(i, j, k, grid, u, v) = (δxᶠᵃᵃ(i, j, k, grid, Δy_qᶜᶠᶜ, v) + δyᵃᶠᵃ(i, j, k, grid, Δx_qᶠᶜᶜ, u)) / Azᶠᶠᶜ(i, j, k, grid) / 2
+@inline strain_rate_xx(i, j, k, grid, u, v) =  ℑyᵃᶜᵃ(i, j, k, grid, δxᶜᵃᵃ, Δy_qᶠᶠᶜ, u) * Az⁻¹ᶜᶜᶜ(i, j, k, grid)
+@inline strain_rate_yy(i, j, k, grid, u, v) =  ℑxᶜᵃᵃ(i, j, k, grid, δyᵃᶜᵃ, Δx_qᶠᶠᶜ, v) * Az⁻¹ᶜᶜᶜ(i, j, k, grid)
+@inline strain_rate_xy(i, j, k, grid, u, v) = (δxᶜᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, Δy_qᶠᶠᶜ, v) + δyᵃᶜᵃ(i, j, k, grid, ℑxᶜᵃᵃ, Δx_qᶠᶠᶜ, u)) * Az⁻¹ᶜᶜᶜ(i, j, k, grid) / 2
+
+@inline ice_pressure(i, j, k, grid, ::IceStrength, r, fields) = @inbounds fields.P[i, j, k]
 
-@kernel function _compute_evp_viscosities!(fields, grid, rheology, u, v)
+@inline function ice_pressure(i, j, k, grid, ::ReplacementPressure, r, fields)
+    Pᶜᶜᶜ = @inbounds fields.P[i, j, k]
+    Δᶜᶜᶜ = @inbounds fields.Δ[i, j, k]
+    Δm   = r.minimum_plastic_stress
+    return Pᶜᶜᶜ * Δᶜᶜᶜ / (Δᶜᶜᶜ + Δm)
+end
+
+@kernel function _compute_evp_viscosities!(fields, grid, rheology, u, v, h, ℵ, ρᵢ, Δt)
     i, j = @index(Global, NTuple)
     kᴺ   = size(grid, 3)
 
@@ -210,41 +219,24 @@ end
     # Strain rates
     ϵ̇₁₁ = strain_rate_xx(i, j, kᴺ, grid, u, v) 
     ϵ̇₂₂ = strain_rate_yy(i, j, kᴺ, grid, u, v) 
-
-    # Center - Center variables:
-    ϵ̇₁₂ᶜᶜᶜ = ℑxyᶜᶜᵃ(i, j, kᴺ, grid, strain_rate_xy, u, v)
+    ϵ̇₁₂ = strain_rate_xy(i, j, kᴺ, grid, u, v)
 
     # Ice divergence 
     δ = ϵ̇₁₁ + ϵ̇₂₂
 
     # Ice shear (at Centers)
-    s = sqrt((ϵ̇₁₁ - ϵ̇₂₂)^2 + 4ϵ̇₁₂ᶜᶜᶜ^2)
+    s = sqrt((ϵ̇₁₁ - ϵ̇₂₂)^2 + 4ϵ̇₁₂^2)
 
     # Visco - Plastic parameter 
     # if Δ is very small we assume a linear viscous response
     # adding a minimum Δ_min (at Centers)
-    Δᶜᶜᶜ = max(sqrt(δ^2 + s^2 * e⁻²), Δm)
-    Pᶜᶜᶜ = @inbounds P[i, j, 1]
+    Δ = max(sqrt(δ^2 + s^2 * e⁻²), Δm) 
+    P = @inbounds P[i, j, kᴺ]
+    ζ = P / 2Δ
 
-    @inbounds fields.ζ[i, j, 1] = Pᶜᶜᶜ / 2Δᶜᶜᶜ
-    @inbounds fields.Δ[i, j, 1] = Δᶜᶜᶜ
-end
-
-@inline ice_pressure(i, j, k, grid, ::IceStrength, r, fields) = @inbounds fields.P[i, j, k]
-
-@inline function ice_pressure(i, j, k, grid, ::ReplacementPressure, r, fields)
-    Pᶜᶜᶜ = @inbounds fields.P[i, j, k]
-    Δᶜᶜᶜ = @inbounds fields.Δ[i, j, k]
-    Δm   = r.minimum_plastic_stress
-    return Pᶜᶜᶜ * Δᶜᶜᶜ / (Δᶜᶜᶜ + Δm)
-end
-
-# Compute the visco-plastic stresses for a slab sea ice model.
-# The function updates the internal stress variables `σ₁₁`, `σ₂₂`, and `σ₁₂` in the `rheology` object
-# following the αEVP formulation of Kimmritz et al (2016).
-@kernel function _compute_evp_stresses!(fields, grid, rheology, u, v, h, ℵ, ρᵢ, Δt)
-    i, j = @index(Global, NTuple)
-    kᴺ   = size(grid, 3)
+    # Store viscosity and the deformation for analysis purposes
+    @inbounds fields.ζ[i, j, kᴺ] = ζ
+    @inbounds fields.Δ[i, j, kᴺ] = Δ
 
     e⁻² = rheology.yield_curve_eccentricity^(-2)
     α⁺  = rheology.max_relaxation_parameter
@@ -257,50 +249,36 @@ end
     σ₁₂ = fields.σ₁₂
     α   = fields.α
 
-    # Strain rates
-    ϵ̇₁₁ = strain_rate_xx(i, j, kᴺ, grid, u, v) 
-    ϵ̇₂₂ = strain_rate_yy(i, j, kᴺ, grid, u, v) 
-    ϵ̇₁₂ = strain_rate_xy(i, j, kᴺ, grid, u, v)
-
-    ζᶜᶜᶜ = @inbounds fields.ζ[i, j, 1]
-    ζᶠᶠᶜ = ℑxyᶠᶠᵃ(i, j, 1, grid, fields.ζ)
-
     # replacement pressure?
     Pᵣ = ice_pressure(i, j, 1, grid, ip, rheology, fields)
-
-    ηᶜᶜᶜ = ζᶜᶜᶜ * e⁻²
-    ηᶠᶠᶜ = ζᶠᶠᶜ * e⁻²
+    η  = ζ * e⁻²
 
     # σ(uᵖ): the tangential stress depends only shear viscosity 
     # while the compressive stresses depend on the bulk viscosity and the ice strength
-    σ₁₁ᵖ⁺¹ = 2 * ηᶜᶜᶜ * ϵ̇₁₁ + ((ζᶜᶜᶜ - ηᶜᶜᶜ) * (ϵ̇₁₁ + ϵ̇₂₂) - Pᵣ / 2) 
-    σ₂₂ᵖ⁺¹ = 2 * ηᶜᶜᶜ * ϵ̇₂₂ + ((ζᶜᶜᶜ - ηᶜᶜᶜ) * (ϵ̇₁₁ + ϵ̇₂₂) - Pᵣ / 2)
-    σ₁₂ᵖ⁺¹ = 2 * ηᶠᶠᶜ * ϵ̇₁₂
-
-    mᵢᶜᶜᶜ = ice_mass(i, j, 1, grid, h, ℵ, ρᵢ) 
-    mᵢᶠᶠᶜ = ℑxyᶠᶠᵃ(i, j, 1, grid, ice_mass, h, ℵ, ρᵢ) 
+    σ₁₁ᵖ⁺¹ = 2 * η * ϵ̇₁₁ + ((ζ - η) * (ϵ̇₁₁ + ϵ̇₂₂) - Pᵣ / 2) 
+    σ₂₂ᵖ⁺¹ = 2 * η * ϵ̇₂₂ + ((ζ - η) * (ϵ̇₁₁ + ϵ̇₂₂) - Pᵣ / 2)
+    σ₁₂ᵖ⁺¹ = 2 * η * ϵ̇₁₂
 
+    mᵢ = ice_mass(i, j, kᴺ, grid, h, ℵ, ρᵢ) 
+
     # Update coefficients for substepping using dynamic substepping
     # with spatially varying coefficients as in Kimmritz et al (2016)
-    γ²ᶜᶜᶜ = ζᶜᶜᶜ * cα * Δt / mᵢᶜᶜᶜ / Azᶜᶜᶜ(i, j, 1, grid)
-    γ²ᶜᶜᶜ = ifelse(isnan(γ²ᶜᶜᶜ), α⁺^2, γ²ᶜᶜᶜ) # In case both ζᶜᶜᶜ and mᵢᶜᶜᶜ are zero
-    γᶜᶜᶜ  = clamp(sqrt(γ²ᶜᶜᶜ), α⁻, α⁺)
-
-    γ²ᶠᶠᶜ = ζᶠᶠᶜ * cα * Δt / mᵢᶠᶠᶜ / Azᶠᶠᶜ(i, j, 1, grid)
-    γ²ᶠᶠᶜ = ifelse(isnan(γ²ᶠᶠᶜ), α⁺^2, γ²ᶠᶠᶜ) # In case both ζᶠᶠᶜ and mᵢᶠᶠᶜ are zero
-    γᶠᶠᶜ  = clamp(sqrt(γ²ᶠᶠᶜ), α⁻, α⁺)
+    γ = ζ * cα * Δt / mᵢ * Az⁻¹ᶜᶜᶜ(i, j, kᴺ, grid)
+    α = clamp(sqrt(γ), α⁻, α⁺)
+    α = ifelse(isnan(α), α⁺, α)
 
     @inbounds begin
         # Compute the new stresses and store the value of the 
         # dynamic substepping coefficient α
-        σ₁₁★ = (σ₁₁ᵖ⁺¹ - σ₁₁[i, j, 1]) / γᶜᶜᶜ
-        σ₂₂★ = (σ₂₂ᵖ⁺¹ - σ₂₂[i, j, 1]) / γᶜᶜᶜ
-        σ₁₂★ = (σ₁₂ᵖ⁺¹ - σ₁₂[i, j, 1]) / γᶠᶠᶜ
-
-        σ₁₁[i, j, 1] += ifelse(mᵢᶜᶜᶜ > 0, σ₁₁★, zero(grid))
-        σ₂₂[i, j, 1] += ifelse(mᵢᶜᶜᶜ > 0, σ₂₂★, zero(grid))
-        σ₁₂[i, j, 1] += ifelse(mᵢᶠᶠᶜ > 0, σ₁₂★, zero(grid))
-          α[i, j, 1]  = γᶜᶜᶜ
+        σ₁₁★ = (σ₁₁ᵖ⁺¹ - σ₁₁[i, j, kᴺ]) / α
+        σ₂₂★ = (σ₂₂ᵖ⁺¹ - σ₂₂[i, j, kᴺ]) / α
+        σ₁₂★ = (σ₁₂ᵖ⁺¹ - σ₁₂[i, j, kᴺ]) / α
+
+        σ₁₁[i, j, kᴺ] += ifelse(mᵢ > 0, σ₁₁★, zero(grid))
+        σ₂₂[i, j, kᴺ] += ifelse(mᵢ > 0, σ₂₂★, zero(grid))
+        σ₁₂[i, j, kᴺ] += ifelse(mᵢ > 0, σ₁₂★, zero(grid))
+
+        fields.α[i, j, kᴺ] = α 
     end
 end
 
@@ -309,27 +287,26 @@ end
 #####
 
 # Here we extend all the functions that a rheology model needs to support:
-@inline ice_stress_ux(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = @inbounds fields.σ₁₁[i, j, k]
-@inline ice_stress_vx(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = @inbounds fields.σ₁₂[i, j, k]
-@inline ice_stress_uy(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = @inbounds fields.σ₁₂[i, j, k]
-@inline ice_stress_vy(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = @inbounds fields.σ₂₂[i, j, k]
+@inline ice_stress_ux(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = ℑyᵃᶠᵃ(i, j, k, grid, fields.σ₁₁)
+@inline ice_stress_uy(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = ℑxᶠᵃᵃ(i, j, k, grid, fields.σ₁₂)
+@inline ice_stress_vx(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = ℑyᵃᶠᵃ(i, j, k, grid, fields.σ₁₂)
+@inline ice_stress_vy(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = ℑxᶠᵃᵃ(i, j, k, grid, fields.σ₂₂)
 
 # To help convergence to the right velocities
-@inline compute_substep_Δtᶠᶜᶜ(i, j, grid, Δt, ::ElastoViscoPlasticRheology, substeps, fields) = Δt / ℑxᶠᵃᵃ(i, j, 1, grid, fields.α)
-@inline compute_substep_Δtᶜᶠᶜ(i, j, grid, Δt, ::ElastoViscoPlasticRheology, substeps, fields) = Δt / ℑyᵃᶠᵃ(i, j, 1, grid, fields.α)
+@inline compute_substep_Δtᶠᶠᶜ(i, j, grid, Δt, ::ElastoViscoPlasticRheology, substeps, fields) = @inbounds Δt / ℑxyᶠᶠᵃ(i, j, 1, grid, fields.α)
 
 #####
 ##### Numerical forcing to help convergence
 #####
 
 @inline function sum_of_forcing_u(i, j, k, grid, ::ElastoViscoPlasticRheology, u_forcing, fields, Δt) 
     user_forcing = u_forcing(i, j, k, grid, fields)
-    rheology_forcing = @inbounds (fields.uⁿ[i, j, k] - fields.u[i, j, k]) / Δt / ℑxᶠᵃᵃ(i, j, k, grid, fields.α)
+    rheology_forcing = @inbounds (fields.uⁿ[i, j, k] - fields.u[i, j, k]) / Δt / ℑxyᶠᶠᵃ(i, j, k, grid, fields.α)
     return user_forcing + rheology_forcing
 end
 
 @inline function sum_of_forcing_v(i, j, k, grid, ::ElastoViscoPlasticRheology, v_forcing, fields, Δt) 
     user_forcing = v_forcing(i, j, k, grid, fields)
-    rheology_forcing = @inbounds (fields.vⁿ[i, j, k] - fields.v[i, j, k]) / Δt / ℑyᵃᶠᵃ(i, j, k, grid, fields.α)
+    rheology_forcing = @inbounds (fields.vⁿ[i, j, k] - fields.v[i, j, k]) / Δt / ℑxyᶠᶠᵃ(i, j, k, grid, fields.α)
     return user_forcing + rheology_forcing
 end