diff --git a/docs/src/APIs/BalanceLaws/BalanceLaws.md b/docs/src/APIs/BalanceLaws/BalanceLaws.md
index e4fd778d231..127bb0d3671 100644
--- a/docs/src/APIs/BalanceLaws/BalanceLaws.md
+++ b/docs/src/APIs/BalanceLaws/BalanceLaws.md
@@ -60,6 +60,7 @@ transform_post_gradient_laplacian!
 ```@docs
 wavespeed
 boundary_state!
+nodal_update_auxiliary_state!
 update_auxiliary_state!
 update_auxiliary_state_gradient!
 ```
diff --git a/experiments/OceanBoxGCM/homogeneous_box.jl b/experiments/OceanBoxGCM/homogeneous_box.jl
index 494832830e7..47835c6a199 100755
--- a/experiments/OceanBoxGCM/homogeneous_box.jl
+++ b/experiments/OceanBoxGCM/homogeneous_box.jl
@@ -8,7 +8,7 @@ using ClimateMachine.Mesh.Filters
 using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
 using ClimateMachine.BalanceLaws: vars_state_conservative
-using ClimateMachine.HydrostaticBoussinesq
+using ClimateMachine.Ocean.HydrostaticBoussinesq
 
 using Test
 using CLIMAParameters
diff --git a/experiments/OceanBoxGCM/ocean_gyre.jl b/experiments/OceanBoxGCM/ocean_gyre.jl
index ae3cd1f125d..bb2b359b1ad 100755
--- a/experiments/OceanBoxGCM/ocean_gyre.jl
+++ b/experiments/OceanBoxGCM/ocean_gyre.jl
@@ -7,7 +7,7 @@ using ClimateMachine.ODESolvers
 using ClimateMachine.Mesh.Filters
 using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
-using ClimateMachine.HydrostaticBoussinesq
+using ClimateMachine.Ocean.HydrostaticBoussinesq
 
 using Test
 
diff --git a/src/BalanceLaws/BalanceLaws.jl b/src/BalanceLaws/BalanceLaws.jl
index 6d57fbc9b6b..9a9ec4d068e 100644
--- a/src/BalanceLaws/BalanceLaws.jl
+++ b/src/BalanceLaws/BalanceLaws.jl
@@ -18,19 +18,22 @@ export BalanceLaw,
     vars_state_gradient_flux,
     init_state_conservative!,
     init_state_auxiliary!,
-    flux_first_order!,
-    flux_second_order!,
     compute_gradient_flux!,
     compute_gradient_argument!,
-    source!,
     transform_post_gradient_laplacian!,
+    flux_first_order!,
+    flux_second_order!,
+    source!,
     wavespeed,
     boundary_state!,
+    nodal_update_auxiliary_state!,
     update_auxiliary_state!,
     update_auxiliary_state_gradient!,
+    vars_integrals,
     integral_load_auxiliary_state!,
     integral_set_auxiliary_state!,
     indefinite_stack_integral!,
+    vars_reverse_integrals,
     reverse_indefinite_stack_integral!,
     reverse_integral_load_auxiliary_state!,
     reverse_integral_set_auxiliary_state!
diff --git a/src/ClimateMachine.jl b/src/ClimateMachine.jl
index 0a858964cb5..0b517da61b5 100644
--- a/src/ClimateMachine.jl
+++ b/src/ClimateMachine.jl
@@ -26,6 +26,8 @@ include(joinpath("BalanceLaws", "BalanceLaws.jl"))
 include(joinpath("Numerics", "DGMethods", "DGMethods.jl"))
 include(joinpath("Common", "Orientations", "Orientations.jl"))
 include(joinpath("Utilities", "SingleStackUtils", "SingleStackUtils.jl"))
+include(joinpath("Ocean/Ocean.jl"))
+#=
 include(joinpath("Ocean", "ShallowWater", "ShallowWaterModel.jl"))
 include(joinpath(
     "Ocean",
@@ -33,6 +35,7 @@ include(joinpath(
     "HydrostaticBoussinesqModel.jl",
 ))
 include(joinpath("Ocean", "SplitExplicit", "SplitExplicitModel.jl"))
+=#
 include(joinpath("Numerics", "SystemSolvers", "SystemSolvers.jl"))
 include(joinpath("Numerics", "ODESolvers", "GenericCallbacks.jl"))
 include(joinpath("Numerics", "ODESolvers", "ODESolvers.jl"))
diff --git a/src/Driver/Driver.jl b/src/Driver/Driver.jl
index c4713e7dab3..019c66f675a 100644
--- a/src/Driver/Driver.jl
+++ b/src/Driver/Driver.jl
@@ -24,7 +24,7 @@ using ..BalanceLaws:
     update_auxiliary_state!
 using ..DGMethods: remainder_DGModel
 using ..DGMethods.NumericalFluxes
-using ..HydrostaticBoussinesq
+using ..Ocean.HydrostaticBoussinesq
 using ..Mesh.Grids
 using ..Mesh.Topologies
 using ..Thermodynamics
diff --git a/src/Ocean/HydrostaticBoussinesq/BoundaryConditions.jl b/src/Ocean/HydrostaticBoussinesq/BoundaryConditions.jl
index ca22ab71e2d..8e255f2de12 100644
--- a/src/Ocean/HydrostaticBoussinesq/BoundaryConditions.jl
+++ b/src/Ocean/HydrostaticBoussinesq/BoundaryConditions.jl
@@ -1,7 +1,6 @@
-using ..DGMethods.NumericalFluxes:
-    NumericalFluxFirstOrder, NumericalFluxGradient, NumericalFluxSecondOrder
+using ...DGMethods.NumericalFluxes
 
-import ..BalanceLaws: boundary_state!
+import ...BalanceLaws: boundary_state!
 
 export OceanBC
 
diff --git a/src/Ocean/HydrostaticBoussinesq/Courant.jl b/src/Ocean/HydrostaticBoussinesq/Courant.jl
index 70d758d0767..cbb8f551ab2 100644
--- a/src/Ocean/HydrostaticBoussinesq/Courant.jl
+++ b/src/Ocean/HydrostaticBoussinesq/Courant.jl
@@ -1,14 +1,14 @@
 using Logging, Printf
 using LinearAlgebra: norm
 
-using ..Mesh.Grids:
+using ...Mesh.Grids:
     VerticalDirection, HorizontalDirection, EveryDirection, min_node_distance
-using ..DGMethods: courant
+using ...DGMethods: courant
 
-import ..Courant:
+import ...Courant:
     advective_courant, nondiffusive_courant, diffusive_courant, viscous_courant
 
-import ..DGMethods: calculate_dt
+import ...DGMethods: calculate_dt
 
 """
     advective_courant(::HBModel)
diff --git a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
index b0e896643e3..044c56c558b 100644
--- a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
+++ b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
@@ -6,50 +6,44 @@ using StaticArrays
 using LinearAlgebra: dot, Diagonal
 using CLIMAParameters.Planet: grav
 
-using ..VariableTemplates
-using ..MPIStateArrays
-using ..Mesh.Filters: apply!
-using ..Mesh.Grids: VerticalDirection
-using ..BalanceLaws: BalanceLaw
-import ..BalanceLaws: nodal_update_auxiliary_state!
-using ..DGMethods.NumericalFluxes: RusanovNumericalFlux
-
-import ..DGMethods.NumericalFluxes: update_penalty!
-import ..DGMethods:
+using ...VariableTemplates
+using ...MPIStateArrays
+using ...Mesh.Filters: apply!
+using ...Mesh.Grids: VerticalDirection
+using ...Mesh.Geometry
+using ...DGMethods: DGModel, copy_stack_field_down!
+using ...DGMethods.NumericalFluxes
+using ...BalanceLaws
+
+import ...DGMethods.NumericalFluxes: update_penalty!
+import ...BalanceLaws:
     vars_state_conservative,
     vars_state_auxiliary,
     vars_state_gradient,
-    vars_integrals,
     vars_state_gradient_flux,
-    vars_reverse_integrals,
     init_state_conservative!,
     init_state_auxiliary!,
     compute_gradient_argument!,
     compute_gradient_flux!,
-    integral_load_auxiliary_state!,
-    integral_set_auxiliary_state!,
-    indefinite_stack_integral!,
-    reverse_indefinite_stack_integral!,
-    indefinite_stack_integral!,
-    reverse_indefinite_stack_integral!,
-    reverse_integral_load_auxiliary_state!,
-    reverse_integral_set_auxiliary_state!,
     flux_first_order!,
     flux_second_order!,
     source!,
     wavespeed,
+    boundary_state!,
     update_auxiliary_state!,
-    update_auxiliary_state_gradient!
-
-using ..DGMethods: DGModel, copy_stack_field_down!
-using ..Mesh.Geometry: LocalGeometry
-
-using ..DGMethods.NumericalFluxes: RusanovNumericalFlux
-
-import ..DGMethods.NumericalFluxes: update_penalty!
+    update_auxiliary_state_gradient!,
+    vars_integrals,
+    integral_load_auxiliary_state!,
+    integral_set_auxiliary_state!,
+    indefinite_stack_integral!,
+    vars_reverse_integrals,
+    reverse_indefinite_stack_integral!,
+    reverse_integral_load_auxiliary_state!,
+    reverse_integral_set_auxiliary_state!
 
 ×(a::SVector, b::SVector) = StaticArrays.cross(a, b)
 ⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b)
+⊗(a::SVector, b::SVector) = a * b'
 
 abstract type AbstractHydrostaticBoussinesqProblem end
 
@@ -296,7 +290,7 @@ location to store integrands for bottom up integrals
 """
 function vars_integrals(m::HBModel, T)
     @vars begin
-        ∇hu::T
+        ∇ʰu::T
         αᵀθ::T
     end
 end
@@ -319,7 +313,7 @@ A -> array of aux variables
     Q::Vars,
     A::Vars,
 )
-    I.∇hu = A.w # borrow the w value from A...
+    I.∇ʰu = A.w # borrow the w value from A...
     I.αᵀθ = -m.αᵀ * Q.θ # integral will be reversed below
 
     return nothing
@@ -337,7 +331,7 @@ A -> array of aux variables
 I -> array of integrand variables
 """
 @inline function integral_set_auxiliary_state!(m::HBModel, A::Vars, I::Vars)
-    A.w = I.∇hu
+    A.w = I.∇ʰu
     A.pkin = I.αᵀθ
 
     return nothing
diff --git a/src/Ocean/Ocean.jl b/src/Ocean/Ocean.jl
new file mode 100644
index 00000000000..9dba464b5f5
--- /dev/null
+++ b/src/Ocean/Ocean.jl
@@ -0,0 +1,8 @@
+module Ocean
+
+include("HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl")
+include("ShallowWater/ShallowWaterModel.jl")
+include("SplitExplicit/SplitExplicitModel.jl")
+# include("OceanProblems/SimpleBoxProblem.jl")
+
+end
diff --git a/src/Ocean/ShallowWater/ShallowWaterModel.jl b/src/Ocean/ShallowWater/ShallowWaterModel.jl
index a45efb93d09..d207fb3609d 100644
--- a/src/Ocean/ShallowWater/ShallowWaterModel.jl
+++ b/src/Ocean/ShallowWater/ShallowWaterModel.jl
@@ -3,30 +3,30 @@ module ShallowWater
 export ShallowWaterModel, ShallowWaterProblem
 
 using StaticArrays
-using ..VariableTemplates
-using LinearAlgebra: Diagonal, dot
+using LinearAlgebra: dot, Diagonal
 using CLIMAParameters.Planet: grav
 
-import ..BalanceLaws:
-    BalanceLaw,
-    vars_state_auxiliary,
+using ...VariableTemplates
+using ...Mesh.Geometry
+using ...DGMethods
+using ...DGMethods.NumericalFluxes
+using ...BalanceLaws
+
+import ...DGMethods.NumericalFluxes: update_penalty!
+import ...BalanceLaws:
     vars_state_conservative,
+    vars_state_auxiliary,
     vars_state_gradient,
     vars_state_gradient_flux,
-    vars_integrals,
+    init_state_conservative!,
+    init_state_auxiliary!,
+    compute_gradient_argument!,
+    compute_gradient_flux!,
     flux_first_order!,
     flux_second_order!,
     source!,
     wavespeed,
-    boundary_state!,
-    compute_gradient_argument!,
-    init_state_auxiliary!,
-    init_state_conservative!,
-    compute_gradient_flux!
-
-import ...Mesh.Geometry: LocalGeometry
-
-using ..DGMethods.NumericalFluxes
+    boundary_state!
 
 ×(a::SVector, b::SVector) = StaticArrays.cross(a, b)
 ⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b)
diff --git a/src/Ocean/SplitExplicit/SplitExplicitModel.jl b/src/Ocean/SplitExplicit/SplitExplicitModel.jl
index 334a8e597be..153344ddd39 100644
--- a/src/Ocean/SplitExplicit/SplitExplicitModel.jl
+++ b/src/Ocean/SplitExplicit/SplitExplicitModel.jl
@@ -3,7 +3,7 @@ module SplitExplicit
 using ..HydrostaticBoussinesq
 using ..ShallowWater
 
-import ..BalanceLaws:
+import ...BalanceLaws:
     initialize_states!,
     tendency_from_slow_to_fast!,
     cummulate_fast_solution!,
diff --git a/test/Ocean/HydrostaticBoussinesq/3D_hydrostatic_spindown.jl b/test/Ocean/HydrostaticBoussinesq/3D_hydrostatic_spindown.jl
index e0b432404fe..323c9dcc30a 100644
--- a/test/Ocean/HydrostaticBoussinesq/3D_hydrostatic_spindown.jl
+++ b/test/Ocean/HydrostaticBoussinesq/3D_hydrostatic_spindown.jl
@@ -8,7 +8,7 @@ using ClimateMachine.ODESolvers
 using ClimateMachine.Mesh.Filters
 using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
-using ClimateMachine.HydrostaticBoussinesq
+using ClimateMachine.Ocean.HydrostaticBoussinesq
 
 using Test
 
diff --git a/test/Ocean/HydrostaticBoussinesq/test_divergence_free.jl b/test/Ocean/HydrostaticBoussinesq/test_divergence_free.jl
index a44d7b47cac..c42c1f3660e 100644
--- a/test/Ocean/HydrostaticBoussinesq/test_divergence_free.jl
+++ b/test/Ocean/HydrostaticBoussinesq/test_divergence_free.jl
@@ -7,7 +7,7 @@ using ClimateMachine.Mesh.Filters
 using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
 using ClimateMachine.BalanceLaws: vars_state_conservative
-using ClimateMachine.HydrostaticBoussinesq
+using ClimateMachine.Ocean.HydrostaticBoussinesq
 
 using CLIMAParameters
 using CLIMAParameters.Planet: grav
diff --git a/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl b/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl
index f0554c25fc0..91d6a0ca8aa 100644
--- a/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl
+++ b/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl
@@ -6,7 +6,7 @@ using ClimateMachine.ODESolvers
 using ClimateMachine.Mesh.Filters
 using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
-using ClimateMachine.HydrostaticBoussinesq
+using ClimateMachine.Ocean.HydrostaticBoussinesq
 
 using CLIMAParameters
 using CLIMAParameters.Planet: grav
diff --git a/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl b/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl
index 5902b9954f7..8a96af6f1a3 100644
--- a/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl
+++ b/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl
@@ -7,8 +7,7 @@ using ClimateMachine.ODESolvers
 using ClimateMachine.Mesh.Filters
 using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
-using ClimateMachine.HydrostaticBoussinesq
-using ClimateMachine.ShallowWater
+using ClimateMachine.Ocean.ShallowWater
 
 using ClimateMachine.Mesh.Topologies
 using ClimateMachine.Mesh.Grids
@@ -23,7 +22,7 @@ using LinearAlgebra
 using StaticArrays
 using Logging, Printf, Dates
 
-import ClimateMachine.ShallowWater: shallow_init_state!, shallow_init_aux!
+import ClimateMachine.Ocean.ShallowWater: shallow_init_state!, shallow_init_aux!
 
 using CLIMAParameters
 using CLIMAParameters.Planet: grav
diff --git a/test/Ocean/ShallowWater/GyreInABox.jl b/test/Ocean/ShallowWater/GyreInABox.jl
index 10cdcf4ca8d..fc3716922dc 100644
--- a/test/Ocean/ShallowWater/GyreInABox.jl
+++ b/test/Ocean/ShallowWater/GyreInABox.jl
@@ -8,12 +8,12 @@ using ClimateMachine.MPIStateArrays
 using ClimateMachine.ODESolvers
 using ClimateMachine.GenericCallbacks
 using ClimateMachine.VariableTemplates: flattenednames
-using ClimateMachine.ShallowWater
+using ClimateMachine.Ocean.ShallowWater
 using LinearAlgebra
 using StaticArrays
 using Logging, Printf, Dates
 using ClimateMachine.VTK
-import ClimateMachine.ShallowWater:
+import ClimateMachine.Ocean.ShallowWater:
     shallow_init_state!,
     shallow_init_aux!,
     vars_state_conservative,
diff --git a/test/Ocean/SplitExplicit/hydrostatic_spindown.jl b/test/Ocean/SplitExplicit/hydrostatic_spindown.jl
index 34a38342aa1..aaaf1432d78 100644
--- a/test/Ocean/SplitExplicit/hydrostatic_spindown.jl
+++ b/test/Ocean/SplitExplicit/hydrostatic_spindown.jl
@@ -7,8 +7,8 @@ using ClimateMachine.ODESolvers
 using ClimateMachine.Mesh.Filters
 using ClimateMachine.VariableTemplates
 using ClimateMachine.Mesh.Grids: polynomialorder
-using ClimateMachine.HydrostaticBoussinesq
-using ClimateMachine.ShallowWater
+using ClimateMachine.Ocean.HydrostaticBoussinesq
+using ClimateMachine.Ocean.ShallowWater
 
 using ClimateMachine.Mesh.Topologies
 using ClimateMachine.Mesh.Grids
@@ -23,7 +23,7 @@ using LinearAlgebra
 using StaticArrays
 using Logging, Printf, Dates
 
-import ClimateMachine.ShallowWater: shallow_init_state!, shallow_init_aux!
+import ClimateMachine.Ocean.ShallowWater: shallow_init_state!, shallow_init_aux!
 
 using CLIMAParameters
 using CLIMAParameters.Planet: grav