diff --git a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
index ca884961772..46c6ab8f253 100644
--- a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
+++ b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
@@ -14,9 +14,9 @@ using ..DGMethods:
     BalanceLaw,
     LocalGeometry,
     DGModel,
+    nodal_update_auxiliary_state!,
     indefinite_stack_integral!,
     reverse_indefinite_stack_integral!,
-    nodal_update_auxiliary_state!,
     copy_stack_field_down!
 using ..DGMethods.NumericalFluxes: RusanovNumericalFlux
 
diff --git a/src/Ocean/SplitExplicit/Communication.jl b/src/Ocean/SplitExplicit/Communication.jl
new file mode 100644
index 00000000000..73ca2a0f96a
--- /dev/null
+++ b/src/Ocean/SplitExplicit/Communication.jl
@@ -0,0 +1,139 @@
+import ..DGMethods:
+    initialize_states!,
+    tendency_from_slow_to_fast!,
+    cummulate_fast_solution!,
+    reconcile_from_fast_to_slow!
+
+@inline function initialize_states!(
+    baroclinic::HydrostaticBoussinesqModel,
+    barotropic::ShallowWaterModel,
+    model_bc,
+    model_bt,
+    state_bc,
+    state_bt,
+)
+    model_bt.state_auxiliary.ηᶜ .= -0
+    model_bt.state_auxiliary.Uᶜ .= (@SVector [-0, -0, -0])'
+
+    state_bt.η .= model_bt.state_auxiliary.ηˢ
+    state_bt.U .= model_bt.state_auxiliary.Uˢ
+
+    model_bc.state_auxiliary.ΔGᵘ .= -0
+
+    return nothing
+end
+
+@inline function tendency_from_slow_to_fast!(
+    baroclinic::HydrostaticBoussinesqModel,
+    barotropic::ShallowWaterModel,
+    model_bc,
+    model_bt,
+    state_bc,
+    state_bt,
+    forcing_tendency,
+)
+    # integrate the tendency
+    tendency_model = model_bc.modeldata.tendency_model
+    update_auxiliary_state!(
+        tendency_model,
+        tendency_model.balance_law,
+        forcing_tendency,
+        0,
+    )
+
+    Nq, Nqk, _, _, nelemv, _, nelemh, _ = basic_grid_info(model_bc)
+
+    ## get top value (=integral over full depth) of ∫du
+    ∫du = tendency_model.state_auxiliary.∫du
+    boxy_∫du = reshape(∫du, Nq^2, Nqk, 2, nelemv, nelemh)
+    flat_∫du = @view boxy_∫du[:, end, :, end, :]
+
+    ## copy into Gᵁ of dgFast
+    model_bt.state_auxiliary.Gᵁ .= reshape(flat_∫du, Nq^2, 2, nelemh)
+
+    ## scale by -1/H and copy back to ΔGu
+    # note: since tendency_dg.auxstate.∫du is not used after this, could be
+    # re-used to store a 3-D copy of "-Gu"
+    ΔGᵘ = model_bc.state_auxiliary.ΔGᵘ
+    boxy_ΔGᵘ = reshape(ΔGᵘ, Nq^2, Nqk, 2, nelemv, nelemh)
+    boxy_ΔGᵘ .= -reshape(flat_∫du, Nq^2, 1, 2, 1, nelemh) / baroclinic.problem.H
+
+    return nothing
+end
+
+@inline function cummulate_fast_solution!(
+    baroclinic::HydrostaticBoussinesqModel,
+    barotropic::ShallowWaterModel,
+    model_bt,
+    state_bt,
+    fast_time,
+    fast_dt,
+    substep,
+)
+    ### for now no averaging, just save the final time step
+    model_bt.state_auxiliary.ηᶜ .= state_bt.η
+    model_bt.state_auxiliary.Uᶜ .= state_bt.U
+
+    ### technically unnecessary right now, but for future changes
+    model_bt.state_auxiliary.ηˢ .= state_bt.η
+    model_bt.state_auxiliary.Uˢ .= state_bt.U
+
+    return nothing
+end
+
+@inline function reconcile_from_fast_to_slow!(
+    baroclinic::HydrostaticBoussinesqModel,
+    barotropic::ShallowWaterModel,
+    model_bc,
+    model_bt,
+    state_bc,
+    state_bt,
+)
+    Nq, Nqk, _, _, nelemv, _, nelemh, _ = basic_grid_info(model_bc)
+
+    # need to calculate int_u using integral kernels
+    # u_slow := u_slow + (1/H) * (u_fast - \int_{-H}^{0} u_slow)
+
+    # Compute: \int_{-H}^{0} u_slow)
+    ### need to make sure this is stored into aux.∫u
+
+    # integrate vertically horizontal velocity
+    flow_model = model_bc.modeldata.flow_model
+    update_auxiliary_state!(flow_model, flow_model.balance_law, state_bc, 0)
+
+    ## get top value (=integral over full depth)
+    ∫u = flow_model.state_auxiliary.∫u
+    boxy_∫u = reshape(∫u, Nq^2, Nqk, 2, nelemv, nelemh)
+    flat_∫u = @view boxy_∫u[:, end, :, end, :]
+
+    ## get time weighted averaged out of cumulative arrays
+    ## no weights right now
+    state_bt.state_auxiliary.Uᶜ .*= 1 # / fast_time_rec[2]
+    state_bt.state_auxiliary.ηᶜ .*= 1 # / fast_time_rec[2]
+
+    ### substract ∫u from U and divide by H
+
+    ### Δu is a place holder for 1/H * (Ū - ∫u)
+    Δu = model_bt.state_auxiliary.Δu
+    Δu .= 1 / baroclinic.problem.H * (model_bt.state_auxiliary.Uᶜ - flat_∫u)
+
+    ### copy the 2D contribution down the 3D solution
+    ### need to reshape these things for the broadcast
+    boxy_u = reshape(state_bc.u, Nq^2, Nqk, 2, nelemv, nelemh)
+    boxy_Δu = reshape(Δu, Nq^2, 1, 3, 1, nelemh)
+
+    ### copy 2D eta over to 3D model
+    boxy_η_3D = reshape(state_bc.η, Nq^2, Nqk, nelemv, nelemh)
+    boxy_η_2D = reshape(model_bt.state_auxiliary.ηᶜ, Nq^2, 1, 1, nelemh)
+
+    ### this works, we tested it
+    # currently SWModel has a 3-vector for U
+    # need to shrink to 2-vector at some point
+
+    #### turn off broadcast for testing purposes
+    #### probably need to add a flag here
+    # boxy_u .+= boxy_Δu[:, :, 1:2, :, :]
+    # boxy_η_3D .= boxy_η_2D
+
+    return nothing
+end
diff --git a/src/Ocean/SplitExplicit/FlowIntegralModel.jl b/src/Ocean/SplitExplicit/FlowIntegralModel.jl
new file mode 100644
index 00000000000..a9cfc1fa502
--- /dev/null
+++ b/src/Ocean/SplitExplicit/FlowIntegralModel.jl
@@ -0,0 +1,62 @@
+struct FlowIntegralModel{M} <: BalanceLaw
+    ocean::M
+    function FlowIntegralModel(ocean::M) where {M}
+        return new{M}(ocean)
+    end
+end
+
+vars_state_gradient(tm::FlowIntegralModel, FT) = @vars()
+vars_state_gradient_flux(tm::FlowIntegralModel, FT) = @vars()
+
+vars_state_conservative(tm::FlowIntegralModel, FT) =
+    vars_state_conservative(tm.ocean, FT)
+
+function vars_state_auxiliary(m::FlowIntegralModel, T)
+    @vars begin
+        ∫u::SVector{2, T}
+    end
+end
+
+init_state_auxiliary!(tm::FlowIntegralModel, A::Vars, geom::LocalGeometry) =
+    nothing
+
+function vars_integrals(m::FlowIntegralModel, T)
+    @vars begin
+        ∫u::SVector{2, T}
+    end
+end
+
+@inline function integral_load_auxiliary_state!(
+    m::FlowIntegralModel,
+    I::Vars,
+    Q::Vars,
+    A::Vars,
+)
+    I.∫u = Q.u
+
+    return nothing
+end
+
+@inline function integral_set_auxiliary_state!(
+    m::FlowIntegralModel,
+    A::Vars,
+    I::Vars,
+)
+    A.∫u = I.∫u
+
+    return nothing
+end
+
+function update_auxiliary_state!(
+    dg::DGModel,
+    fm::FlowIntegralModel,
+    Q::MPIStateArray,
+    t::Real,
+)
+    A = dg.state_auxiliary
+
+    # compute vertical integral of u
+    indefinite_stack_integral!(dg, fm, Q, A, t) # bottom -> top
+
+    return true
+end
diff --git a/src/Ocean/SplitExplicit/SplitExplicitModel.jl b/src/Ocean/SplitExplicit/SplitExplicitModel.jl
index 3059eb05075..b58c18bdc25 100644
--- a/src/Ocean/SplitExplicit/SplitExplicitModel.jl
+++ b/src/Ocean/SplitExplicit/SplitExplicitModel.jl
@@ -1,144 +1,33 @@
 module SplitExplicit
 
+using StaticArrays
+
 using ..HydrostaticBoussinesq
 using ..ShallowWater
-
-using StaticArrays
+using ..VariableTemplates
+using ..MPIStateArrays
+using ..DGMethods:
+    BalanceLaw,
+    LocalGeometry,
+    DGModel,
+    nodal_update_auxiliary_state!,
+    indefinite_stack_integral!,
+    basic_grid_info
 
 import ..DGMethods:
-    initialize_states!,
-    tendency_from_slow_to_fast!,
-    cummulate_fast_solution!,
-    reconcile_from_fast_to_slow!
-
-@inline function initialize_states!(
-    baroclinic::HydrostaticBoussinesqModel,
-    barotropic::ShallowWaterModel,
-    model_bc,
-    model_bt,
-    state_bc,
-    state_bt,
-)
-    model_bt.state_auxiliary.ηᶜ .= -0
-    model_bt.state_auxiliary.Uᶜ .= (@SVector [-0, -0, -0])'
-
-    state_bt.η .= model_bt.state_auxiliary.ηˢ
-    state_bt.U .= model_bt.state_auxiliary.Uˢ
-
-    model_bc.state_auxiliary.ΔGᵘ .= -0
-
-    return nothing
-end
-
-@inline function tendency_from_slow_to_fast!(
-    baroclinic::HydrostaticBoussinesqModel,
-    barotropic::ShallowWaterModel,
-    model_bc,
-    model_bt,
-    state_bc,
-    state_bt,
-    forcing_tendency,
-)
-    # integrate the tendency
-    tendency_model = model_bc.modeldata.tendency_model
-    update_auxiliary_state!(
-        tendency_model,
-        tendency_model.balance_law,
-        forcing_tendency,
-        0,
-    )
-
-    Nq, Nqk, _, _, nelemv, _, nelemh, _ = basic_grid_info(model_bc)
-
-    ## get top value (=integral over full depth) of ∫du
-    ∫du = tendency_model.state_auxiliary.∫du
-    boxy_∫du = reshape(∫du, Nq^2, Nqk, 2, nelemv, nelemh)
-    flat_∫du = @view boxy_∫du[:, end, :, end, :]
-
-    ## copy into Gᵁ of dgFast
-    model_bt.state_auxiliary.Gᵁ .= reshape(flat_∫du, Nq^2, 2, nelemh)
-
-    ## scale by -1/H and copy back to ΔGu
-    # note: since tendency_dg.auxstate.∫du is not used after this, could be
-    # re-used to store a 3-D copy of "-Gu"
-    ΔGᵘ = model_bc.state_auxiliary.ΔGᵘ
-    boxy_ΔGᵘ = reshape(ΔGᵘ, Nq^2, Nqk, 2, nelemv, nelemh)
-    boxy_ΔGᵘ .= -reshape(flat_∫du, Nq^2, 1, 2, 1, nelemh) / baroclinic.problem.H
-
-    return nothing
-end
-
-@inline function cummulate_fast_solution!(
-    baroclinic::HydrostaticBoussinesqModel,
-    barotropic::ShallowWaterModel,
-    model_bt,
-    state_bt,
-    fast_time,
-    fast_dt,
-    substep,
-)
-    ### for now no averaging, just save the final time step
-    model_bt.state_auxiliary.ηᶜ .= state_bt.η
-    model_bt.state_auxiliary.Uᶜ .= state_bt.U
-
-    ### technically unnecessary right now, but for future changes
-    model_bt.state_auxiliary.ηˢ .= state_bt.η
-    model_bt.state_auxiliary.Uˢ .= state_bt.U
-
-    return nothing
-end
-
-@inline function reconcile_from_fast_to_slow!(
-    baroclinic::HydrostaticBoussinesqModel,
-    barotropic::ShallowWaterModel,
-    model_bc,
-    model_bt,
-    state_bc,
-    state_bt,
-)
-    Nq, Nqk, _, _, nelemv, _, nelemh, _ = basic_grid_info(model_bc)
-
-    # need to calculate int_u using integral kernels
-    # u_slow := u_slow + (1/H) * (u_fast - \int_{-H}^{0} u_slow)
-
-    # Compute: \int_{-H}^{0} u_slow)
-    ### need to make sure this is stored into aux.∫u
-
-    # integrate vertically horizontal velocity
-    flow_model = model_bc.modeldata.flow_model
-    update_auxiliary_state!(flow_model, flow_model.balance_law, state_bc, 0)
-
-    ## get top value (=integral over full depth)
-    ∫u = flow_model.state_auxiliary.∫u
-    boxy_∫u = reshape(∫u, Nq^2, Nqk, 2, nelemv, nelemh)
-    flat_∫u = @view boxy_∫u[:, end, :, end, :]
-
-    ## get time weighted averaged out of cumulative arrays
-    ## no weights right now
-    state_bt.state_auxiliary.Uᶜ .*= 1 # / fast_time_rec[2]
-    state_bt.state_auxiliary.ηᶜ .*= 1 # / fast_time_rec[2]
-
-    ### substract ∫u from U and divide by H
-
-    ### Δu is a place holder for 1/H * (Ū - ∫u)
-    Δu = model_bt.state_auxiliary.Δu
-    Δu .= 1 / baroclinic.problem.H * (model_bt.state_auxiliary.Uᶜ - flat_∫u)
-
-    ### copy the 2D contribution down the 3D solution
-    ### need to reshape these things for the broadcast
-    boxy_u = reshape(state_bc.u, Nq^2, Nqk, 2, nelemv, nelemh)
-    boxy_Δu = reshape(Δu, Nq^2, 1, 3, 1, nelemh)
-    ### this works, we tested it
-    # currently SWModel has a 3-vector for U
-    # need to shrink to 2-vector at some point
-    boxy_u .+= boxy_Δu[:, :, 1:2, :, :]
-
-    ### copy 2D eta over to 3D model
-    boxy_η_3D = reshape(state_bc.η, Nq^2, Nqk, nelemv, nelemh)
-    boxy_η_2D = reshape(model_bt.state_auxiliary.ηᶜ, Nq^2, 1, 1, nelemh)
-    boxy_η_3D .= boxy_η_2D
-
-    return nothing
-end
+    vars_state_conservative,
+    init_state_conservative!,
+    vars_state_auxiliary,
+    init_state_auxiliary!,
+    vars_state_gradient,
+    vars_state_gradient_flux,
+    vars_integrals,
+    integral_load_auxiliary_state!,
+    integral_set_auxiliary_state!,
+    update_auxiliary_state!
+
+include("Communication.jl")
+include("TendencyIntegralModel.jl")
+include("FlowIntegralModel.jl")
 
 end
diff --git a/src/Ocean/SplitExplicit/TendencyIntegralModel.jl b/src/Ocean/SplitExplicit/TendencyIntegralModel.jl
new file mode 100644
index 00000000000..88b1087b1dc
--- /dev/null
+++ b/src/Ocean/SplitExplicit/TendencyIntegralModel.jl
@@ -0,0 +1,72 @@
+struct TendencyIntegralModel{M} <: BalanceLaw
+    ocean::M
+    function TendencyIntegralModel(ocean::M) where {M}
+        return new{M}(ocean)
+    end
+end
+
+vars_state_gradient(tm::TendencyIntegralModel, FT) = @vars()
+vars_state_gradient_flux(tm::TendencyIntegralModel, FT) = @vars()
+
+vars_state_conservative(tm::TendencyIntegralModel, FT) =
+    vars_state_conservative(tm.ocean, FT)
+
+function vars_state_auxiliary(m::TendencyIntegralModel, T)
+    @vars begin
+        ∫du::SVector{2, T}
+    end
+end
+
+init_state_auxiliary!(tm::TendencyIntegralModel, A::Vars, geom::LocalGeometry) =
+    nothing
+
+function vars_integrals(m::TendencyIntegralModel, T)
+    @vars begin
+        ∫du::SVector{2, T}
+    end
+end
+
+@inline function integral_load_auxiliary_state!(
+    m::TendencyIntegralModel,
+    I::Vars,
+    Q::Vars,
+    A::Vars,
+)
+    I.∫du = A.∫du
+
+    return nothing
+end
+
+@inline function integral_set_auxiliary_state!(
+    m::TendencyIntegralModel,
+    A::Vars,
+    I::Vars,
+)
+    A.∫du = I.∫du
+
+    return nothing
+end
+
+function update_auxiliary_state!(
+    dg::DGModel,
+    tm::TendencyIntegralModel,
+    dQ::MPIStateArray,
+    t::Real,
+)
+    A = dg.state_auxiliary
+
+    # copy tendency vector to aux state for integration
+    function f!(::TendencyIntegralModel, dQ, A, t)
+        @inbounds begin
+            A.∫du = @SVector [dQ.u[1], dQ.u[2]]
+        end
+
+        return nothing
+    end
+    nodal_update_auxiliary_state!(f!, dg, tm, dQ, t)
+
+    # compute integral for Gᵁ
+    indefinite_stack_integral!(dg, tm, dQ, A, t) # bottom -> top
+
+    return true
+end
diff --git a/test/Ocean/SplitExplicit/hydrostatic_spindown.jl b/test/Ocean/SplitExplicit/hydrostatic_spindown.jl
index cbd36e41d11..a33f29c9722 100644
--- a/test/Ocean/SplitExplicit/hydrostatic_spindown.jl
+++ b/test/Ocean/SplitExplicit/hydrostatic_spindown.jl
@@ -9,6 +9,7 @@ using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
 using ClimateMachine.HydrostaticBoussinesq
 using ClimateMachine.ShallowWater
+using ClimateMachine.SplitExplicit: FlowIntegralModel, TendencyIntegralModel
 
 using ClimateMachine.Mesh.Topologies
 using ClimateMachine.Mesh.Grids
@@ -148,7 +149,29 @@ function run_hydrostatic_spindown(; coupled = true, refDat = ())
 
     vert_filter = CutoffFilter(grid_3D, polynomialorder(grid_3D) - 1)
     exp_filter = ExponentialFilter(grid_3D, 1, 8)
-    modeldata = (vert_filter = vert_filter, exp_filter = exp_filter)
+
+    flow_model = DGModel(
+        FlowIntegralModel(model_3D),
+        grid_3D,
+        RusanovNumericalFlux(),
+        CentralNumericalFluxSecondOrder(),
+        CentralNumericalFluxGradient(),
+    )
+
+    tendency_model = DGModel(
+        TendencyIntegralModel(model_3D),
+        grid_3D,
+        RusanovNumericalFlux(),
+        CentralNumericalFluxSecondOrder(),
+        CentralNumericalFluxGradient(),
+    )
+
+    modeldata = (
+        vert_filter = vert_filter,
+        exp_filter = exp_filter,
+        flow_model = flow_model,
+        tendency_model = tendency_model,
+    )
 
     dg_3D = DGModel(
         model_3D,
@@ -333,4 +356,5 @@ const cᶻ = 0
     include("../refvals/hydrostatic_spindown_refvals.jl")
 
     run_hydrostatic_spindown(coupled = false, refDat = refVals.uncoupled)
+    # run_hydrostatic_spindown(coupled = true)
 end