diff --git a/experiments/OceanBoxGCM/refvals/ocean_gyre_refvals.jl b/experiments/OceanBoxGCM/refvals/ocean_gyre_refvals.jl
index 16ab65f501b..f08a5e2e7e8 100644
--- a/experiments/OceanBoxGCM/refvals/ocean_gyre_refvals.jl
+++ b/experiments/OceanBoxGCM/refvals/ocean_gyre_refvals.jl
@@ -105,6 +105,38 @@ varr = [
         3.78116102589123325624e-10,
         1.07073256826409470244e-05,
     ],
+    [
+        "s_aux",
+        "uᵈ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "uᵈ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "ΔGᵘ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "ΔGᵘ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
 ]
 parr = [
     ["Q", "u[1]", 12, 12, 12, 12],
@@ -115,6 +147,10 @@ parr = [
     ["s_aux", :w, 12, 12, 12, 12],
     ["s_aux", :pkin, 12, 12, 12, 12],
     ["s_aux", :wz0, 12, 12, 8, 12],
+    ["s_aux", "uᵈ[1]", 12, 12, 12, 12],
+    ["s_aux", "uᵈ[2]", 12, 12, 12, 12],
+    ["s_aux", "ΔGᵘ[1]", 12, 12, 12, 12],
+    ["s_aux", "ΔGᵘ[2]", 12, 12, 12, 12],
 ]
 # END SCPRINT
 # SC ====================================================================================
@@ -197,6 +233,38 @@ varr = [
         4.05688221636862481177e-10,
         1.10561449663397378925e-05,
     ],
+    [
+        "s_aux",
+        "uᵈ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "uᵈ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "ΔGᵘ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "ΔGᵘ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
 ]
 parr = [
     ["Q", "u[1]", 12, 12, 12, 12],
@@ -207,6 +275,10 @@ parr = [
     ["s_aux", :w, 12, 12, 12, 12],
     ["s_aux", :pkin, 12, 12, 12, 12],
     ["s_aux", :wz0, 12, 12, 8, 12],
+    ["s_aux", "uᵈ[1]", 12, 12, 12, 12],
+    ["s_aux", "uᵈ[2]", 12, 12, 12, 12],
+    ["s_aux", "ΔGᵘ[1]", 12, 12, 12, 12],
+    ["s_aux", "ΔGᵘ[2]", 12, 12, 12, 12],
 ]
 # END SCPRINT
 # SC ====================================================================================
diff --git a/src/Numerics/ODESolvers/SplitExplicitMethod.jl b/src/Numerics/ODESolvers/SplitExplicitMethod.jl
index 5f071ce6f80..fa369e3bb8c 100644
--- a/src/Numerics/ODESolvers/SplitExplicitMethod.jl
+++ b/src/Numerics/ODESolvers/SplitExplicitMethod.jl
@@ -36,8 +36,6 @@ mutable struct SplitExplicitSolver{SS, FS, RT, MSA} <: AbstractODESolver
     dt::RT
     "time"
     t::RT
-    "flag to couple the models"
-    coupled
     "storage for transfer tendency"
     dQ2fast::MSA
 
@@ -47,20 +45,20 @@ mutable struct SplitExplicitSolver{SS, FS, RT, MSA} <: AbstractODESolver
         Q = nothing;
         dt = getdt(slow_solver),
         t0 = slow_solver.t,
-        coupled = true,
     ) where {AT <: AbstractArray}
         SS = typeof(slow_solver)
         FS = typeof(fast_solver)
         RT = real(eltype(slow_solver.dQ))
 
         dQ2fast = similar(slow_solver.dQ)
+        dQ2fast .= -0
         MSA = typeof(dQ2fast)
+
         return new{SS, FS, RT, MSA}(
             slow_solver,
             fast_solver,
             RT(dt),
             RT(t0),
-            coupled,
             dQ2fast,
         )
     end
@@ -95,53 +93,28 @@ function dostep!(
 
         # Initialize fast model and tendency adjustment
         # before evalution of slow mode
-        if split.coupled
-            initialize_states!(
-                slow_bl,
-                fast_bl,
-                slow.rhs!,
-                fast.rhs!,
-                Qslow,
-                Qfast,
-            )
-
-            # Evaluate the slow mode
-            # --> save tendency for the fast
-            slow.rhs!(dQ2fast, Qslow, param, slow_stage_time, increment = false)
-
-            # vertically integrate slow tendency to advance fast equation
-            # and use vertical mean for slow model (negative source)
-            # ---> work with dQ2fast as input
-            tendency_from_slow_to_fast!(
-                slow_bl,
-                fast_bl,
-                slow.rhs!,
-                fast.rhs!,
-                Qslow,
-                Qfast,
-                dQ2fast,
-            )
-        end
+        initialize_states!(slow_bl, fast_bl, slow.rhs!, fast.rhs!, Qslow, Qfast)
+
+        # Evaluate the slow mode
+        # --> save tendency for the fast
+        slow.rhs!(dQ2fast, Qslow, param, slow_stage_time, increment = false)
+
+        # vertically integrate slow tendency to advance fast equation
+        # and use vertical mean for slow model (negative source)
+        # ---> work with dQ2fast as input
+        tendency_from_slow_to_fast!(
+            slow_bl,
+            fast_bl,
+            slow.rhs!,
+            fast.rhs!,
+            Qslow,
+            Qfast,
+            dQ2fast,
+        )
 
         # Compute (and RK update) slow tendency
         slow.rhs!(dQslow, Qslow, param, slow_stage_time, increment = true)
 
-        # Update (RK-stage) slow state
-        event = Event(array_device(Qslow))
-        event = update!(array_device(Qslow), groupsize)(
-            realview(dQslow),
-            realview(Qslow),
-            slow.RKA[slow_s % length(slow.RKA) + 1],
-            slow.RKB[slow_s],
-            slow_dt,
-            nothing,
-            nothing,
-            nothing;
-            ndrange = length(realview(Qslow)),
-            dependencies = (event,),
-        )
-        wait(array_device(Qslow), event)
-
         # Fractional time for slow stage
         if slow_s == length(slow.RKA)
             γ = 1 - slow.RKC[slow_s]
@@ -159,30 +132,42 @@ function dostep!(
         for substep in 1:nsubsteps
             fast_time = slow_stage_time + (substep - 1) * fast_dt
             dostep!(Qfast, fast, param, fast_time)
-            if split.coupled
-                cummulate_fast_solution!(
-                    slow_bl,
-                    fast_bl,
-                    fast.rhs!,
-                    Qfast,
-                    fast_time,
-                    fast_dt,
-                    substep,
-                )
-            end
-        end
-
-        # reconcile slow equation using fast equation
-        if split.coupled
-            reconcile_from_fast_to_slow!(
+            cummulate_fast_solution!(
                 slow_bl,
                 fast_bl,
-                slow.rhs!,
                 fast.rhs!,
-                Qslow,
                 Qfast,
+                fast_time,
+                fast_dt,
+                substep,
             )
         end
+
+        # Update (RK-stage) slow state
+        event = Event(array_device(Qslow))
+        event = update!(array_device(Qslow), groupsize)(
+            realview(dQslow),
+            realview(Qslow),
+            slow.RKA[slow_s % length(slow.RKA) + 1],
+            slow.RKB[slow_s],
+            slow_dt,
+            nothing,
+            nothing,
+            nothing;
+            ndrange = length(realview(Qslow)),
+            dependencies = (event,),
+        )
+        wait(array_device(Qslow), event)
+
+        # reconcile slow equation using fast equation
+        reconcile_from_fast_to_slow!(
+            slow_bl,
+            fast_bl,
+            slow.rhs!,
+            fast.rhs!,
+            Qslow,
+            Qfast,
+        )
     end
     updatedt!(fast, fast_dt_in)
 
diff --git a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
index 1b95b393d7d..53defec666d 100644
--- a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
+++ b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
@@ -6,6 +6,7 @@ using StaticArrays
 using LinearAlgebra: dot, Diagonal
 using CLIMAParameters.Planet: grav
 
+using ..Ocean
 using ...VariableTemplates
 using ...MPIStateArrays
 using ...Mesh.Filters: apply!
@@ -13,11 +14,11 @@ using ...Mesh.Grids: VerticalDirection
 using ...Mesh.Geometry
 using ...DGMethods
 using ...DGMethods.NumericalFluxes
-import ...DGMethods.NumericalFluxes: update_penalty!
 using ...DGMethods.NumericalFluxes: RusanovNumericalFlux
 using ...BalanceLaws
-using ...BalanceLaws: number_states
 
+import ..Ocean: coriolis_parameter
+import ...DGMethods.NumericalFluxes: update_penalty!
 import ...BalanceLaws:
     vars_state,
     init_state_prognostic!,
@@ -67,9 +68,10 @@ fₒ = first coriolis parameter (constant term)
     HydrostaticBoussinesqModel(problem)
 
 """
-struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw
+struct HydrostaticBoussinesqModel{C, PS, P, T} <: BalanceLaw
     param_set::PS
     problem::P
+    coupling::C
     ρₒ::T
     cʰ::T
     cᶻ::T
@@ -83,7 +85,8 @@ struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw
     β::T
     function HydrostaticBoussinesqModel{FT}(
         param_set::PS,
-        problem;
+        problem::P;
+        coupling::C = Uncoupled(),
         ρₒ = FT(1000),  # kg / m^3
         cʰ = FT(0),     # m/s
         cᶻ = FT(0),     # m/s
@@ -95,10 +98,11 @@ struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw
         κᶜ = FT(1e-1),  # m^2 / s # diffusivity for convective adjustment
         fₒ = FT(1e-4),  # Hz
         β = FT(1e-11), # Hz / m
-    ) where {FT <: AbstractFloat, PS}
-        return new{PS, typeof(problem), FT}(
+    ) where {FT <: AbstractFloat, PS, P, C}
+        return new{C, PS, P, FT}(
             param_set,
             problem,
+            coupling,
             ρₒ,
             cʰ,
             cᶻ,
@@ -164,6 +168,8 @@ function vars_state(m::HBModel, ::Auxiliary, T)
         w::T     # ∫(-∇⋅u)
         pkin::T  # ∫(-αᵀθ)
         wz0::T   # w at z=0
+        uᵈ::SVector{2, T}    # velocity deviation from vertical mean
+        ΔGᵘ::SVector{2, T}   # vertically averaged tendency
     end
 end
 
@@ -187,6 +193,7 @@ these are just copies in our model
 function vars_state(m::HBModel, ::Gradient, T)
     @vars begin
         ∇u::SVector{2, T}
+        ∇uᵈ::SVector{2, T}
         ∇θ::T
     end
 end
@@ -205,9 +212,23 @@ this computation is done pointwise at each nodal point
 - `t`: time, not used
 """
 @inline function compute_gradient_argument!(m::HBModel, G::Vars, Q::Vars, A, t)
-    G.∇u = Q.u
     G.∇θ = Q.θ
 
+    velocity_gradient_argument!(m, m.coupling, G, Q, A, t)
+
+    return nothing
+end
+
+@inline function velocity_gradient_argument!(
+    m::HBModel,
+    ::Uncoupled,
+    G,
+    Q,
+    A,
+    t,
+)
+    G.∇u = Q.u
+
     return nothing
 end
 
@@ -250,8 +271,7 @@ this computation is done pointwise at each nodal point
     # store ∇ʰu for continuity equation (convert gradient to divergence)
     D.∇ʰu = G.∇u[1, 1] + G.∇u[2, 2]
 
-    ν = viscosity_tensor(m)
-    D.ν∇u = -ν * G.∇u
+    velocity_gradient_flux!(m, m.coupling, D, G, Q, A, t)
 
     κ = diffusivity_tensor(m, G.∇θ[3])
     D.κ∇θ = -κ * G.∇θ
@@ -259,6 +279,13 @@ this computation is done pointwise at each nodal point
     return nothing
 end
 
+@inline function velocity_gradient_flux!(m::HBModel, ::Uncoupled, D, G, Q, A, t)
+    ν = viscosity_tensor(m)
+    D.ν∇u = -ν * G.∇u
+
+    return nothing
+end
+
 """
     viscosity_tensor(::HBModel)
 
@@ -420,38 +447,44 @@ t -> time, not used
     t::Real,
     direction,
 )
-    FT = eltype(Q)
-    _grav::FT = grav(m.param_set)
     @inbounds begin
-        u = Q.u # Horizontal components of velocity
-        η = Q.η
-        θ = Q.θ
-        w = A.w   # vertical velocity
-        pkin = A.pkin
+        # ∇h • (g η)
+        hydrostatic_pressure!(m, m.coupling, F, Q, A, t)
 
-        v = @SVector [u[1], u[2], w]
+        # ∇h • (- ∫(αᵀ θ))
+        pkin = A.pkin
         Iʰ = @SMatrix [
             1 -0
             -0 1
             -0 -0
         ]
-
-        # ∇h • (g η)
-        F.u += _grav * η * Iʰ
-
-        # ∇h • (- ∫(αᵀ θ))
-        F.u += _grav * pkin * Iʰ
+        F.u += grav(m.param_set) * pkin * Iʰ
 
         # ∇h • (v ⊗ u)
         # F.u += v * u'
 
         # ∇ • (u θ)
+        θ = Q.θ
+        v = @SVector [Q.u[1], Q.u[2], A.w]
         F.θ += v * θ
     end
 
     return nothing
 end
 
+@inline function hydrostatic_pressure!(m::HBModel, ::Uncoupled, F, Q, A, t)
+    η = Q.η
+    Iʰ = @SMatrix [
+        1 -0
+        -0 1
+        -0 -0
+    ]
+
+    F.u += grav(m.param_set) * η * Iʰ
+
+    return nothing
+end
+
 """
     flux_second_order!(::HBModel)
 
@@ -510,22 +543,25 @@ end
     t::Real,
     direction,
 )
-    @inbounds begin
-        u, v = Q.u # Horizontal components of velocity
-        wz0 = A.wz0
+    wz0 = A.wz0
+    S.η += wz0
 
-        # f × u
-        f = coriolis_force(m, A.y)
-        S.u -= @SVector [-f * v, f * u]
+    coriolis_force!(m, m.coupling, S, Q, A, t)
 
-        S.η += wz0
-    end
+    return nothing
+end
+
+@inline function coriolis_force!(m::HBModel, ::Uncoupled, S, Q, A, t)
+    # f × u
+    f = coriolis_parameter(m, A.y)
+    u, v = Q.u # Horizontal components of velocity
+    S.u -= @SVector [-f * v, f * u]
 
     return nothing
 end
 
 """
-    coriolis_force(::HBModel)
+    coriolis_parameter(::HBModel)
 
 northern hemisphere coriolis
 
@@ -533,7 +569,7 @@ northern hemisphere coriolis
 - `m`: model object to dispatch on and get coriolis parameters
 - `y`: y-coordinate in the box
 """
-@inline coriolis_force(m::HBModel, y) = m.fₒ + m.β * y
+@inline coriolis_parameter(m::HBModel, y) = m.fₒ + m.β * y
 
 """
     wavespeed(::HBModel)
@@ -593,9 +629,13 @@ function update_auxiliary_state!(
         apply!(Q, (:θ,), dg.grid, exp_filter, direction = VerticalDirection())
     end
 
+    compute_flow_deviation!(dg, m, m.coupling, Q, t)
+
     return true
 end
 
+@inline compute_flow_deviation!(dg, ::HBModel, ::Uncoupled, _...) = nothing
+
 """
     update_auxiliary_state_gradient!(::HBModel)
 
@@ -634,13 +674,13 @@ function update_auxiliary_state_gradient!(
     # We are unable to use vars (ie A.w) for this because this operation will
     # return a SubArray, and adapt (used for broadcasting along reshaped arrays)
     # has a limited recursion depth for the types allowed.
-    number_auxiliary = number_states(m, Auxiliary(), FT)
+    number_aux = number_states(m, Auxiliary(), FT)
     index_w = varsindex(vars_state(m, Auxiliary(), FT), :w)
     index_wz0 = varsindex(vars_state(m, Auxiliary(), FT), :wz0)
     Nq, Nqk, _, _, nelemv, nelemh, nhorzrealelem, _ = basic_grid_info(dg)
 
     # project w(z=0) down the stack
-    data = reshape(A.data, Nq^2, Nqk, number_auxiliary, nelemv, nelemh)
+    data = reshape(A.data, Nq^2, Nqk, number_aux, nelemv, nelemh)
     flat_wz0 = @view data[:, end:end, index_w, end:end, 1:nhorzrealelem]
     boxy_wz0 = @view data[:, :, index_wz0, :, 1:nhorzrealelem]
     boxy_wz0 .= flat_wz0
diff --git a/src/Ocean/HydrostaticBoussinesq/SimpleBoxProblem.jl b/src/Ocean/HydrostaticBoussinesq/SimpleBoxProblem.jl
index 4d4506fd166..705dc08f67f 100644
--- a/src/Ocean/HydrostaticBoussinesq/SimpleBoxProblem.jl
+++ b/src/Ocean/HydrostaticBoussinesq/SimpleBoxProblem.jl
@@ -52,6 +52,9 @@ function ocean_init_aux!(m::HBModel, p::AbstractSimpleBoxProblem, A, geom)
     A.pkin = 0
     A.wz0 = 0
 
+    A.uᵈ = @SVector [-0, -0]
+    A.ΔGᵘ = @SVector [-0, -0]
+
     return nothing
 end
 
diff --git a/src/Ocean/Ocean.jl b/src/Ocean/Ocean.jl
index 9dba464b5f5..5136bb84cc8 100644
--- a/src/Ocean/Ocean.jl
+++ b/src/Ocean/Ocean.jl
@@ -1,5 +1,13 @@
 module Ocean
 
+export AbstractOceanCoupling, Uncoupled, Coupled
+
+abstract type AbstractOceanCoupling end
+struct Uncoupled <: AbstractOceanCoupling end
+struct Coupled <: AbstractOceanCoupling end
+
+function coriolis_parameter end
+
 include("HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl")
 include("ShallowWater/ShallowWaterModel.jl")
 include("SplitExplicit/SplitExplicitModel.jl")
diff --git a/src/Ocean/ShallowWater/ShallowWaterModel.jl b/src/Ocean/ShallowWater/ShallowWaterModel.jl
index 074df5ae4cb..f5d41070c46 100644
--- a/src/Ocean/ShallowWater/ShallowWaterModel.jl
+++ b/src/Ocean/ShallowWater/ShallowWaterModel.jl
@@ -6,6 +6,7 @@ using StaticArrays
 using LinearAlgebra: dot, Diagonal
 using CLIMAParameters.Planet: grav
 
+using ..Ocean
 using ...VariableTemplates
 using ...Mesh.Geometry
 using ...DGMethods
@@ -42,9 +43,10 @@ end
 abstract type AdvectionTerm end
 struct NonLinearAdvection <: AdvectionTerm end
 
-struct ShallowWaterModel{PS, P, T, A, S} <: BalanceLaw
+struct ShallowWaterModel{C, PS, P, T, A, S} <: BalanceLaw
     param_set::PS
     problem::P
+    coupling::C
     turbulence::T
     advection::A
     c::S
@@ -63,6 +65,8 @@ function vars_state(m::SWModel, ::Auxiliary, T)
     @vars begin
         f::T
         τ::SVector{2, T}  # value includes τₒ, g, and ρ
+        Gᵁ::SVector{2, T} # integral of baroclinic tendency
+        Δu::SVector{2, T} # reconciliation Δu = 1/H * (Ū - ∫u)
     end
 end
 
@@ -154,12 +158,13 @@ function compute_gradient_flux!(
     α::Vars,
     t::Real,
 )
-    compute_gradient_flux!(m.turbulence, σ, δ, q, α, t)
+    compute_gradient_flux!(m, m.turbulence, σ, δ, q, α, t)
 end
 
-compute_gradient_flux!(::LinearDrag, _...) = nothing
+compute_gradient_flux!(::SWModel, ::LinearDrag, _...) = nothing
 
 @inline function compute_gradient_flux!(
+    ::SWModel,
     T::ConstantViscosity,
     σ::Vars,
     δ::Grad,
@@ -184,12 +189,13 @@ function flux_second_order!(
     α::Vars,
     t::Real,
 )
-    flux_second_order!(m.turbulence, G, q, σ, α, t)
+    flux_second_order!(m, m.turbulence, G, q, σ, α, t)
 end
 
-flux_second_order!(::LinearDrag, _...) = nothing
+flux_second_order!(::SWModel, ::LinearDrag, _...) = nothing
 
 @inline function flux_second_order!(
+    ::SWModel,
     ::ConstantViscosity,
     G::Grad,
     q::Vars,
@@ -208,22 +214,28 @@ end
     m::SWModel{P},
     S::Vars,
     q::Vars,
-    diffusive::Vars,
+    d::Vars,
     α::Vars,
     t::Real,
     direction,
 ) where {P}
-    S.U += α.τ
-
+    # f × u
     f = α.f
     U, V = q.U
     S.U -= @SVector [-f * V, f * U]
 
+    forcing_term!(m, m.coupling, S, q, α, t)
     linear_drag!(m.turbulence, S, q, α, t)
 
     return nothing
 end
 
+@inline function forcing_term!(::SWModel, ::Uncoupled, S, Q, A, t)
+    S.U += A.τ
+
+    return nothing
+end
+
 linear_drag!(::ConstantViscosity, _...) = nothing
 
 @inline function linear_drag!(T::LinearDrag, S::Vars, q::Vars, α::Vars, t::Real)
@@ -237,7 +249,7 @@ end
 
 function shallow_init_aux! end
 function init_state_auxiliary!(m::SWModel, aux::Vars, geom::LocalGeometry)
-    shallow_init_aux!(m.problem, aux, geom)
+    shallow_init_aux!(m, m.problem, aux, geom)
 end
 
 function shallow_init_state! end
@@ -245,6 +257,7 @@ function init_state_prognostic!(m::SWModel, state::Vars, aux::Vars, coords, t)
     shallow_init_state!(m, m.problem, state, aux, coords, t)
 end
 
+function shallow_boundary_state! end
 function boundary_state!(
     nf,
     m::SWModel,
@@ -260,8 +273,25 @@ function boundary_state!(
     shallow_boundary_state!(nf, m, m.turbulence, q⁺, a⁺, n⁻, q⁻, a⁻, t)
 end
 
+function boundary_state!(
+    nf,
+    m::SWModel,
+    q⁺::Vars,
+    σ⁺::Vars,
+    α⁺::Vars,
+    n⁻,
+    q⁻::Vars,
+    σ⁻::Vars,
+    α⁻::Vars,
+    bctype,
+    t,
+    _...,
+)
+    shallow_boundary_state!(nf, m, m.turbulence, q⁺, σ⁺, α⁺, n⁻, q⁻, σ⁻, α⁻, t)
+end
+
 @inline function shallow_boundary_state!(
-    ::RusanovNumericalFlux,
+    ::NumericalFluxFirstOrder,
     m::SWModel,
     ::LinearDrag,
     q⁺,
@@ -281,38 +311,21 @@ end
 end
 
 shallow_boundary_state!(
-    ::CentralNumericalFluxGradient,
+    ::NumericalFluxGradient,
     m::SWModel,
     ::LinearDrag,
     _...,
 ) = nothing
 
 shallow_boundary_state!(
-    ::CentralNumericalFluxSecondOrder,
+    ::NumericalFluxSecondOrder,
     m::SWModel,
     ::LinearDrag,
     _...,
 ) = nothing
 
-function boundary_state!(
-    nf,
-    m::SWModel,
-    q⁺::Vars,
-    σ⁺::Vars,
-    α⁺::Vars,
-    n⁻,
-    q⁻::Vars,
-    σ⁻::Vars,
-    α⁻::Vars,
-    bctype,
-    t,
-    _...,
-)
-    shallow_boundary_state!(nf, m, m.turbulence, q⁺, σ⁺, α⁺, n⁻, q⁻, σ⁻, α⁻, t)
-end
-
 @inline function shallow_boundary_state!(
-    ::RusanovNumericalFlux,
+    ::NumericalFluxFirstOrder,
     m::SWModel,
     ::ConstantViscosity,
     q⁺,
@@ -329,7 +342,7 @@ end
 end
 
 @inline function shallow_boundary_state!(
-    ::CentralNumericalFluxGradient,
+    ::NumericalFluxGradient,
     m::SWModel,
     ::ConstantViscosity,
     q⁺,
@@ -346,7 +359,7 @@ end
 end
 
 @inline function shallow_boundary_state!(
-    ::CentralNumericalFluxSecondOrder,
+    ::NumericalFluxSecondOrder,
     m::SWModel,
     ::ConstantViscosity,
     q⁺,
diff --git a/src/Ocean/SplitExplicit/Communication.jl b/src/Ocean/SplitExplicit/Communication.jl
new file mode 100644
index 00000000000..f7166d7ac96
--- /dev/null
+++ b/src/Ocean/SplitExplicit/Communication.jl
@@ -0,0 +1,144 @@
+@inline function initialize_states!(
+    baroclinic::HBModel{C},
+    barotropic::SWModel{C},
+    model_bc,
+    model_bt,
+    state_bc,
+    state_bt,
+) where {C <: Coupled}
+    model_bc.state_auxiliary.ΔGᵘ .= -0
+
+    return nothing
+end
+
+@inline function tendency_from_slow_to_fast!(
+    baroclinic::HBModel{C},
+    barotropic::SWModel{C},
+    model_bc,
+    model_bt,
+    state_bc,
+    state_bt,
+    forcing_tendency,
+) where {C <: Coupled}
+    FT = eltype(state_bc)
+    Nq, Nqk, _, _, nelemv, nelemh, nrealelemh, _ = basic_grid_info(model_bc)
+
+    #### integrate the tendency
+    model_int = model_bc.modeldata.integral_model
+    integral = model_int.balance_law
+    update_auxiliary_state!(model_int, integral, forcing_tendency, 0)
+
+    ### properly shape MPIStateArrays
+    num_aux_int = number_states(integral, Auxiliary(), FT)
+    data_int = model_int.state_auxiliary.data
+    data_int = reshape(data_int, Nq^2, Nqk, num_aux_int, nelemv, nelemh)
+
+    num_aux_bt = number_states(barotropic, Auxiliary(), FT)
+    data_bt = model_bt.state_auxiliary.data
+    data_bt = reshape(data_bt, Nq^2, num_aux_bt, nelemh)
+
+    num_aux_bc = number_states(baroclinic, Auxiliary(), FT)
+    data_bc = model_bc.state_auxiliary.data
+    data_bc = reshape(data_bc, Nq^2, Nqk, num_aux_bc, nelemv, nelemh)
+
+    ### get vars indices
+    index_∫du = varsindex(vars_state(integral, Auxiliary(), FT), :(∫x))
+    index_Gᵁ = varsindex(vars_state(barotropic, Auxiliary(), FT), :Gᵁ)
+    index_ΔGᵘ = varsindex(vars_state(baroclinic, Auxiliary(), FT), :ΔGᵘ)
+
+    ### get top value (=integral over full depth) of ∫du
+    ∫du = @view data_int[:, end, index_∫du, end, 1:nrealelemh]
+
+    ### copy into Gᵁ of barotropic model
+    Gᵁ = @view data_bt[:, index_Gᵁ, 1:nrealelemh]
+    Gᵁ .= ∫du
+
+    ### get top value (=integral over full depth) of ∫du
+    ∫du = @view data_int[:, end:end, index_∫du, end:end, 1:nrealelemh]
+
+    ### save vertically averaged tendency to remove from 3D tendency
+    ### need to reshape for the broadcast
+    ΔGᵘ = @view data_bc[:, :, index_ΔGᵘ, :, 1:nrealelemh]
+    ΔGᵘ .-= ∫du / baroclinic.problem.H
+
+    return nothing
+end
+
+@inline function cummulate_fast_solution!(
+    baroclinic::HBModel{C},
+    barotropic::SWModel{C},
+    model_bt,
+    state_bt,
+    fast_time,
+    fast_dt,
+    substep,
+) where {C <: Coupled}
+    return nothing
+end
+
+@inline function reconcile_from_fast_to_slow!(
+    baroclinic::HBModel{C},
+    barotropic::SWModel{C},
+    model_bc,
+    model_bt,
+    state_bc,
+    state_bt,
+) where {C <: Coupled}
+    FT = eltype(state_bc)
+    Nq, Nqk, _, _, nelemv, nelemh, nrealelemh, _ = basic_grid_info(model_bc)
+
+    ### integrate the horizontal velocity
+    model_int = model_bc.modeldata.integral_model
+    integral = model_int.balance_law
+    update_auxiliary_state!(model_int, integral, state_bc, 0)
+
+    ### properly shape MPIStateArrays
+    num_aux_int = number_states(integral, Auxiliary(), FT)
+    data_int = model_int.state_auxiliary.data
+    data_int = reshape(data_int, Nq^2, Nqk, num_aux_int, nelemv, nelemh)
+
+    num_aux_bt = number_states(barotropic, Auxiliary(), FT)
+    data_bt_aux = model_bt.state_auxiliary.data
+    data_bt_aux = reshape(data_bt_aux, Nq^2, num_aux_bt, nelemh)
+
+    num_state_bt = number_states(barotropic, Prognostic(), FT)
+    data_bt_state = state_bt.data
+    data_bt_state = reshape(data_bt_state, Nq^2, num_state_bt, nelemh)
+
+    num_state_bc = number_states(baroclinic, Prognostic(), FT)
+    data_bc_state = state_bc.data
+    data_bc_state =
+        reshape(data_bc_state, Nq^2, Nqk, num_state_bc, nelemv, nelemh)
+
+    ### get vars indices
+    index_∫u = varsindex(vars_state(integral, Auxiliary(), FT), :(∫x))
+    index_Δu = varsindex(vars_state(barotropic, Auxiliary(), FT), :Δu)
+    index_U = varsindex(vars_state(barotropic, Prognostic(), FT), :U)
+    index_u = varsindex(vars_state(baroclinic, Prognostic(), FT), :u)
+    index_η_3D = varsindex(vars_state(baroclinic, Prognostic(), FT), :η)
+    index_η_2D = varsindex(vars_state(barotropic, Prognostic(), FT), :η)
+
+    ### get top value (=integral over full depth)
+    ∫u = @view data_int[:, end, index_∫u, end, 1:nrealelemh]
+
+    ### Δu is a place holder for 1/H * (Ū - ∫u)
+    Δu = @view data_bt_aux[:, index_Δu, 1:nrealelemh]
+    U = @view data_bt_state[:, index_U, 1:nrealelemh]
+    Δu .= 1 / baroclinic.problem.H * (U - ∫u)
+
+    ### copy the 2D contribution down the 3D solution
+    ### need to reshape for the broadcast
+    data_bt_aux = reshape(data_bt_aux, Nq^2, 1, num_aux_bt, 1, nelemh)
+    Δu = @view data_bt_aux[:, :, index_Δu, :, 1:nrealelemh]
+    u = @view data_bc_state[:, :, index_u, :, 1:nrealelemh]
+    u .+= Δu
+
+    ### copy η from barotropic mode to baroclinic mode
+    ### need to reshape for the broadcast
+    data_bt_state = reshape(data_bt_state, Nq^2, 1, num_state_bt, 1, nelemh)
+    η_2D = @view data_bt_state[:, :, index_η_2D, :, 1:nrealelemh]
+    η_3D = @view data_bc_state[:, :, index_η_3D, :, 1:nrealelemh]
+    η_3D .= η_2D
+
+    return nothing
+end
diff --git a/src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl b/src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl
new file mode 100644
index 00000000000..a6e4420b586
--- /dev/null
+++ b/src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl
@@ -0,0 +1,82 @@
+using ..Ocean: coriolis_parameter
+using ..HydrostaticBoussinesq
+using ...DGMethods
+
+import ...BalanceLaws
+import ..HydrostaticBoussinesq:
+    viscosity_tensor,
+    coriolis_force!,
+    velocity_gradient_argument!,
+    velocity_gradient_flux!,
+    hydrostatic_pressure!,
+    compute_flow_deviation!
+
+@inline function velocity_gradient_argument!(m::HBModel, ::Coupled, G, Q, A, t)
+    G.∇u = Q.u
+    G.∇uᵈ = A.uᵈ
+
+    return nothing
+end
+
+@inline function velocity_gradient_flux!(m::HBModel, ::Coupled, D, G, Q, A, t)
+    ν = viscosity_tensor(m)
+    ∇u = @SMatrix [
+        G.∇uᵈ[1, 1] G.∇uᵈ[1, 2]
+        G.∇uᵈ[2, 1] G.∇uᵈ[2, 2]
+        G.∇u[3, 1] G.∇u[3, 2]
+    ]
+    D.ν∇u = -ν * ∇u
+
+    return nothing
+end
+
+@inline hydrostatic_pressure!(::HBModel, ::Coupled, _...) = nothing
+
+@inline function coriolis_force!(m::HBModel, ::Coupled, S, Q, A, t)
+    # f × u
+    f = coriolis_parameter(m, A.y)
+    uᵈ, vᵈ = A.uᵈ # Horizontal components of velocity
+    S.u -= @SVector [-f * vᵈ, f * uᵈ]
+
+    return nothing
+end
+
+# Compute Horizontal Flow deviation from vertical mean
+@inline function compute_flow_deviation!(dg, m::HBModel, ::Coupled, Q, t)
+    FT = eltype(Q)
+    Nq, Nqk, _, _, nelemv, nelemh, nrealelemh, _ = basic_grid_info(dg)
+
+    #### integrate the tendency
+    model_int = dg.modeldata.integral_model
+    integral = model_int.balance_law
+    update_auxiliary_state!(model_int, integral, Q, 0)
+
+    ### properly shape MPIStateArrays
+    num_int = number_states(integral, Auxiliary(), FT)
+    data_int = model_int.state_auxiliary.data
+    data_int = reshape(data_int, Nq^2, Nqk, num_int, nelemv, nelemh)
+
+    num_aux = number_states(m, Auxiliary(), FT)
+    data_aux = dg.state_auxiliary.data
+    data_aux = reshape(data_aux, Nq^2, Nqk, num_aux, nelemv, nelemh)
+
+    num_state = number_states(m, Prognostic(), FT)
+    data_state = reshape(Q.data, Nq^2, Nqk, num_state, nelemv, nelemh)
+
+    ### get vars indices
+    index_∫u = varsindex(vars_state(integral, Auxiliary(), FT), :(∫x))
+    index_uᵈ = varsindex(vars_state(m, Auxiliary(), FT), :uᵈ)
+    index_u = varsindex(vars_state(m, Prognostic(), FT), :u)
+
+    ### get top value (=integral over full depth)
+    ∫u = @view data_int[:, end:end, index_∫u, end:end, 1:nrealelemh]
+    uᵈ = @view data_aux[:, :, index_uᵈ, :, 1:nrealelemh]
+    u = @view data_state[:, :, index_u, :, 1:nrealelemh]
+
+    ## make a copy of horizontal velocity
+    ## and remove vertical mean velocity
+    uᵈ .= u
+    uᵈ .-= ∫u / m.problem.H
+
+    return nothing
+end
diff --git a/src/Ocean/SplitExplicit/ShallowWaterCoupling.jl b/src/Ocean/SplitExplicit/ShallowWaterCoupling.jl
new file mode 100644
index 00000000000..9c19e161bf7
--- /dev/null
+++ b/src/Ocean/SplitExplicit/ShallowWaterCoupling.jl
@@ -0,0 +1,7 @@
+import ..ShallowWater: forcing_term!
+
+@inline function forcing_term!(::SWModel, ::Coupled, S, Q, A, t)
+    S.U += A.Gᵁ
+
+    return nothing
+end
diff --git a/src/Ocean/SplitExplicit/SplitExplicitModel.jl b/src/Ocean/SplitExplicit/SplitExplicitModel.jl
index dc285ffcafe..46c73fb0f83 100644
--- a/src/Ocean/SplitExplicit/SplitExplicitModel.jl
+++ b/src/Ocean/SplitExplicit/SplitExplicitModel.jl
@@ -2,6 +2,7 @@ module SplitExplicit
 
 using StaticArrays
 
+using ..Ocean
 using ..HydrostaticBoussinesq
 using ..ShallowWater
 
@@ -11,36 +12,47 @@ using ...Mesh.Geometry
 using ...DGMethods
 using ...BalanceLaws
 
-include("VerticalIntegralModel.jl")
-
 import ...BalanceLaws:
     initialize_states!,
     tendency_from_slow_to_fast!,
     cummulate_fast_solution!,
     reconcile_from_fast_to_slow!
 
-@inline initialize_states!(
-    slow::HydrostaticBoussinesqModel,
-    fast::ShallowWaterModel,
-    _...,
-) = nothing
+HBModel = HydrostaticBoussinesqModel
+SWModel = ShallowWaterModel
 
-@inline tendency_from_slow_to_fast!(
-    slow::HydrostaticBoussinesqModel,
-    fast::ShallowWaterModel,
+function initialize_states!(
+    ::HBModel{C},
+    ::SWModel{C},
     _...,
-) = nothing
-
-@inline cummulate_fast_solution!(
-    slow::HydrostaticBoussinesqModel,
-    fast::ShallowWaterModel,
+) where {C <: Uncoupled}
+    return nothing
+end
+function tendency_from_slow_to_fast!(
+    ::HBModel{C},
+    ::SWModel{C},
     _...,
-) = nothing
-
-@inline reconcile_from_fast_to_slow!(
-    slow::HydrostaticBoussinesqModel,
-    fast::ShallowWaterModel,
+) where {C <: Uncoupled}
+    return nothing
+end
+function cummulate_fast_solution!(
+    ::HBModel{C},
+    ::SWModel{C},
+    _...,
+) where {C <: Uncoupled}
+    return nothing
+end
+function reconcile_from_fast_to_slow!(
+    ::HBModel{C},
+    ::SWModel{C},
     _...,
-) = nothing
+) where {C <: Uncoupled}
+    return nothing
+end
+
+include("VerticalIntegralModel.jl")
+include("Communication.jl")
+include("ShallowWaterCoupling.jl")
+include("HydrostaticBoussinesqCoupling.jl")
 
 end
diff --git a/test/Ocean/HydrostaticBoussinesq/test_divergence_free.jl b/test/Ocean/HydrostaticBoussinesq/test_divergence_free.jl
index a9e26ad6454..3240e1cfd97 100644
--- a/test/Ocean/HydrostaticBoussinesq/test_divergence_free.jl
+++ b/test/Ocean/HydrostaticBoussinesq/test_divergence_free.jl
@@ -102,14 +102,12 @@ function test_divergence_free(; imex::Bool = false, BC = nothing)
 end
 
 @testset "$(@__FILE__)" begin
-    boundary_conditions = [(
+    BC = (
         OceanBC(Impenetrable(NoSlip()), Insulating()),
         OceanBC(Impenetrable(NoSlip()), Insulating()),
         OceanBC(Penetrable(KinematicStress()), Insulating()),
-    )]
+    )
 
-    for BC in boundary_conditions
-        test_divergence_free(imex = false, BC = BC)
-        test_divergence_free(imex = true, BC = BC)
-    end
+    test_divergence_free(imex = false, BC = BC)
+    test_divergence_free(imex = true, BC = BC)
 end
diff --git a/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl b/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl
index a2a051dae7a..05734beecb1 100644
--- a/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl
+++ b/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl
@@ -39,7 +39,6 @@ function test_ocean_gyre(;
     imex::Bool = false,
     BC = nothing,
     Δt = 60,
-    nt = 0,
     refDat = (),
 )
     FT = Float64
@@ -109,11 +108,7 @@ function test_ocean_gyre(;
     regenRefVals = false
     if regenRefVals
         ## Print state statistics in format for use as reference values
-        println(
-            "# SC ========== Test number ",
-            nt,
-            " reference values and precision match template. =======",
-        )
+        println("# SC ========== Test number reference values and precision match template. =======",)
         println("# SC ========== $(@__FILE__) test reference values ======================================")
         ClimateMachine.StateCheck.scprintref(cb)
         println("# SC ====================================================================================")
@@ -134,31 +129,13 @@ end
 
     include("../refvals/test_ocean_gyre_refvals.jl")
 
-    nt = 1
-    boundary_conditions = [
-        (
-            OceanBC(Impenetrable(NoSlip()), Insulating()),
-            OceanBC(Impenetrable(NoSlip()), Insulating()),
-            OceanBC(Penetrable(KinematicStress()), TemperatureFlux()),
-        ),
-    ]
-
-    for BC in boundary_conditions
-        test_ocean_gyre(
-            imex = false,
-            BC = BC,
-            Δt = 600,
-            nt = nt,
-            refDat = (refVals[nt], refPrecs[nt]),
-        )
-        nt = nt + 1
-        test_ocean_gyre(
-            imex = true,
-            BC = BC,
-            Δt = 150,
-            nt = nt,
-            refDat = (refVals[nt], refPrecs[nt]),
-        )
-        nt = nt + 1
-    end
+    BC = (
+        OceanBC(Impenetrable(NoSlip()), Insulating()),
+        OceanBC(Impenetrable(NoSlip()), Insulating()),
+        OceanBC(Penetrable(KinematicStress()), TemperatureFlux()),
+    )
+
+    test_ocean_gyre(imex = false, BC = BC, Δt = 600, refDat = refVals.explicit)
+
+    test_ocean_gyre(imex = true, BC = BC, Δt = 150, refDat = refVals.imex)
 end
diff --git a/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl b/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl
index 150ff8a164c..3a7e9ae24b7 100644
--- a/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl
+++ b/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl
@@ -7,6 +7,7 @@ using ClimateMachine.ODESolvers
 using ClimateMachine.Mesh.Filters
 using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
+using ClimateMachine.Ocean
 using ClimateMachine.Ocean.ShallowWater
 
 using ClimateMachine.Mesh.Topologies
@@ -62,7 +63,7 @@ function shallow_init_state!(
     return nothing
 end
 
-function shallow_init_aux!(p::SimpleBox, A, geom)
+function shallow_init_aux!(m::ShallowWaterModel, p::SimpleBox, A, geom)
     @inbounds y = geom.coord[2]
 
     FT = eltype(A)
@@ -102,6 +103,7 @@ function run_hydrostatic_spindown(; refDat = ())
     model_2D = ShallowWaterModel(
         param_set,
         prob_2D,
+        Uncoupled(),
         ShallowWater.ConstantViscosity{FT}(5e3),
         nothing,
         FT(1),
@@ -114,7 +116,7 @@ function run_hydrostatic_spindown(; refDat = ())
     dg_2D = DGModel(
         model_2D,
         grid_2D,
-        RusanovNumericalFlux(),
+        CentralNumericalFluxFirstOrder(),
         CentralNumericalFluxSecondOrder(),
         CentralNumericalFluxGradient(),
     )
@@ -202,7 +204,7 @@ function make_callbacks(
         if s
             starttime[] = now()
         else
-            energy = norm(Q_slow)
+            energy = norm(Q_fast)
             @info @sprintf(
                 """Update
                 simtime = %8.2f / %8.2f
diff --git a/test/Ocean/ShallowWater/GyreDriver.jl b/test/Ocean/ShallowWater/GyreDriver.jl
index 710c7d2b6d8..f736854c3ef 100644
--- a/test/Ocean/ShallowWater/GyreDriver.jl
+++ b/test/Ocean/ShallowWater/GyreDriver.jl
@@ -58,7 +58,14 @@ function setup_model(FT, stommel, linear, τₒ, fₒ, β, γ, ν, Lˣ, Lʸ, H)
         advection = NonLinearAdvection()
     end
 
-    model = ShallowWaterModel(param_set, problem, turbulence, advection, c)
+    model = ShallowWaterModel(
+        param_set,
+        problem,
+        Uncoupled(),
+        turbulence,
+        advection,
+        c,
+    )
 end
 
 function shallow_init_state!(
@@ -76,7 +83,7 @@ function shallow_init_state!(
     end
 end
 
-function shallow_init_aux!(p::GyreInABox, aux, geom)
+function shallow_init_aux!(m::ShallowWaterModel, p::GyreInABox, aux, geom)
     @inbounds y = geom.coord[2]
 
     Lʸ = p.Lʸ
diff --git a/test/Ocean/ShallowWater/GyreInABox.jl b/test/Ocean/ShallowWater/GyreInABox.jl
index bf3d81d0139..90d19f0d39d 100644
--- a/test/Ocean/ShallowWater/GyreInABox.jl
+++ b/test/Ocean/ShallowWater/GyreInABox.jl
@@ -8,6 +8,7 @@ using ClimateMachine.MPIStateArrays
 using ClimateMachine.ODESolvers
 using ClimateMachine.GenericCallbacks
 using ClimateMachine.VariableTemplates: flattenednames
+using ClimateMachine.Ocean
 using ClimateMachine.Ocean.ShallowWater
 using LinearAlgebra
 using StaticArrays
diff --git a/test/Ocean/SplitExplicit/hydrostatic_spindown.jl b/test/Ocean/SplitExplicit/hydrostatic_spindown.jl
index e8809de02dd..39ddfc8168a 100644
--- a/test/Ocean/SplitExplicit/hydrostatic_spindown.jl
+++ b/test/Ocean/SplitExplicit/hydrostatic_spindown.jl
@@ -7,6 +7,7 @@ using ClimateMachine.ODESolvers
 using ClimateMachine.Mesh.Filters
 using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
+using ClimateMachine.Ocean
 using ClimateMachine.Ocean.HydrostaticBoussinesq
 using ClimateMachine.Ocean.ShallowWater
 using ClimateMachine.Ocean.SplitExplicit: VerticalIntegralModel
@@ -67,7 +68,8 @@ function shallow_init_state!(
     return nothing
 end
 
-function shallow_init_aux!(p::GyreInABox, aux, geom)
+function shallow_init_aux!(m::ShallowWaterModel, p::GyreInABox, A, geom)
+    @inbounds x = geom.coord[1]
     @inbounds y = geom.coord[2]
 
     Lʸ = p.Lʸ
@@ -75,13 +77,20 @@ function shallow_init_aux!(p::GyreInABox, aux, geom)
     fₒ = p.fₒ
     β = p.β
 
-    aux.τ = @SVector [-τₒ * cos(π * y / Lʸ), 0]
-    aux.f = fₒ + β * (y - Lʸ / 2)
+    A.τ = @SVector [-τₒ * cos(π * y / Lʸ), 0]
+    A.f = fₒ + β * (y - Lʸ / 2)
+
+    A.Gᵁ = @SVector [-0, -0]
+    A.Δu = @SVector [-0, -0]
 
     return nothing
 end
 
-function run_hydrostatic_spindown(; coupled = true, refDat = ())
+function run_hydrostatic_spindown(;
+    coupling = Uncoupled(),
+    dt_slow = 300,
+    refDat = (),
+)
     mpicomm = MPI.COMM_WORLD
     ArrayType = ClimateMachine.array_type()
 
@@ -126,6 +135,7 @@ function run_hydrostatic_spindown(; coupled = true, refDat = ())
     model_3D = HydrostaticBoussinesqModel{FT}(
         param_set,
         prob_3D;
+        coupling = coupling,
         cʰ = FT(1),
         αᵀ = FT(0),
         κʰ = FT(0),
@@ -137,13 +147,13 @@ function run_hydrostatic_spindown(; coupled = true, refDat = ())
     model_2D = ShallowWaterModel(
         param_set,
         prob_2D,
+        coupling,
         ShallowWater.ConstantViscosity{FT}(model_3D.νʰ),
         nothing,
         FT(1),
     )
 
     dt_fast = 300
-    dt_slow = 300
     nout = ceil(Int64, tout / dt_slow)
     dt_slow = tout / nout
 
@@ -176,7 +186,7 @@ function run_hydrostatic_spindown(; coupled = true, refDat = ())
     dg_2D = DGModel(
         model_2D,
         grid_2D,
-        RusanovNumericalFlux(),
+        CentralNumericalFluxFirstOrder(),
         CentralNumericalFluxSecondOrder(),
         CentralNumericalFluxGradient(),
     )
@@ -186,8 +196,7 @@ function run_hydrostatic_spindown(; coupled = true, refDat = ())
 
     lsrk_3D = LSRK54CarpenterKennedy(dg_3D, Q_3D, dt = dt_slow, t0 = 0)
     lsrk_2D = LSRK54CarpenterKennedy(dg_2D, Q_2D, dt = dt_fast, t0 = 0)
-
-    odesolver = SplitExplicitSolver(lsrk_3D, lsrk_2D; coupled = coupled)
+    odesolver = SplitExplicitSolver(lsrk_3D, lsrk_2D)
 
     step = [0, 0]
     cbvector = make_callbacks(
@@ -308,6 +317,7 @@ function make_callbacks(
             (Q_slow, "3D state"),
             (dg_slow.state_auxiliary, "3D aux"),
             (Q_fast, "2D state"),
+            (dg_fast.state_auxiliary, "2D aux"),
         ],
         nout;
         prec = 12,
@@ -323,9 +333,9 @@ FT = Float64
 vtkpath = "vtk_split"
 
 const timeend = 24 * 3600 # s
-const tout = 2 * 3600 # s
+const tout = 3 * 3600 # s
 # const timeend = 1200 # s
-# const tout = 600 # s
+# const tout = 300 # s
 
 const N = 4
 const Nˣ = 5
@@ -346,5 +356,49 @@ const cᶻ = 0
 @testset "$(@__FILE__)" begin
     include("../refvals/hydrostatic_spindown_refvals.jl")
 
-    run_hydrostatic_spindown(coupled = false, refDat = refVals.uncoupled)
+    @testset "Multi-rate" begin
+        @testset "Δt = 30 mins" begin
+            run_hydrostatic_spindown(
+                coupling = Coupled(),
+                dt_slow = 30 * 60,
+                refDat = refVals.thirty_minutes,
+            )
+        end
+
+        @testset "Δt = 60 mins" begin
+            run_hydrostatic_spindown(
+                coupling = Coupled(),
+                dt_slow = 60 * 60,
+                refDat = refVals.sixty_minutes,
+            )
+        end
+
+        @testset "Δt = 90 mins" begin
+            run_hydrostatic_spindown(
+                coupling = Coupled(),
+                dt_slow = 90 * 60,
+                refDat = refVals.ninety_minutes,
+            )
+        end
+    end
+
+    if ClimateMachine.Settings.integration_testing
+        @testset "Single-Rate" begin
+            @testset "Not Coupled" begin
+                run_hydrostatic_spindown(
+                    coupling = Uncoupled(),
+                    dt_slow = 300,
+                    refDat = refVals.uncoupled,
+                )
+            end
+
+            @testset "Fully Coupled" begin
+                run_hydrostatic_spindown(
+                    coupling = Coupled(),
+                    dt_slow = 300,
+                    refDat = refVals.coupled,
+                )
+            end
+        end
+    end
 end
diff --git a/test/Ocean/SplitExplicit/test_vertical_integral_model.jl b/test/Ocean/SplitExplicit/test_vertical_integral_model.jl
index 2c388f2e0f6..d5c07492faf 100644
--- a/test/Ocean/SplitExplicit/test_vertical_integral_model.jl
+++ b/test/Ocean/SplitExplicit/test_vertical_integral_model.jl
@@ -7,6 +7,7 @@ using ClimateMachine.ODESolvers
 using ClimateMachine.Mesh.Filters
 using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
+using ClimateMachine.Ocean
 using ClimateMachine.Ocean.HydrostaticBoussinesq
 using ClimateMachine.Ocean.ShallowWater
 using ClimateMachine.Ocean.SplitExplicit: VerticalIntegralModel
@@ -67,7 +68,7 @@ function shallow_init_state!(
     return nothing
 end
 
-function shallow_init_aux!(p::GyreInABox, aux, geom)
+function shallow_init_aux!(m::ShallowWaterModel, p::GyreInABox, aux, geom)
     @inbounds y = geom.coord[2]
 
     Lʸ = p.Lʸ
@@ -131,6 +132,7 @@ function test_vertical_integral_model(time; refDat = ())
     model_2D = ShallowWaterModel(
         param_set,
         prob_2D,
+        Uncoupled(),
         ShallowWater.ConstantViscosity{FT}(model_3D.νʰ),
         nothing,
         FT(1),
@@ -157,7 +159,7 @@ function test_vertical_integral_model(time; refDat = ())
     dg_2D = DGModel(
         model_2D,
         grid_2D,
-        RusanovNumericalFlux(),
+        CentralNumericalFluxFirstOrder(),
         CentralNumericalFluxSecondOrder(),
         CentralNumericalFluxGradient(),
     )
diff --git a/test/Ocean/refvals/2D_hydrostatic_spindown_refvals.jl b/test/Ocean/refvals/2D_hydrostatic_spindown_refvals.jl
index ddddc247c8c..94ab62facdb 100644
--- a/test/Ocean/refvals/2D_hydrostatic_spindown_refvals.jl
+++ b/test/Ocean/refvals/2D_hydrostatic_spindown_refvals.jl
@@ -13,26 +13,26 @@ explicit = [
     [
         "2D state",
         :η,
-        -8.52722863677453757347e-01,
-        8.52828186384614328475e-01,
-        -2.32169838909612754587e-16,
-        6.03454974395413956678e-01,
+        -8.52722969951589915283e-01,
+        8.52846676313531282254e-01,
+        -2.49578135935735214742e-16,
+        6.03454239990563690021e-01,
     ],
     [
         "2D state",
         "U[1]",
-        -3.15402615071803076319e+01,
-        3.15402615071801157853e+01,
-        6.16398876385204608714e-15,
-        2.24269405191038053715e+01,
+        -3.15431401945821825450e+01,
+        3.15431401945818628008e+01,
+        6.11504145930918957291e-15,
+        2.24273815174625497093e+01,
     ],
     [
         "2D state",
         "U[2]",
-        -4.60204855830197230363e-13,
-        5.29714995818968166418e-13,
-        1.33362479920413591115e-14,
-        1.36912259776863805290e-13,
+        -7.62224398365580242501e-13,
+        9.72156930292624284356e-13,
+        1.39269607441935025982e-14,
+        1.95606703846656748360e-13,
     ],
 ]
 
diff --git a/test/Ocean/refvals/3D_hydrostatic_spindown_refvals.jl b/test/Ocean/refvals/3D_hydrostatic_spindown_refvals.jl
index 365b1e493e3..3f39f153b24 100644
--- a/test/Ocean/refvals/3D_hydrostatic_spindown_refvals.jl
+++ b/test/Ocean/refvals/3D_hydrostatic_spindown_refvals.jl
@@ -11,6 +11,10 @@ parr = [
     ["s_aux", :w, 12, 12, 0, 12],
     ["s_aux", :pkin, 15, 15, 15, 15],
     ["s_aux", :wz0, 12, 12, 0, 12],
+    ["aux", "uᵈ[1]", 15, 15, 15, 15],
+    ["aux", "uᵈ[2]", 15, 15, 15, 15],
+    ["aux", "ΔGᵘ[1]", 15, 15, 15, 15],
+    ["aux", "ΔGᵘ[2]", 15, 15, 15, 15],
 ]
 
 ### fully explicit
@@ -18,26 +22,26 @@ explicit = [
     [
         "state",
         "u[1]",
-        -9.58544066049459297929e-01,
-        9.58544066049460963264e-01,
-        3.63797880709171324926e-17,
-        4.45400263687296016357e-01,
+        -9.58544066049463849843e-01,
+        9.58544066049465071089e-01,
+        -6.13908923696726568442e-17,
+        4.45400263687296238402e-01,
     ],
     [
         "state",
         "u[2]",
-        -4.23271791659704063816e-14,
-        2.42232537903204158032e-14,
-        1.59523604488079777837e-16,
-        2.92759167504422236236e-15,
+        -4.33260855197780914020e-14,
+        1.99397359064900580713e-14,
+        -1.30343101521973575132e-16,
+        2.95525689323845820606e-15,
     ],
     [
         "state",
         :η,
-        -8.52732886154671576584e-01,
-        8.52845586939208977206e-01,
-        2.12048689718358208845e-14,
-        6.02992088522926961147e-01,
+        -8.52732886154656810618e-01,
+        8.52845586939211197652e-01,
+        2.20052243093959998331e-14,
+        6.02992088522925295813e-01,
     ],
     [
         "state",
@@ -58,10 +62,10 @@ explicit = [
     [
         "aux",
         :w,
-        -4.04553460063794881850e-04,
-        4.04714358463234059361e-04,
-        -1.55819801506140727111e-18,
-        1.63958655681888468020e-04,
+        -4.04553460063758398447e-04,
+        4.04714358463272711169e-04,
+        4.75730566051879549438e-19,
+        1.63958655681888576441e-04,
     ],
     [
         "aux",
@@ -74,10 +78,42 @@ explicit = [
     [
         "aux",
         :wz0,
-        -2.01164684799357424945e-04,
-        2.01041968159478722693e-04,
-        -4.04787314778332067465e-18,
-        1.42228420244457093172e-04,
+        -2.01164684799271339293e-04,
+        2.01041968159484089494e-04,
+        -2.10942374678779754294e-20,
+        1.42228420244455277133e-04,
+    ],
+    [
+        "aux",
+        "uᵈ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "aux",
+        "uᵈ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "aux",
+        "ΔGᵘ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "aux",
+        "ΔGᵘ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
     ],
 ]
 
@@ -85,26 +121,26 @@ imex = [
     [
         "state",
         "u[1]",
-        -9.57705359685807167125e-01,
-        9.57705359685806945080e-01,
-        1.59161572810262427692e-17,
-        4.45323727947873060362e-01,
+        -9.57705359685807500192e-01,
+        9.57705359685806500991e-01,
+        4.54747350886464140750e-17,
+        4.45323727947873004851e-01,
     ],
     [
         "state",
         "u[2]",
-        -4.14590421224419817494e-14,
-        2.60670000979756158115e-14,
-        1.43782448712293940635e-17,
-        2.67938064945467351808e-15,
+        -4.07596731389788922865e-14,
+        2.53992382675900717845e-14,
+        9.88087105439151860723e-18,
+        2.63742927063789745855e-15,
     ],
     [
         "state",
         :η,
-        -8.55941716778508276953e-01,
-        8.56014908798836349213e-01,
-        2.18005880014970914289e-14,
-        6.05260814296588622874e-01,
+        -8.55941716778505390373e-01,
+        8.56014908798837681481e-01,
+        2.16732587432488803528e-14,
+        6.05260814296588400829e-01,
     ],
     [
         "state",
@@ -125,10 +161,10 @@ imex = [
     [
         "aux",
         :w,
-        -4.03419301824631281773e-04,
-        4.03626959596156903790e-04,
-        2.77555756156289140884e-20,
-        1.63816237485524170961e-04,
+        -4.03419301824637895407e-04,
+        4.03626959596159180614e-04,
+        -3.33066907387546970565e-21,
+        1.63816237485524198066e-04,
     ],
     [
         "aux",
@@ -141,10 +177,42 @@ imex = [
     [
         "aux",
         :wz0,
-        -1.97016630922658897082e-04,
-        1.97066086934920417029e-04,
-        -5.08482145278321681062e-19,
-        1.39393318574472491710e-04,
+        -1.97016630922629948884e-04,
+        1.97066086934914155761e-04,
+        -7.49400541621980656312e-19,
+        1.39393318574472925390e-04,
+    ],
+    [
+        "aux",
+        "uᵈ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "aux",
+        "uᵈ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "aux",
+        "ΔGᵘ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "aux",
+        "ΔGᵘ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
     ],
 ]
 
diff --git a/test/Ocean/refvals/hydrostatic_spindown_refvals.jl b/test/Ocean/refvals/hydrostatic_spindown_refvals.jl
index c21b414d7f6..84d27edff8f 100644
--- a/test/Ocean/refvals/hydrostatic_spindown_refvals.jl
+++ b/test/Ocean/refvals/hydrostatic_spindown_refvals.jl
@@ -2,44 +2,56 @@
 #  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
 #  [         :                :          :        :      :          :           ],
 # ]
+
 parr = [
-    ["Q", "u[1]", 12, 12, 0, 12],
-    ["Q", "u[2]", 0, 0, 0, 0],
-    ["Q", :η, 12, 12, 0, 12],
-    ["Q", :θ, 15, 15, 15, 15],
-    ["s_aux", :y, 15, 15, 15, 15],
-    ["s_aux", :w, 12, 12, 0, 12],
-    ["s_aux", :pkin, 15, 15, 15, 15],
-    ["s_aux", :wz0, 12, 12, 0, 12],
+    ["3D state", "u[1]", 12, 12, 0, 12],
+    ["3D state", "u[2]", 0, 0, 0, 0],
+    ["3D state", :η, 12, 12, 0, 12],
+    ["3D state", :θ, 15, 15, 15, 15],
+    ["3D aux", :y, 15, 15, 15, 15],
+    ["3D aux", :w, 12, 12, 0, 12],
+    ["3D aux", :pkin, 15, 15, 15, 15],
+    ["3D aux", :wz0, 12, 12, 0, 12],
+    ["3D aux", "uᵈ[1]", 12, 12, 0, 12],
+    ["3D aux", "uᵈ[2]", 0, 0, 0, 0],
+    ["3D aux", "ΔGᵘ[1]", 11, 12, 0, 12],
+    ["3D aux", "ΔGᵘ[2]", 0, 0, 0, 0],
     ["2D state", :η, 12, 12, 0, 12],
     ["2D state", "U[1]", 12, 12, 0, 12],
     ["2D state", "U[2]", 0, 0, 0, 0],
+    ["2D aux", :f, 15, 15, 15, 15],
+    ["2D aux", "τ[1]", 15, 15, 15, 15],
+    ["2D aux", "τ[2]", 15, 15, 15, 15],
+    ["2D aux", "Gᵁ[1]", 12, 11, 0, 12],
+    ["2D aux", "Gᵁ[2]", 0, 0, 0, 0],
+    ["2D aux", "Δu[1]", 12, 12, 0, 12],
+    ["2D aux", "Δu[2]", 0, 0, 0, 0],
 ]
 
 uncoupled = [
     [
         "3D state",
         "u[1]",
-        -9.58547428246829258391e-01,
-        9.58547428246829924525e-01,
-        -2.72848410531878478287e-17,
+        -9.58547428246830479637e-01,
+        9.58547428246829258391e-01,
+        -6.82121026329696195717e-18,
         4.45400655325317695876e-01,
     ],
     [
         "3D state",
         "u[2]",
-        -4.17338353173569118678e-14,
-        1.83725063946779860804e-14,
-        2.10967530406246663799e-17,
-        2.12845310693392360741e-15,
+        -4.13574641907576004779e-14,
+        1.84325494853789287544e-14,
+        1.94567057477529826177e-17,
+        2.12815907156924710074e-15,
     ],
     [
         "3D state",
         :η,
-        -8.52721734395888830704e-01,
+        -8.52721734395889496838e-01,
         8.52834259382379111791e-01,
-        2.16959961107932042100e-14,
-        6.02983573168680564436e-01,
+        2.18869899981655176688e-14,
+        6.02983573168680453414e-01,
     ],
     [
         "3D state",
@@ -60,10 +72,10 @@ uncoupled = [
     [
         "3D aux",
         :w,
-        -4.04489222828921936130e-04,
-        4.04649979048339884964e-04,
-        6.73350264435157450488e-19,
-        1.63949585575251874240e-04,
+        -4.04489222828922911912e-04,
+        4.04649979048336849198e-04,
+        7.79376563286859913379e-19,
+        1.63949585575251901345e-04,
     ],
     [
         "3D aux",
@@ -76,35 +88,845 @@ uncoupled = [
     [
         "3D aux",
         :wz0,
-        -2.00933099660539595454e-04,
-        2.00809859469747143366e-04,
-        9.99200722162640859033e-20,
-        1.42064682697776525650e-04,
+        -2.00933099660541601228e-04,
+        2.00809859469744893647e-04,
+        2.62012633811536956698e-19,
+        1.42064682697776552755e-04,
+    ],
+    [
+        "3D aux",
+        "uᵈ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        "uᵈ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        "ΔGᵘ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        "ΔGᵘ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
     ],
     [
         "2D state",
         :η,
-        -8.52735654368523721125e-01,
-        8.52839610988160190530e-01,
-        -5.46940270851337130538e-16,
-        6.03463847107759909782e-01,
+        -8.52715894965927923010e-01,
+        8.52861218874714666072e-01,
+        -7.70938868299708696148e-16,
+        6.03458382295490314284e-01,
     ],
     [
         "2D state",
         "U[1]",
-        -3.15391722007734536248e+01,
-        3.15391722007730805899e+01,
-        4.59412508035939049964e-15,
-        2.24261366371601056358e+01,
+        -3.15423945317149438949e+01,
+        3.15423945317147342848e+01,
+        1.42724412149908289268e-14,
+        2.24262219407601044452e+01,
     ],
     [
         "2D state",
         "U[2]",
-        -4.37584346168357661464e-13,
-        5.56621403621539302858e-13,
-        1.62907035110957408041e-14,
-        1.42652658850586700348e-13,
+        -9.38673083374999972374e-13,
+        1.16534781393478256382e-12,
+        1.93385973223830695204e-14,
+        2.16634164433978071172e-13,
+    ],
+    [
+        "2D aux",
+        :f,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "τ[1]",
+        -0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "τ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "Gᵁ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "Gᵁ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "Δu[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "Δu[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+]
+
+coupled = [
+    [
+        "3D state",
+        "u[1]",
+        -9.58544854429326798062e-01,
+        9.58544854429326798062e-01,
+        -4.09272615797817732838e-17,
+        4.45400401208546903309e-01,
+    ],
+    [
+        "3D state",
+        "u[2]",
+        -4.09119364541065205420e-15,
+        3.41954410407034240802e-15,
+        1.57855184364031418231e-19,
+        7.69982545727116047928e-16,
+    ],
+    [
+        "3D state",
+        :η,
+        -8.52733075123627615177e-01,
+        8.52843573070954374948e-01,
+        7.27595761418342649852e-17,
+        6.02992194663175995473e-01,
+    ],
+    [
+        "3D state",
+        :θ,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        :y,
+        0.00000000000000000000e+00,
+        1.00000000000000011642e+06,
+        5.00000000000000000000e+05,
+        2.92775877460665535182e+05,
+    ],
+    [
+        "3D aux",
+        :w,
+        -4.04456725919283483547e-04,
+        4.04616610112968017130e-04,
+        1.33892896769793873666e-18,
+        1.63945200913093809460e-04,
+    ],
+    [
+        "3D aux",
+        :pkin,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        :wz0,
+        -2.00813442787548966564e-04,
+        2.00686991710668984502e-04,
+        1.68642877440561282675e-18,
+        1.41979510156587984725e-04,
+    ],
+    [
+        "3D aux",
+        "uᵈ[1]",
+        -8.79713597641291755735e-01,
+        8.79713597641291200624e-01,
+        -9.26547727431170607429e-17,
+        4.41871673548253185437e-01,
+    ],
+    [
+        "3D aux",
+        "uᵈ[2]",
+        -5.88884665747485312740e-28,
+        6.47654803186450692243e-28,
+        1.40221555861782666734e-30,
+        4.22950073030020620100e-29,
+    ],
+    [
+        "3D aux",
+        "ΔGᵘ[1]",
+        -1.50292385646252865781e-09,
+        1.50292385640701573972e-09,
+        -3.03178502810731714103e-21,
+        6.72141424487058536370e-10,
+    ],
+    [
+        "3D aux",
+        "ΔGᵘ[2]",
+        -2.50330922252795548175e-19,
+        2.79421319642066487612e-19,
+        8.87468518373638336981e-36,
+        5.19465794381941235415e-20,
+    ],
+    [
+        "2D state",
+        :η,
+        -8.52733075123627615177e-01,
+        8.52843573070954374948e-01,
+        1.61115565333602722195e-16,
+        6.03463098440266798583e-01,
+    ],
+    [
+        "2D state",
+        "U[1]",
+        -3.15402749945522060671e+01,
+        3.15402749945527389741e+01,
+        2.12974776703234167348e-14,
+        2.24262621677375086904e+01,
+    ],
+    [
+        "2D state",
+        "U[2]",
+        -1.63647745816416607002e-12,
+        1.36781764162808076476e-12,
+        6.31420737450488849608e-17,
+        3.08233543917742288569e-13,
+    ],
+    [
+        "2D aux",
+        :f,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "τ[1]",
+        -0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "τ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "Gᵁ[1]",
+        -6.01169542562806285963e-07,
+        6.01169542585011466433e-07,
+        1.21270571032004377760e-18,
+        2.69066532005499711718e-07,
+    ],
+    [
+        "2D aux",
+        "Gᵁ[2]",
+        -1.11768527856826595815e-16,
+        1.00132368901118217729e-16,
+        -3.50057026691823957451e-33,
+        2.07948587451660775232e-17,
+    ],
+    [
+        "2D aux",
+        "Δu[1]",
+        -6.53936429129693404076e-04,
+        6.53936429129569046087e-04,
+        -2.70378337774435083119e-18,
+        4.64975945389895016328e-04,
+    ],
+    [
+        "2D aux",
+        "Δu[2]",
+        -3.76049176138251460297e-17,
+        3.57641761985933683411e-17,
+        -4.12443309266433447246e-19,
+        7.09162842149286857838e-18,
+    ],
+]
+
+thirty_minutes = [
+    [
+        "3D state",
+        "u[1]",
+        -9.58527737426365766815e-01,
+        9.58527737426364989659e-01,
+        -2.95585778076201651428e-17,
+        4.45397585011924779241e-01,
+    ],
+    [
+        "3D state",
+        "u[2]",
+        -3.91776258115697290002e-15,
+        3.71573344368902785991e-15,
+        9.99453214450878842160e-18,
+        7.55533641450465537744e-16,
+    ],
+    [
+        "3D state",
+        :η,
+        -8.52729582431237531637e-01,
+        8.52830944991449846349e-01,
+        -2.54658516496419884307e-16,
+        6.02990141900936249542e-01,
+    ],
+    [
+        "3D state",
+        :θ,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        :y,
+        0.00000000000000000000e+00,
+        1.00000000000000011642e+06,
+        5.00000000000000000000e+05,
+        2.92775877460665535182e+05,
+    ],
+    [
+        "3D aux",
+        :w,
+        -4.06712181304362465229e-04,
+        4.06869755985816939844e-04,
+        9.07052211118752913371e-19,
+        1.64290385445136678539e-04,
+    ],
+    [
+        "3D aux",
+        :pkin,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        :wz0,
+        -2.09003099782708142351e-04,
+        2.08868143992394606845e-04,
+        1.33781874467331362536e-18,
+        1.47768608464728741676e-04,
+    ],
+    [
+        "3D aux",
+        "uᵈ[1]",
+        -8.79792191807299950312e-01,
+        8.79792191807298729067e-01,
+        -9.54969436861574689412e-17,
+        4.41911149114936674387e-01,
+    ],
+    [
+        "3D aux",
+        "uᵈ[2]",
+        -4.49650715975976729085e-28,
+        3.22939933074851707840e-28,
+        1.28559371735991545201e-30,
+        3.12162656228539824003e-29,
+    ],
+    [
+        "3D aux",
+        "ΔGᵘ[1]",
+        -1.52356619979685413105e-09,
+        1.52356619973800994984e-09,
+        -3.00415481336788178489e-21,
+        6.81373145706173027634e-10,
+    ],
+    [
+        "3D aux",
+        "ΔGᵘ[2]",
+        -2.17115289587171316651e-19,
+        2.99271337844227832879e-19,
+        3.57452597678270976423e-36,
+        5.23982958192948877053e-20,
+    ],
+    [
+        "2D state",
+        :η,
+        -8.52729582431237531637e-01,
+        8.52830944991449846349e-01,
+        -3.20454773827805192348e-16,
+        6.03461044074932395631e-01,
+    ],
+    [
+        "2D state",
+        "U[1]",
+        -3.15407643291233448224e+01,
+        3.15407643291237640426e+01,
+        2.30678864898692384736e-14,
+        2.24265574817805202201e+01,
+    ],
+    [
+        "2D state",
+        "U[2]",
+        -1.56710503246274425174e-12,
+        1.48629337747559491236e-12,
+        3.99781285780299048032e-15,
+        3.02449468686897625244e-13,
+    ],
+    [
+        "2D aux",
+        :f,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "τ[1]",
+        -0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "τ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "Gᵁ[1]",
+        -6.09426479895203986555e-07,
+        6.09426479918741655731e-07,
+        1.20165633492970096897e-18,
+        2.72762104279987029506e-07,
+    ],
+    [
+        "2D aux",
+        "Gᵁ[2]",
+        -1.19708535137691129300e-16,
+        8.68461158348685235787e-17,
+        -1.50869648123518501517e-33,
+        2.09756864038773538444e-17,
+    ],
+    [
+        "2D aux",
+        "Δu[1]",
+        -3.89449877638695009935e-03,
+        3.89449877638627485824e-03,
+        -2.02396849009311295664e-17,
+        2.76903473119681107703e-03,
+    ],
+    [
+        "2D aux",
+        "Δu[2]",
+        -2.13991366909841745553e-16,
+        2.02991354131742782465e-16,
+        -1.26781868574011981865e-18,
+        4.34606820130057938426e-17,
+    ],
+]
+
+sixty_minutes = [
+    [
+        "3D state",
+        "u[1]",
+        -9.58507705628042327994e-01,
+        9.58507705628042550039e-01,
+        -5.22959453519433752618e-17,
+        4.45394100308371843067e-01,
+    ],
+    [
+        "3D state",
+        "u[2]",
+        -4.09893498131451878178e-15,
+        3.66595063838075513859e-15,
+        3.52058149555652693147e-18,
+        7.65972478919525466697e-16,
+    ],
+    [
+        "3D state",
+        :η,
+        -8.52715753705294066123e-01,
+        8.52819886810596838878e-01,
+        0.00000000000000000000e+00,
+        6.02986351734545844572e-01,
+    ],
+    [
+        "3D state",
+        :θ,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        :y,
+        0.00000000000000000000e+00,
+        1.00000000000000011642e+06,
+        5.00000000000000000000e+05,
+        2.92775877460665535182e+05,
+    ],
+    [
+        "3D aux",
+        :w,
+        -4.09355343664679748733e-04,
+        4.09526445665209497954e-04,
+        8.01581023779363006131e-19,
+        1.64782726593626366188e-04,
+    ],
+    [
+        "3D aux",
+        :pkin,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        :wz0,
+        -2.18597874938533713813e-04,
+        2.18512725277962983460e-04,
+        1.16684439888103955116e-18,
+        1.54581525183602542031e-04,
+    ],
+    [
+        "3D aux",
+        "uᵈ[1]",
+        -8.79886113674857694988e-01,
+        8.79886113674856806810e-01,
+        -8.81072992342524199517e-17,
+        4.41958323665245178535e-01,
+    ],
+    [
+        "3D aux",
+        "uᵈ[2]",
+        -5.04476548888837049561e-28,
+        4.09418809809705127009e-28,
+        1.53142189509630563206e-30,
+        3.27579623353337388378e-29,
+    ],
+    [
+        "3D aux",
+        "ΔGᵘ[1]",
+        -1.12461707299293835875e-09,
+        1.12461707291060059378e-09,
+        -2.89803852573586336546e-21,
+        5.02954103857782471789e-10,
+    ],
+    [
+        "3D aux",
+        "ΔGᵘ[2]",
+        -2.29220407034248886795e-19,
+        3.23795504633739099693e-19,
+        4.19082355898662551731e-36,
+        5.50021671320874203932e-20,
+    ],
+    [
+        "2D state",
+        :η,
+        -8.52715753705294066123e-01,
+        8.52819886810596838878e-01,
+        -2.04281036531028799987e-17,
+        6.03457250948630341547e-01,
+    ],
+    [
+        "2D state",
+        "U[1]",
+        -3.15415180428826182890e+01,
+        3.15415180428829948767e+01,
+        1.77902803599749859788e-14,
+        2.24265293235027023400e+01,
+    ],
+    [
+        "2D state",
+        "U[2]",
+        -1.63957399252576030728e-12,
+        1.46638025535227397199e-12,
+        1.40823259822193190848e-15,
+        3.06628264537953242560e-13,
+    ],
+    [
+        "2D aux",
+        :f,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "τ[1]",
+        -0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "τ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "Gᵁ[1]",
+        -4.49846829164240224277e-07,
+        4.49846829197175353425e-07,
+        1.15921837490966064206e-18,
+        2.01338753352722549629e-07,
+    ],
+    [
+        "2D aux",
+        "Gᵁ[2]",
+        -1.29518201853495637566e-16,
+        9.16881628136995516365e-17,
+        -2.08062063752041870574e-33,
+        2.20180483211723137251e-17,
+    ],
+    [
+        "2D aux",
+        "Δu[1]",
+        -7.71677156339491132631e-03,
+        7.71677156339359640591e-03,
+        -2.97775205101297201914e-17,
+        5.48673197909435757247e-03,
+    ],
+    [
+        "2D aux",
+        "Δu[2]",
+        -5.01794629524824969143e-16,
+        5.51655281127317078031e-16,
+        3.20748218455203547025e-18,
+        9.36761976343254920594e-17,
+    ],
+]
+
+ninety_minutes = [
+    [
+        "3D state",
+        "u[1]",
+        -9.58488051495743564878e-01,
+        9.58488051495744008967e-01,
+        -2.72848410531878478287e-17,
+        4.45390687454859546257e-01,
+    ],
+    [
+        "3D state",
+        "u[2]",
+        -4.05044403917298945083e-15,
+        3.37351039462521829098e-15,
+        2.85851376495618499592e-18,
+        7.53889134369023261669e-16,
+    ],
+    [
+        "3D state",
+        :η,
+        -8.52731405452108237597e-01,
+        8.52809353848989259994e-01,
+        -1.36424205265939223736e-16,
+        6.02988963894262375298e-01,
+    ],
+    [
+        "3D state",
+        :θ,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        :y,
+        0.00000000000000000000e+00,
+        1.00000000000000011642e+06,
+        5.00000000000000000000e+05,
+        2.92775877460665535182e+05,
+    ],
+    [
+        "3D aux",
+        :w,
+        -4.11959288170114497438e-04,
+        4.12126623290113712657e-04,
+        8.97615315409439082099e-19,
+        1.65355949740610678371e-04,
+    ],
+    [
+        "3D aux",
+        :pkin,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "3D aux",
+        :wz0,
+        -2.28048874140054377109e-04,
+        2.27949849600823879541e-04,
+        1.39666056497844683559e-18,
+        1.61271245461141397555e-04,
+    ],
+    [
+        "3D aux",
+        "uᵈ[1]",
+        -8.79979607451059520073e-01,
+        8.79979607451058964962e-01,
+        -9.32232069317251423826e-17,
+        4.42005283454923125763e-01,
+    ],
+    [
+        "3D aux",
+        "uᵈ[2]",
+        -5.48652759581213710664e-28,
+        3.92852730800063879095e-28,
+        1.36055159257466689218e-30,
+        3.58669418355622340063e-29,
+    ],
+    [
+        "3D aux",
+        "ΔGᵘ[1]",
+        -8.86991425948811828542e-10,
+        8.86991425888052000430e-10,
+        -2.87311034609816924955e-21,
+        3.96682558396691740283e-10,
+    ],
+    [
+        "3D aux",
+        "ΔGᵘ[2]",
+        -2.66880587297174792697e-19,
+        3.32595187063539502433e-19,
+        3.45126646034192634633e-36,
+        5.71491963745139874533e-20,
+    ],
+    [
+        "2D state",
+        :η,
+        -8.52731405452108237597e-01,
+        8.52809353848989259994e-01,
+        -1.20259358027396949872e-16,
+        6.03459865148300300675e-01,
+    ],
+    [
+        "2D state",
+        "U[1]",
+        -3.15423820746578229546e+01,
+        3.15423820746582954655e+01,
+        2.75259190440912945341e-14,
+        2.24266774102465760166e+01,
+    ],
+    [
+        "2D state",
+        "U[2]",
+        -1.62017761566911967103e-12,
+        1.34940415784996204528e-12,
+        1.14340550598186518880e-15,
+        3.01791152146446424566e-13,
+    ],
+    [
+        "2D aux",
+        :f,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "τ[1]",
+        -0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "τ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "2D aux",
+        "Gᵁ[1]",
+        -3.54796570355220815061e-07,
+        3.54796570379524731417e-07,
+        1.14925108410943513278e-18,
+        1.58796938272805492033e-07,
+    ],
+    [
+        "2D aux",
+        "Gᵁ[2]",
+        -1.33038074825415802514e-16,
+        1.06752234918869909375e-16,
+        -1.45939267465887190460e-33,
+        2.28775307028972280128e-17,
+    ],
+    [
+        "2D aux",
+        "Δu[1]",
+        -1.14629123365995507638e-02,
+        1.14629123365980502280e-02,
+        -3.95057006856089328375e-17,
+        8.15149713200660763768e-03,
+    ],
+    [
+        "2D aux",
+        "Δu[2]",
+        -7.52467366304891760665e-16,
+        5.94961001279223574512e-16,
+        6.25849006198159897312e-19,
+        1.42516297748358802032e-16,
     ],
 ]
 
-refVals = (uncoupled = (uncoupled, parr),)
+refVals = (
+    uncoupled = (uncoupled, parr),
+    coupled = (coupled, parr),
+    thirty_minutes = (thirty_minutes, parr),
+    sixty_minutes = (sixty_minutes, parr),
+    ninety_minutes = (ninety_minutes, parr),
+)
diff --git a/test/Ocean/refvals/test_ocean_gyre_refvals.jl b/test/Ocean/refvals/test_ocean_gyre_refvals.jl
index adc58466ea6..e3b2bb6bde7 100644
--- a/test/Ocean/refvals/test_ocean_gyre_refvals.jl
+++ b/test/Ocean/refvals/test_ocean_gyre_refvals.jl
@@ -1,45 +1,19 @@
-# Testing reference values and precisions
-# Each test block of varr and parr should be followed by an append to refVals, refPrecs arrays.
-# e.g.
-#   refVals=[]
-#   refPrecs=[]
-#
-#   varr = ..........
-#   par  = ..........
-#
-#   append!(refVals ,[ varr ] )
-#   append!(refPrecs,[ parr ] )
-#
-#   varr = ..........
-#   par  = ..........
-#
-#   append!(refVals ,[ varr ] )
-#   append!(refPrecs,[ parr ] )
-#
-#   varr = ..........
-#   par  = ..........
-#
-#   append!(refVals ,[ varr ] )
-#   append!(refPrecs,[ parr ] )
-#
-#   etc.....
-#
-#   Now for real!
-#
+parr = [
+    ["Q", "u[1]", 12, 12, 12, 12],
+    ["Q", "u[2]", 12, 12, 12, 12],
+    ["Q", :η, 12, 12, 12, 12],
+    ["Q", :θ, 12, 12, 12, 12],
+    ["s_aux", :y, 12, 12, 12, 12],
+    ["s_aux", :w, 12, 12, 12, 12],
+    ["s_aux", :pkin, 12, 12, 12, 12],
+    ["s_aux", :wz0, 12, 12, 12, 12],
+    ["s_aux", "uᵈ[1]", 12, 12, 12, 12],
+    ["s_aux", "uᵈ[2]", 12, 12, 12, 12],
+    ["s_aux", "ΔGᵘ[1]", 12, 12, 12, 12],
+    ["s_aux", "ΔGᵘ[2]", 12, 12, 12, 12],
+]
 
-refVals = []
-refPrecs = []
-# SC ========== Test number 1 reference values and precision match template. =======
-# SC ========== /Users/chrishill/projects/clima/cm/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl test reference values ======================================
-# BEGIN SCPRINT
-# varr - reference values (from reference run)
-# parr - digits match precision (hand edit as needed)
-#
-# [
-#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
-#  [         :                :          :        :      :          :           ],
-# ]
-varr = [
+explicit = [
     [
         "Q",
         "u[1]",
@@ -104,33 +78,41 @@ varr = [
         -2.78589672590168561413e-08,
         7.72520757253795920232e-06,
     ],
+    [
+        "s_aux",
+        "uᵈ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "uᵈ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "ΔGᵘ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "ΔGᵘ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
 ]
-parr = [
-    ["Q", "u[1]", 12, 12, 12, 12],
-    ["Q", "u[2]", 12, 12, 12, 12],
-    ["Q", :η, 12, 12, 12, 12],
-    ["Q", :θ, 12, 12, 12, 12],
-    ["s_aux", :y, 12, 12, 12, 12],
-    ["s_aux", :w, 12, 12, 12, 12],
-    ["s_aux", :pkin, 12, 12, 12, 12],
-    ["s_aux", :wz0, 12, 12, 12, 12],
-]
-# END SCPRINT
-
-append!(refVals, [varr])
-append!(refPrecs, [parr])
 
-# SC ========== Test number 2 reference values and precision match template. =======
-# SC ========== /Users/chrishill/projects/clima/cm/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl test reference values ======================================
-# BEGIN SCPRINT
-# varr - reference values (from reference run)
-# parr - digits match precision (hand edit as needed)
-#
-# [
-#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
-#  [         :                :          :        :      :          :           ],
-# ]
-varr = [
+imex = [
     [
         "Q",
         "u[1]",
@@ -195,17 +177,38 @@ varr = [
         -2.71999565492576589218e-08,
         7.76164707621753630128e-06,
     ],
-]
-parr = [
-    ["Q", "u[1]", 12, 12, 12, 12],
-    ["Q", "u[2]", 12, 12, 12, 12],
-    ["Q", :η, 12, 12, 12, 12],
-    ["Q", :θ, 12, 12, 12, 12],
-    ["s_aux", :y, 12, 12, 12, 12],
-    ["s_aux", :w, 12, 12, 12, 12],
-    ["s_aux", :pkin, 12, 12, 12, 12],
-    ["s_aux", :wz0, 12, 12, 12, 12],
+    [
+        "s_aux",
+        "uᵈ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "uᵈ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "ΔGᵘ[1]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
+    [
+        "s_aux",
+        "ΔGᵘ[2]",
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+        0.00000000000000000000e+00,
+    ],
 ]
 
-append!(refVals, [varr])
-append!(refPrecs, [parr])
+refVals = (explicit = (explicit, parr), imex = (imex, parr))