From 2d8308edcdf8eda4cbfc13db93fd3ae64d984e02 Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Tue, 29 Jul 2025 11:55:26 -0400
Subject: [PATCH 01/17] starting to add the complete Coriolis for the
 nonhydrostatic model on a sphere

---
 .../nonhydrostatic_spherical_coriolis.jl      | 102 ++++++++++++++++++
 1 file changed, 102 insertions(+)
 create mode 100644 src/Coriolis/nonhydrostatic_spherical_coriolis.jl

diff --git a/src/Coriolis/nonhydrostatic_spherical_coriolis.jl b/src/Coriolis/nonhydrostatic_spherical_coriolis.jl
new file mode 100644
index 0000000000..8b51413e97
--- /dev/null
+++ b/src/Coriolis/nonhydrostatic_spherical_coriolis.jl
@@ -0,0 +1,102 @@
+using Oceananigans.Grids: LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, peripheral_node, φnode
+using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind, hack_cosd
+using Oceananigans.Advection: EnergyConserving, EnstrophyConserving
+using Oceananigans.BoundaryConditions
+using Oceananigans.Fields
+using Oceananigans.AbstractOperations: KernelFunctionOperation
+using Oceananigans.ImmersedBoundaries
+
+"""
+    struct NonhydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
+
+A parameter object for constant rotation around a vertical axis on the sphere.
+"""
+struct NonhydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
+    rotation_rate :: FT
+    scheme :: S
+end
+
+"""
+    NonhydrostaticSphericalCoriolis([FT=Float64;]
+                                 rotation_rate = Ω_Earth,
+                                 scheme = EnstrophyConserving())
+
+Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`.
+
+Keyword arguments
+=================
+
+- `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref).
+- `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `EnstrophyConserving()` (default).
+"""
+function NonhydrostaticSphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
+                                      rotation_rate = Ω_Earth,
+                                      scheme :: S = EnstrophyConserving(FT)) where S
+
+    return HydrostaticSphericalCoriolis{S, FT}(rotation_rate, scheme)
+end
+
+Adapt.adapt_structure(to, coriolis::NonhydrostaticSphericalCoriolis) =
+    NonhydrostaticSphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
+                                    Adapt.adapt(to, coriolis.scheme))
+
+@inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
+@inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
+@inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)
+
+@inline fᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis) =
+    2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
+
+@inline fntᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis) =
+    2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid))
+
+#FJP: make this non-zero
+@inline z_f_cross_U(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis, U) = zero(grid)
+
+#####
+##### Active Point Enstrophy-conserving scheme
+#####
+
+# It might happen that a cell is active but all the neighbouring staggered nodes are inactive,
+# (an example is a 1-cell large channel)
+# In that case the Coriolis force is equal to zero
+
+const CoriolisEnstrophyConserving = NonhydrostaticSphericalCoriolis{<:EnstrophyConserving}
+
+@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
+    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
+                active_weighted_ℑxyᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+
+@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
+    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
+                active_weighted_ℑxyᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+
+#####
+##### Energy-conserving scheme
+#####
+
+const CoriolisEnergyConserving = NonhydrostaticSphericalCoriolis{<:EnergyConserving}
+
+@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
+@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
+
+@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
+    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+
+@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
+    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+
+#####
+##### Show
+#####
+
+function Base.show(io::IO, nonhydrostatic_spherical_coriolis::NonhydrostaticSphericalCoriolis)
+    coriolis_scheme = nonhydrostatic_spherical_coriolis.scheme
+    rotation_rate   = nonhydrostatic_spherical_coriolis.rotation_rate
+    rotation_rate_Earth = rotation_rate / Ω_Earth
+    rotation_rate_str = @sprintf("%.2e s⁻¹ = %.2e Ω_Earth", rotation_rate, rotation_rate_Earth)
+
+    return print(io, "NonhydrostaticSphericalCoriolis", '\n',
+                 "├─ rotation rate: ", rotation_rate_str, '\n',
+                 "└─ scheme: ", summary(coriolis_scheme))
+end

From 5d958326ca8db40a23f527cd90afbe7763f65c9d Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Tue, 29 Jul 2025 16:51:01 -0400
Subject: [PATCH 02/17] First attempt to include the complete Coriolis force

---
 .../nonhydrostatic_spherical_coriolis.jl      | 73 ++++++-------------
 1 file changed, 21 insertions(+), 52 deletions(-)

diff --git a/src/Coriolis/nonhydrostatic_spherical_coriolis.jl b/src/Coriolis/nonhydrostatic_spherical_coriolis.jl
index 8b51413e97..b0a7fc2e51 100644
--- a/src/Coriolis/nonhydrostatic_spherical_coriolis.jl
+++ b/src/Coriolis/nonhydrostatic_spherical_coriolis.jl
@@ -10,16 +10,18 @@ using Oceananigans.ImmersedBoundaries
     struct NonhydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
 
 A parameter object for constant rotation around a vertical axis on the sphere.
+
+This is only designed to work for a LatitudeLongitudeGrid.
+
+There is no attempt to conserve energy or enstrophy at this early stage.
 """
 struct NonhydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
     rotation_rate :: FT
-    scheme :: S
 end
 
 """
     NonhydrostaticSphericalCoriolis([FT=Float64;]
                                  rotation_rate = Ω_Earth,
-                                 scheme = EnstrophyConserving())
 
 Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`.
 
@@ -27,76 +29,43 @@ Keyword arguments
 =================
 
 - `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref).
-- `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `EnstrophyConserving()` (default).
 """
 function NonhydrostaticSphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
                                       rotation_rate = Ω_Earth,
-                                      scheme :: S = EnstrophyConserving(FT)) where S
 
-    return HydrostaticSphericalCoriolis{S, FT}(rotation_rate, scheme)
+    return NonhydrostaticSphericalCoriolis{FT}(rotation_rate)
 end
 
 Adapt.adapt_structure(to, coriolis::NonhydrostaticSphericalCoriolis) =
-    NonhydrostaticSphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
-                                    Adapt.adapt(to, coriolis.scheme))
+    NonhydrostaticSphericalCoriolis(Adapt.adapt(to))
 
 @inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
-@inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
-@inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)
-
-@inline fᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis) =
-    2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
-
-@inline fntᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis) =
-    2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid))
-
-#FJP: make this non-zero
-@inline z_f_cross_U(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis, U) = zero(grid)
-
-#####
-##### Active Point Enstrophy-conserving scheme
-#####
-
-# It might happen that a cell is active but all the neighbouring staggered nodes are inactive,
-# (an example is a 1-cell large channel)
-# In that case the Coriolis force is equal to zero
 
-const CoriolisEnstrophyConserving = NonhydrostaticSphericalCoriolis{<:EnstrophyConserving}
-
-@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
-    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
-                active_weighted_ℑxyᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
-
-@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
-    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
-                active_weighted_ℑxyᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
-
-#####
-##### Energy-conserving scheme
-#####
-
-const CoriolisEnergyConserving = NonhydrostaticSphericalCoriolis{<:EnergyConserving}
+@inline fᶠᶠᵃ(i, j, k, grid) = 2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
+@inline f̃ᶠᶠᵃ(i, j, k, grid) = 2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid))
 
 @inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
+@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis, v) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
+
 @inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
 
-@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
-    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u)
+
+@inline x_f_cross_U(i, j, k, grid, U) =
+    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)  
+              + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
 
-@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
-    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+@inline y_f_cross_U(i, j, k, grid, U) =
+    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
 
-#####
-##### Show
-#####
+@inline z_f_cross_U(i, j, k, grid, U) =
+    @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f̃_ℑz_uᶠᶠᵃ, U[1]) * Δz⁻¹ᶜᶜᶠ(i, j, k, grid)
 
 function Base.show(io::IO, nonhydrostatic_spherical_coriolis::NonhydrostaticSphericalCoriolis)
-    coriolis_scheme = nonhydrostatic_spherical_coriolis.scheme
-    rotation_rate   = nonhydrostatic_spherical_coriolis.rotation_rate
+    rotation_rate   = hydrostatic_spherical_coriolis.rotation_rate
     rotation_rate_Earth = rotation_rate / Ω_Earth
     rotation_rate_str = @sprintf("%.2e s⁻¹ = %.2e Ω_Earth", rotation_rate, rotation_rate_Earth)
 
     return print(io, "NonhydrostaticSphericalCoriolis", '\n',
-                 "├─ rotation rate: ", rotation_rate_str, '\n',
-                 "└─ scheme: ", summary(coriolis_scheme))
+                 "├─ rotation rate: ", rotation_rate_str, '\n')
 end

From 01fdbc8bad4ce9764f244cd0398de74c2dbd486c Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Thu, 31 Jul 2025 09:03:22 -0400
Subject: [PATCH 03/17] Starting to build SphericalCoriolis

---
 src/Coriolis/Coriolis.jl                                      | 4 ++--
 ...drostatic_spherical_coriolis.jl => spherical_corioilis.jl} | 0
 2 files changed, 2 insertions(+), 2 deletions(-)
 rename src/Coriolis/{hydrostatic_spherical_coriolis.jl => spherical_corioilis.jl} (100%)

diff --git a/src/Coriolis/Coriolis.jl b/src/Coriolis/Coriolis.jl
index 651abac6eb..1dd49bbaeb 100644
--- a/src/Coriolis/Coriolis.jl
+++ b/src/Coriolis/Coriolis.jl
@@ -2,7 +2,7 @@ module Coriolis
 
 export
     FPlane, ConstantCartesianCoriolis, BetaPlane, NonTraditionalBetaPlane,
-    HydrostaticSphericalCoriolis, ActiveCellEnstrophyConserving,
+    SphericalCoriolis, ActiveCellEnstrophyConserving,
     x_f_cross_U, y_f_cross_U, z_f_cross_U
 
 using Printf
@@ -33,6 +33,6 @@ include("f_plane.jl")
 include("constant_cartesian_coriolis.jl")
 include("beta_plane.jl")
 include("non_traditional_beta_plane.jl")
-include("hydrostatic_spherical_coriolis.jl")
+include("spherical_coriolis.jl")
 
 end # module
diff --git a/src/Coriolis/hydrostatic_spherical_coriolis.jl b/src/Coriolis/spherical_corioilis.jl
similarity index 100%
rename from src/Coriolis/hydrostatic_spherical_coriolis.jl
rename to src/Coriolis/spherical_corioilis.jl

From 5b275e982d08086cfb4e39d5123b9c5049566b71 Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Sat, 2 Aug 2025 15:34:03 -0400
Subject: [PATCH 04/17] fix silly typo

---
 src/Coriolis/spherical_coriolis.jl | 98 ++++++++++++++++++++++++++++++
 1 file changed, 98 insertions(+)
 create mode 100644 src/Coriolis/spherical_coriolis.jl

diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
new file mode 100644
index 0000000000..2c80696058
--- /dev/null
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -0,0 +1,98 @@
+using Oceananigans.Grids: LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, peripheral_node, φnode
+using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind
+using Oceananigans.Advection: EnergyConserving, EnstrophyConserving
+using Oceananigans.BoundaryConditions
+using Oceananigans.Fields
+using Oceananigans.AbstractOperations: KernelFunctionOperation
+using Oceananigans.ImmersedBoundaries
+
+"""
+    struct HydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
+
+A parameter object for constant rotation around a vertical axis on the sphere.
+"""
+struct HydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
+    rotation_rate :: FT
+    scheme :: S
+end
+
+"""
+    HydrostaticSphericalCoriolis([FT=Float64;]
+                                 rotation_rate = Ω_Earth,
+                                 scheme = EnstrophyConserving())
+
+Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`.
+
+Keyword arguments
+=================
+
+- `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref).
+- `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `EnstrophyConserving()` (default).
+"""
+function HydrostaticSphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
+                                      rotation_rate = Ω_Earth,
+                                      scheme :: S = EnstrophyConserving(FT)) where S
+
+    return HydrostaticSphericalCoriolis{S, FT}(rotation_rate, scheme)
+end
+
+Adapt.adapt_structure(to, coriolis::HydrostaticSphericalCoriolis) =
+    HydrostaticSphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
+                                 Adapt.adapt(to, coriolis.scheme))
+
+@inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
+@inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
+@inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)
+
+@inline fᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis) =
+    2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
+
+@inline z_f_cross_U(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, U) = zero(grid)
+
+#####
+##### Active Point Enstrophy-conserving scheme
+#####
+
+# It might happen that a cell is active but all the neighbouring staggered nodes are inactive,
+# (an example is a 1-cell large channel)
+# In that case the Coriolis force is equal to zero
+
+const CoriolisEnstrophyConserving = HydrostaticSphericalCoriolis{<:EnstrophyConserving}
+
+@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
+    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
+                active_weighted_ℑxyᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+
+@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
+    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
+                active_weighted_ℑxyᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+
+#####
+##### Energy-conserving scheme
+#####
+
+const CoriolisEnergyConserving = HydrostaticSphericalCoriolis{<:EnergyConserving}
+
+@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
+@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
+
+@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
+    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+
+@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
+    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+
+#####
+##### Show
+#####
+
+function Base.show(io::IO, hydrostatic_spherical_coriolis::HydrostaticSphericalCoriolis)
+    coriolis_scheme = hydrostatic_spherical_coriolis.scheme
+    rotation_rate   = hydrostatic_spherical_coriolis.rotation_rate
+    rotation_rate_Earth = rotation_rate / Ω_Earth
+    rotation_rate_str = @sprintf("%.2e s⁻¹ = %.2e Ω_Earth", rotation_rate, rotation_rate_Earth)
+
+    return print(io, "HydrostaticSphericalCoriolis", '\n',
+                 "├─ rotation rate: ", rotation_rate_str, '\n',
+                 "└─ scheme: ", summary(coriolis_scheme))
+end

From d4ad5bcb439b216d182467bd2e70e40dac9d15f7 Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Sat, 2 Aug 2025 16:15:08 -0400
Subject: [PATCH 05/17] first draft of SphericalCoriolis

---
 src/Coriolis/spherical_corioilis.jl |  98 --------------------------
 src/Coriolis/spherical_coriolis.jl  | 103 ++++++++++++++++++----------
 2 files changed, 66 insertions(+), 135 deletions(-)
 delete mode 100644 src/Coriolis/spherical_corioilis.jl

diff --git a/src/Coriolis/spherical_corioilis.jl b/src/Coriolis/spherical_corioilis.jl
deleted file mode 100644
index 2c80696058..0000000000
--- a/src/Coriolis/spherical_corioilis.jl
+++ /dev/null
@@ -1,98 +0,0 @@
-using Oceananigans.Grids: LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, peripheral_node, φnode
-using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind
-using Oceananigans.Advection: EnergyConserving, EnstrophyConserving
-using Oceananigans.BoundaryConditions
-using Oceananigans.Fields
-using Oceananigans.AbstractOperations: KernelFunctionOperation
-using Oceananigans.ImmersedBoundaries
-
-"""
-    struct HydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
-
-A parameter object for constant rotation around a vertical axis on the sphere.
-"""
-struct HydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
-    rotation_rate :: FT
-    scheme :: S
-end
-
-"""
-    HydrostaticSphericalCoriolis([FT=Float64;]
-                                 rotation_rate = Ω_Earth,
-                                 scheme = EnstrophyConserving())
-
-Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`.
-
-Keyword arguments
-=================
-
-- `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref).
-- `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `EnstrophyConserving()` (default).
-"""
-function HydrostaticSphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
-                                      rotation_rate = Ω_Earth,
-                                      scheme :: S = EnstrophyConserving(FT)) where S
-
-    return HydrostaticSphericalCoriolis{S, FT}(rotation_rate, scheme)
-end
-
-Adapt.adapt_structure(to, coriolis::HydrostaticSphericalCoriolis) =
-    HydrostaticSphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
-                                 Adapt.adapt(to, coriolis.scheme))
-
-@inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
-@inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
-@inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)
-
-@inline fᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis) =
-    2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
-
-@inline z_f_cross_U(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, U) = zero(grid)
-
-#####
-##### Active Point Enstrophy-conserving scheme
-#####
-
-# It might happen that a cell is active but all the neighbouring staggered nodes are inactive,
-# (an example is a 1-cell large channel)
-# In that case the Coriolis force is equal to zero
-
-const CoriolisEnstrophyConserving = HydrostaticSphericalCoriolis{<:EnstrophyConserving}
-
-@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
-    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
-                active_weighted_ℑxyᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
-
-@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
-    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
-                active_weighted_ℑxyᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
-
-#####
-##### Energy-conserving scheme
-#####
-
-const CoriolisEnergyConserving = HydrostaticSphericalCoriolis{<:EnergyConserving}
-
-@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
-@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
-
-@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
-    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
-
-@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
-    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
-
-#####
-##### Show
-#####
-
-function Base.show(io::IO, hydrostatic_spherical_coriolis::HydrostaticSphericalCoriolis)
-    coriolis_scheme = hydrostatic_spherical_coriolis.scheme
-    rotation_rate   = hydrostatic_spherical_coriolis.rotation_rate
-    rotation_rate_Earth = rotation_rate / Ω_Earth
-    rotation_rate_str = @sprintf("%.2e s⁻¹ = %.2e Ω_Earth", rotation_rate, rotation_rate_Earth)
-
-    return print(io, "HydrostaticSphericalCoriolis", '\n',
-                 "├─ rotation rate: ", rotation_rate_str, '\n',
-                 "└─ scheme: ", summary(coriolis_scheme))
-end
diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index 2c80696058..398095dffc 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -1,5 +1,5 @@
 using Oceananigans.Grids: LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, peripheral_node, φnode
-using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind
+using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind, hack_cosd
 using Oceananigans.Advection: EnergyConserving, EnstrophyConserving
 using Oceananigans.BoundaryConditions
 using Oceananigans.Fields
@@ -7,17 +7,17 @@ using Oceananigans.AbstractOperations: KernelFunctionOperation
 using Oceananigans.ImmersedBoundaries
 
 """
-    struct HydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
+    struct SphericalCoriolis{S, FT} <: AbstractRotation
 
 A parameter object for constant rotation around a vertical axis on the sphere.
 """
-struct HydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
+struct SphericalCoriolis{S, FT} <: AbstractRotation
     rotation_rate :: FT
     scheme :: S
 end
 
 """
-    HydrostaticSphericalCoriolis([FT=Float64;]
+    SphericalCoriolis([FT=Float64;]
                                  rotation_rate = Ω_Earth,
                                  scheme = EnstrophyConserving())
 
@@ -29,70 +29,99 @@ Keyword arguments
 - `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref).
 - `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `EnstrophyConserving()` (default).
 """
-function HydrostaticSphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
+function SphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
                                       rotation_rate = Ω_Earth,
                                       scheme :: S = EnstrophyConserving(FT)) where S
 
-    return HydrostaticSphericalCoriolis{S, FT}(rotation_rate, scheme)
+    return SphericalCoriolis{S, FT}(rotation_rate, scheme)
 end
 
-Adapt.adapt_structure(to, coriolis::HydrostaticSphericalCoriolis) =
-    HydrostaticSphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
+Adapt.adapt_structure(to, coriolis::SphericalCoriolis) =
+    SphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
                                  Adapt.adapt(to, coriolis.scheme))
 
 @inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
 @inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
 @inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)
 
-@inline fᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis) =
+@inline fᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) =
     2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
+@inline f̃ᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) =
+    2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid))
 
-@inline z_f_cross_U(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, U) = zero(grid)
+# Check if the model is Nonhydrostatic
+# I know that want to do something different but this is place holder until
+# I have a better idea as to what people would like to do here.
 
-#####
-##### Active Point Enstrophy-conserving scheme
-#####
+nonhydrostatic = false #true
 
-# It might happen that a cell is active but all the neighbouring staggered nodes are inactive,
-# (an example is a 1-cell large channel)
-# In that case the Coriolis force is equal to zero
+if nonhydrostatic
+    @inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
+    @inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis, v) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
 
-const CoriolisEnstrophyConserving = HydrostaticSphericalCoriolis{<:EnstrophyConserving}
+    @inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
 
-@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
-    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
-                active_weighted_ℑxyᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+    @inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u)
 
-@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
-    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
-                active_weighted_ℑxyᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+    @inline x_f_cross_U(i, j, k, grid, U) =
+        @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)  
+                  + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
 
-#####
-##### Energy-conserving scheme
-#####
+    @inline y_f_cross_U(i, j, k, grid, U) =
+        @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+
+    @inline z_f_cross_U(i, j, k, grid, U) =
+        @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f̃_ℑz_uᶠᶠᵃ, U[1]) * Δz⁻¹ᶜᶜᶠ(i, j, k, grid)
+else
+    @inline z_f_cross_U(i, j, k, grid, coriolis::SphericalCoriolis, U) = zero(grid)
+
+    #####
+    ##### Active Point Enstrophy-conserving scheme
+    #####
+
+    # It might happen that a cell is active but all the neighbouring staggered nodes are inactive,
+    # (an example is a 1-cell large channel)
+    # In that case the Coriolis force is equal to zero
 
-const CoriolisEnergyConserving = HydrostaticSphericalCoriolis{<:EnergyConserving}
+    const CoriolisEnstrophyConserving = SphericalCoriolis{<:EnstrophyConserving}
 
-@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
-@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
+    @inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
+        @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
+                    active_weighted_ℑxyᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
 
-@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
-    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+    @inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
+        @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
+                    active_weighted_ℑxyᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+
+    #####
+    ##### Energy-conserving scheme
+    #####
+
+    const CoriolisEnergyConserving = SphericalCoriolis{<:EnergyConserving}
+
+    @inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
+    @inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
+
+    @inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
+        @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+
+    @inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
+        @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+
+end
 
-@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
-    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
 
 #####
 ##### Show
 #####
 
-function Base.show(io::IO, hydrostatic_spherical_coriolis::HydrostaticSphericalCoriolis)
-    coriolis_scheme = hydrostatic_spherical_coriolis.scheme
-    rotation_rate   = hydrostatic_spherical_coriolis.rotation_rate
+function Base.show(io::IO, spherical_coriolis::SphericalCoriolis)
+    coriolis_scheme = spherical_coriolis.scheme
+    rotation_rate   = spherical_coriolis.rotation_rate
     rotation_rate_Earth = rotation_rate / Ω_Earth
     rotation_rate_str = @sprintf("%.2e s⁻¹ = %.2e Ω_Earth", rotation_rate, rotation_rate_Earth)
 
-    return print(io, "HydrostaticSphericalCoriolis", '\n',
+    return print(io, "SphericalCoriolis", '\n',
                  "├─ rotation rate: ", rotation_rate_str, '\n',
                  "└─ scheme: ", summary(coriolis_scheme))
 end

From fe1044d458dac802745bd98430e18153dcaad460 Mon Sep 17 00:00:00 2001
From: "Francis J. Poulin" <fpoulin@uwaterloo.ca>
Date: Tue, 5 Aug 2025 15:29:26 -0400
Subject: [PATCH 06/17] Update src/Coriolis/spherical_coriolis.jl

Co-authored-by: Simone Silvestri <silvestri.simone0@gmail.com>
---
 src/Coriolis/spherical_coriolis.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index 398095dffc..a9ce6428bd 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -11,7 +11,7 @@ using Oceananigans.ImmersedBoundaries
 
 A parameter object for constant rotation around a vertical axis on the sphere.
 """
-struct SphericalCoriolis{S, FT} <: AbstractRotation
+struct SphericalCoriolis{H, S, FT} <: AbstractRotation
     rotation_rate :: FT
     scheme :: S
 end

From d4fc356b7adaf10c355f96ee7e6ce9aa8212d4e2 Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Tue, 5 Aug 2025 16:00:49 -0400
Subject: [PATCH 07/17] Trying to create a dispatch between hydrostatic and
 nonhydrostatic

---
 src/Coriolis/Coriolis.jl           |  7 +++
 src/Coriolis/spherical_coriolis.jl | 74 ++++++++----------------------
 2 files changed, 25 insertions(+), 56 deletions(-)

diff --git a/src/Coriolis/Coriolis.jl b/src/Coriolis/Coriolis.jl
index 1dd49bbaeb..c33b5bf29e 100644
--- a/src/Coriolis/Coriolis.jl
+++ b/src/Coriolis/Coriolis.jl
@@ -25,9 +25,16 @@ Abstract supertype for parameters related to background rotation rates.
 """
 abstract type AbstractRotation end
 
+abstract type HydrostaticCoriolis end
+
+abstract type NonhydrostaticCoriolis end
+
 const face = Face()
 const center = Center()
 
+const HydrostaticSphericalCoriolis{S, FT} = SphericalCoriolis{HydrostaticCoriolis, S, FT}
+const NonhydrostaticSphericalCoriolis{S, FT} =  SphericalCoriolis{NonhydrostaticCoriolis, S, FT}
+
 include("no_rotation.jl")
 include("f_plane.jl")
 include("constant_cartesian_coriolis.jl")
diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index a9ce6428bd..51e836519b 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -14,12 +14,14 @@ A parameter object for constant rotation around a vertical axis on the sphere.
 struct SphericalCoriolis{H, S, FT} <: AbstractRotation
     rotation_rate :: FT
     scheme :: S
+    HydrostaticCoriolis :: H
 end
 
 """
     SphericalCoriolis([FT=Float64;]
                                  rotation_rate = Ω_Earth,
                                  scheme = EnstrophyConserving())
+                                 HydrostaticCoriolis = 
 
 Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`.
 
@@ -28,7 +30,9 @@ Keyword arguments
 
 - `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref).
 - `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `EnstrophyConserving()` (default).
+- 'HydrostaticCoriolis`: `HydrostaticCoriolis` or `NonhydrostaticCoriolis`.
 """
+
 function SphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
                                       rotation_rate = Ω_Earth,
                                       scheme :: S = EnstrophyConserving(FT)) where S
@@ -44,72 +48,30 @@ Adapt.adapt_structure(to, coriolis::SphericalCoriolis) =
 @inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
 @inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)
 
-@inline fᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) =
+@inline fᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) = 
     2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
 @inline f̃ᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) =
     2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid))
 
-# Check if the model is Nonhydrostatic
-# I know that want to do something different but this is place holder until
-# I have a better idea as to what people would like to do here.
-
-nonhydrostatic = false #true
-
-if nonhydrostatic
-    @inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
-    @inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis, v) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
-
-    @inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
-
-    @inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u)
-
-    @inline x_f_cross_U(i, j, k, grid, U) =
-        @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)  
-                  + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
-
-    @inline y_f_cross_U(i, j, k, grid, U) =
-        @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
-
-    @inline z_f_cross_U(i, j, k, grid, U) =
-        @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f̃_ℑz_uᶠᶠᵃ, U[1]) * Δz⁻¹ᶜᶜᶠ(i, j, k, grid)
-else
-    @inline z_f_cross_U(i, j, k, grid, coriolis::SphericalCoriolis, U) = zero(grid)
+@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
+@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u)
+@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, u) = zero(grid)
 
-    #####
-    ##### Active Point Enstrophy-conserving scheme
-    #####
+@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
 
-    # It might happen that a cell is active but all the neighbouring staggered nodes are inactive,
-    # (an example is a 1-cell large channel)
-    # In that case the Coriolis force is equal to zero
+@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis, v) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
+@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, v) = zero(grid)
 
-    const CoriolisEnstrophyConserving = SphericalCoriolis{<:EnstrophyConserving}
+@inline x_f_cross_U(i, j, k, grid, U) =
+    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)  
+                + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
 
-    @inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
-        @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
-                    active_weighted_ℑxyᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+@inline y_f_cross_U(i, j, k, grid, U) =
+    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
 
-    @inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
-        @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
-                    active_weighted_ℑxyᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
-
-    #####
-    ##### Energy-conserving scheme
-    #####
-
-    const CoriolisEnergyConserving = SphericalCoriolis{<:EnergyConserving}
-
-    @inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
-    @inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
-
-    @inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
-        @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
-
-    @inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
-        @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
-
-end
 
+@inline z_f_cross_U(i, j, k, grid, U) =
+    @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f̃_ℑz_uᶠᶠᵃ, U[1]) * Δz⁻¹ᶜᶜᶠ(i, j, k, grid)
 
 #####
 ##### Show

From 960cf687f2716f712abf2578a888a87476929191 Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Tue, 5 Aug 2025 16:11:03 -0400
Subject: [PATCH 08/17] fix the coriolis term in the z equation

---
 src/Coriolis/spherical_coriolis.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index 51e836519b..acf7940c4b 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -59,8 +59,8 @@ Adapt.adapt_structure(to, coriolis::SphericalCoriolis) =
 
 @inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
 
-@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis, v) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
-@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, v) = zero(grid)
+@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis, w) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
+@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, w) = zero(grid)
 
 @inline x_f_cross_U(i, j, k, grid, U) =
     @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)  

From de89b49149b3bf8733a4dbfe112f047fb8e2f6c3 Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Wed, 6 Aug 2025 10:47:41 -0400
Subject: [PATCH 09/17] removed unnecessary file

---
 .../nonhydrostatic_spherical_coriolis.jl      | 71 -------------------
 1 file changed, 71 deletions(-)
 delete mode 100644 src/Coriolis/nonhydrostatic_spherical_coriolis.jl

diff --git a/src/Coriolis/nonhydrostatic_spherical_coriolis.jl b/src/Coriolis/nonhydrostatic_spherical_coriolis.jl
deleted file mode 100644
index b0a7fc2e51..0000000000
--- a/src/Coriolis/nonhydrostatic_spherical_coriolis.jl
+++ /dev/null
@@ -1,71 +0,0 @@
-using Oceananigans.Grids: LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, peripheral_node, φnode
-using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind, hack_cosd
-using Oceananigans.Advection: EnergyConserving, EnstrophyConserving
-using Oceananigans.BoundaryConditions
-using Oceananigans.Fields
-using Oceananigans.AbstractOperations: KernelFunctionOperation
-using Oceananigans.ImmersedBoundaries
-
-"""
-    struct NonhydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
-
-A parameter object for constant rotation around a vertical axis on the sphere.
-
-This is only designed to work for a LatitudeLongitudeGrid.
-
-There is no attempt to conserve energy or enstrophy at this early stage.
-"""
-struct NonhydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
-    rotation_rate :: FT
-end
-
-"""
-    NonhydrostaticSphericalCoriolis([FT=Float64;]
-                                 rotation_rate = Ω_Earth,
-
-Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`.
-
-Keyword arguments
-=================
-
-- `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref).
-"""
-function NonhydrostaticSphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
-                                      rotation_rate = Ω_Earth,
-
-    return NonhydrostaticSphericalCoriolis{FT}(rotation_rate)
-end
-
-Adapt.adapt_structure(to, coriolis::NonhydrostaticSphericalCoriolis) =
-    NonhydrostaticSphericalCoriolis(Adapt.adapt(to))
-
-@inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
-
-@inline fᶠᶠᵃ(i, j, k, grid) = 2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
-@inline f̃ᶠᶠᵃ(i, j, k, grid) = 2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid))
-
-@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
-@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis, v) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
-
-@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
-
-@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u)
-
-@inline x_f_cross_U(i, j, k, grid, U) =
-    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)  
-              + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
-
-@inline y_f_cross_U(i, j, k, grid, U) =
-    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
-
-@inline z_f_cross_U(i, j, k, grid, U) =
-    @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f̃_ℑz_uᶠᶠᵃ, U[1]) * Δz⁻¹ᶜᶜᶠ(i, j, k, grid)
-
-function Base.show(io::IO, nonhydrostatic_spherical_coriolis::NonhydrostaticSphericalCoriolis)
-    rotation_rate   = hydrostatic_spherical_coriolis.rotation_rate
-    rotation_rate_Earth = rotation_rate / Ω_Earth
-    rotation_rate_str = @sprintf("%.2e s⁻¹ = %.2e Ω_Earth", rotation_rate, rotation_rate_Earth)
-
-    return print(io, "NonhydrostaticSphericalCoriolis", '\n',
-                 "├─ rotation rate: ", rotation_rate_str, '\n')
-end

From 43d400dc80e106f9ead872284bcb64f3aee934e0 Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Wed, 6 Aug 2025 10:49:21 -0400
Subject: [PATCH 10/17] move definitions to the correct file

---
 src/Coriolis/Coriolis.jl           | 3 ---
 src/Coriolis/spherical_coriolis.jl | 3 +++
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/Coriolis/Coriolis.jl b/src/Coriolis/Coriolis.jl
index c33b5bf29e..a53fd681e3 100644
--- a/src/Coriolis/Coriolis.jl
+++ b/src/Coriolis/Coriolis.jl
@@ -32,9 +32,6 @@ abstract type NonhydrostaticCoriolis end
 const face = Face()
 const center = Center()
 
-const HydrostaticSphericalCoriolis{S, FT} = SphericalCoriolis{HydrostaticCoriolis, S, FT}
-const NonhydrostaticSphericalCoriolis{S, FT} =  SphericalCoriolis{NonhydrostaticCoriolis, S, FT}
-
 include("no_rotation.jl")
 include("f_plane.jl")
 include("constant_cartesian_coriolis.jl")
diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index acf7940c4b..bae1220e67 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -6,6 +6,9 @@ using Oceananigans.Fields
 using Oceananigans.AbstractOperations: KernelFunctionOperation
 using Oceananigans.ImmersedBoundaries
 
+const HydrostaticSphericalCoriolis{S, FT} = SphericalCoriolis{HydrostaticCoriolis, S, FT}
+const NonhydrostaticSphericalCoriolis{S, FT} =  SphericalCoriolis{NonhydrostaticCoriolis, S, FT}
+
 """
     struct SphericalCoriolis{S, FT} <: AbstractRotation
 

From 8707be52604455478566e6ef7c444ed2feec5ef5 Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Wed, 6 Aug 2025 10:57:56 -0400
Subject: [PATCH 11/17] move where const are defined

---
 src/Coriolis/spherical_coriolis.jl | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index bae1220e67..c7fab922b5 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -6,9 +6,6 @@ using Oceananigans.Fields
 using Oceananigans.AbstractOperations: KernelFunctionOperation
 using Oceananigans.ImmersedBoundaries
 
-const HydrostaticSphericalCoriolis{S, FT} = SphericalCoriolis{HydrostaticCoriolis, S, FT}
-const NonhydrostaticSphericalCoriolis{S, FT} =  SphericalCoriolis{NonhydrostaticCoriolis, S, FT}
-
 """
     struct SphericalCoriolis{S, FT} <: AbstractRotation
 
@@ -47,6 +44,10 @@ Adapt.adapt_structure(to, coriolis::SphericalCoriolis) =
     SphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
                                  Adapt.adapt(to, coriolis.scheme))
 
+const HydrostaticSphericalCoriolis{S, FT} = SphericalCoriolis{HydrostaticCoriolis, S, FT}
+const NonhydrostaticSphericalCoriolis{S, FT} =  SphericalCoriolis{NonhydrostaticCoriolis, S, FT}
+
+
 @inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
 @inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
 @inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)

From 0117481bf04fb2bd2d65ea4dd41e3cb2e14ab410 Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Wed, 6 Aug 2025 11:00:08 -0400
Subject: [PATCH 12/17] fix small error with i

---
 src/Coriolis/spherical_coriolis.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index c7fab922b5..2cc28f1b57 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -67,8 +67,7 @@ const NonhydrostaticSphericalCoriolis{S, FT} =  SphericalCoriolis{Nonhydrostatic
 @inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, w) = zero(grid)
 
 @inline x_f_cross_U(i, j, k, grid, U) =
-    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)  
-                + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid) + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
 
 @inline y_f_cross_U(i, j, k, grid, U) =
     @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)

From cb1fa35a2914e2f58caaf9393102553c7bb120c1 Mon Sep 17 00:00:00 2001
From: Francis Poulin <fpoulin@uwaterloo.ca>
Date: Thu, 14 Aug 2025 14:49:36 -0400
Subject: [PATCH 13/17] trying to set up the structure as suggested

---
 src/Coriolis/Coriolis.jl           |   4 --
 src/Coriolis/spherical_coriolis.jl | 110 ++++++++++++++++++++++++-----
 2 files changed, 94 insertions(+), 20 deletions(-)

diff --git a/src/Coriolis/Coriolis.jl b/src/Coriolis/Coriolis.jl
index a53fd681e3..1dd49bbaeb 100644
--- a/src/Coriolis/Coriolis.jl
+++ b/src/Coriolis/Coriolis.jl
@@ -25,10 +25,6 @@ Abstract supertype for parameters related to background rotation rates.
 """
 abstract type AbstractRotation end
 
-abstract type HydrostaticCoriolis end
-
-abstract type NonhydrostaticCoriolis end
-
 const face = Face()
 const center = Center()
 
diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index 2cc28f1b57..81e9d0044a 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -6,22 +6,36 @@ using Oceananigans.Fields
 using Oceananigans.AbstractOperations: KernelFunctionOperation
 using Oceananigans.ImmersedBoundaries
 
+
+"""
+Will use the following syntax:
+
+hydrostatic_coriolis = HydrostaticSphericalCoriolis()
+full_coriolis = SphericalCoriolis()
+
 """
-    struct SphericalCoriolis{S, FT} <: AbstractRotation
+
+struct HydrostaticFormulation end
+struct NonhydrostaticFormulation end
+
+
+"""
+    struct SphericalCoriolis{S, FT, F} 
 
 A parameter object for constant rotation around a vertical axis on the sphere.
 """
-struct SphericalCoriolis{H, S, FT} <: AbstractRotation
+
+struct SphericalCoriolis{S, FT, F}
     rotation_rate :: FT
     scheme :: S
-    HydrostaticCoriolis :: H
+    formulation :: F
 end
 
 """
     SphericalCoriolis([FT=Float64;]
                                  rotation_rate = Ω_Earth,
                                  scheme = EnstrophyConserving())
-                                 HydrostaticCoriolis = 
+                                 formulation = HydrostaticFormulation/NonhydrostaticFormulation 
 
 Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`.
 
@@ -30,12 +44,16 @@ Keyword arguments
 
 - `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref).
 - `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `EnstrophyConserving()` (default).
-- 'HydrostaticCoriolis`: `HydrostaticCoriolis` or `NonhydrostaticCoriolis`.
+- 'formulation`: `HydrostaticFormulation` or `NonhydrostaticFormulation`.
 """
 
+const HydrostaticSphericalCoriolis{S, FT} = SphericalCoriolis{S, FT, <:HydrostaticFormulation} where {S, FT}
+const SphericalCoriolis{S, FT} = SphericalCoriolis{S, FT, <:NonhydrostaticFormulation} where {S, FT}
+
 function SphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
                                       rotation_rate = Ω_Earth,
-                                      scheme :: S = EnstrophyConserving(FT)) where S
+                                      scheme :: S = EnstrophyConserving(FT),
+                                      ::NonhydrostaticFormulation) where {S, FT}
 
     return SphericalCoriolis{S, FT}(rotation_rate, scheme)
 end
@@ -44,10 +62,6 @@ Adapt.adapt_structure(to, coriolis::SphericalCoriolis) =
     SphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
                                  Adapt.adapt(to, coriolis.scheme))
 
-const HydrostaticSphericalCoriolis{S, FT} = SphericalCoriolis{HydrostaticCoriolis, S, FT}
-const NonhydrostaticSphericalCoriolis{S, FT} =  SphericalCoriolis{NonhydrostaticCoriolis, S, FT}
-
-
 @inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
 @inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
 @inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)
@@ -57,14 +71,16 @@ const NonhydrostaticSphericalCoriolis{S, FT} =  SphericalCoriolis{Nonhydrostatic
 @inline f̃ᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) =
     2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid))
 
-@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
-@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u)
-@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, u) = zero(grid)
 
-@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
 
-@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::NonhydrostaticSphericalCoriolis, w) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
-@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis, w) = zero(grid)
+@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
+@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, formulation::NonhydrostaticFormulation, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u)
+@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, formulation::HydrostaticFormulation, u) = zero(grid)
+
+@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
+
+@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, formulation::NonhydrostaticFormulation, w) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
+@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, formulation::HydrostaticFormulation, w) = zero(grid)
 
 @inline x_f_cross_U(i, j, k, grid, U) =
     @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid) + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
@@ -76,17 +92,79 @@ const NonhydrostaticSphericalCoriolis{S, FT} =  SphericalCoriolis{Nonhydrostatic
 @inline z_f_cross_U(i, j, k, grid, U) =
     @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f̃_ℑz_uᶠᶠᵃ, U[1]) * Δz⁻¹ᶜᶜᶠ(i, j, k, grid)
 
+
+
+
+
+
+function SphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
+                                      rotation_rate = Ω_Earth,
+                                      scheme :: S = EnstrophyConserving(FT),
+                                      formulation::HydrostaticFormulation) where S
+
+    return HydrostaticSphericalCoriolis{S, FT}(rotation_rate, scheme)
+end
+
+Adapt.adapt_structure(to, coriolis, formulation::HydrostaticFormulation) =
+    HydrostaticSphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
+                                 Adapt.adapt(to, coriolis.scheme))
+
+@inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
+@inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
+@inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)
+
+@inline fᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, formation::HydrostaticFormulation) =
+    2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
+
+@inline z_f_cross_U(i, j, k, grid, coriolis::SphericalCoriolis, formulation::HydrostaticFormulationU) = zero(grid)
+
+#####
+##### Active Point Enstrophy-conserving scheme
+#####
+
+# It might happen that a cell is active but all the neighbouring staggered nodes are inactive,
+# (an example is a 1-cell large channel)
+# In that case the Coriolis force is equal to zero
+
+const CoriolisEnstrophyConserving = SphericalCoriolis{<:HydrostaticFormulation, <:EnstrophyConserving}
+
+@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
+    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
+                active_weighted_ℑxyᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+
+@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) =
+    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
+                active_weighted_ℑxyᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+
+#####
+##### Energy-conserving scheme
+#####
+
+const CoriolisEnergyConserving = SphericalCoriolis{<:HydrostaticFormulation, <:EnergyConserving}
+
+@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
+@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
+
+@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
+    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+
+@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
+    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+
+
 #####
 ##### Show
 #####
 
 function Base.show(io::IO, spherical_coriolis::SphericalCoriolis)
     coriolis_scheme = spherical_coriolis.scheme
+    coriolis_formulation = spherical_coriolis.formulation
     rotation_rate   = spherical_coriolis.rotation_rate
     rotation_rate_Earth = rotation_rate / Ω_Earth
     rotation_rate_str = @sprintf("%.2e s⁻¹ = %.2e Ω_Earth", rotation_rate, rotation_rate_Earth)
 
     return print(io, "SphericalCoriolis", '\n',
                  "├─ rotation rate: ", rotation_rate_str, '\n',
+                 "├─ formulation: ", summary(coriolis_formulation), '\n',
                  "└─ scheme: ", summary(coriolis_scheme))
 end

From d2a78295a55804d37e40e7534672328a1265268b Mon Sep 17 00:00:00 2001
From: Simone Silvestri <silvestri.simone0@gmail.com>
Date: Mon, 18 Aug 2025 22:50:17 +0200
Subject: [PATCH 14/17] Update spherical_coriolis.jl

I have corrected a couple of methods, now it should at least run, we have to check that the location of the operations are correct
---
 src/Coriolis/spherical_coriolis.jl | 70 ++++++++----------------------
 1 file changed, 18 insertions(+), 52 deletions(-)

diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index 81e9d0044a..500ae78320 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -48,75 +48,41 @@ Keyword arguments
 """
 
 const HydrostaticSphericalCoriolis{S, FT} = SphericalCoriolis{S, FT, <:HydrostaticFormulation} where {S, FT}
-const SphericalCoriolis{S, FT} = SphericalCoriolis{S, FT, <:NonhydrostaticFormulation} where {S, FT}
 
 function SphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
-                                      rotation_rate = Ω_Earth,
-                                      scheme :: S = EnstrophyConserving(FT),
-                                      ::NonhydrostaticFormulation) where {S, FT}
+                           rotation_rate = Ω_Earth,
+                           scheme = EnstrophyConserving(FT),
+                           formulation = HydrostaticFormulation()) 
 
-    return SphericalCoriolis{S, FT}(rotation_rate, scheme)
+    return SphericalCoriolis(rotation_rate, scheme, formulation)
 end
 
 Adapt.adapt_structure(to, coriolis::SphericalCoriolis) =
     SphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
-                                 Adapt.adapt(to, coriolis.scheme))
+                      Adapt.adapt(to, coriolis.scheme),
+                      Adapt.adapt(to, coriolis.formulation))
 
 @inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
 @inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
 @inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)
 
-@inline fᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) = 
-    2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
-@inline f̃ᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) =
-    2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid))
+@inline fᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) = 2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
+@inline f̃ᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis) = 2 * coriolis.rotation_rate * hack_cosd(φᶠᶠᵃ(i, j, k, grid)
 
+@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
+@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
 
+@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u)
+@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis, w) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
 
-@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
-@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, formulation::NonhydrostaticFormulation, u) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑzᵃᶠᵃ(i, j, k, grid, Δz_qᶜᶜᶠ, u)
-@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, formulation::HydrostaticFormulation, u) = zero(grid)
-
-@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
-
-@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, formulation::NonhydrostaticFormulation, w) = f̃ᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, w)
-@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, formulation::HydrostaticFormulation, w) = zero(grid)
-
-@inline x_f_cross_U(i, j, k, grid, U) =
-    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid) + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
-
-@inline y_f_cross_U(i, j, k, grid, U) =
-    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
-
-
-@inline z_f_cross_U(i, j, k, grid, U) =
-    @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f̃_ℑz_uᶠᶠᵃ, U[1]) * Δz⁻¹ᶜᶜᶠ(i, j, k, grid)
-
-
-
-
-
-
-function SphericalCoriolis(FT::DataType=Oceananigans.defaults.FloatType;
-                                      rotation_rate = Ω_Earth,
-                                      scheme :: S = EnstrophyConserving(FT),
-                                      formulation::HydrostaticFormulation) where S
-
-    return HydrostaticSphericalCoriolis{S, FT}(rotation_rate, scheme)
-end
-
-Adapt.adapt_structure(to, coriolis, formulation::HydrostaticFormulation) =
-    HydrostaticSphericalCoriolis(Adapt.adapt(to, coriolis.rotation_rate),
-                                 Adapt.adapt(to, coriolis.scheme))
-
-@inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid)        = φnode(j, grid, face)
-@inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, face, face)
-@inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid)         = φᶠᶠᵃ(i, j, k, grid.underlying_grid)
+@inline f̃_ℑz_uᶠᶠᵃ(i, j, k, grid, ::HydrostaticSphericalCoriolis, u) = zero(grid)
+@inline f̃_ℑx_wᶠᶠᵃ(i, j, k, grid, ::HydrostaticSphericalCoriolis, w) = zero(grid)
 
-@inline fᶠᶠᵃ(i, j, k, grid, coriolis::SphericalCoriolis, formation::HydrostaticFormulation) =
-    2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
+@inline x_f_cross_U(i, j, k, grid, U) = @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, U[2]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid) + ℑzᵃᶜᵃ(i, j, k, grid, f̃_ℑx_wᶠᶠᵃ, U[3]) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+@inline y_f_cross_U(i, j, k, grid, U) = @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, U[1]) * Δy⁻¹ᶜᶠᶜ(i, j, k, grid)
+@inline z_f_cross_U(i, j, k, grid, U) = @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f̃_ℑz_uᶠᶠᵃ, U[1]) * Δz⁻¹ᶜᶜᶠ(i, j, k, grid)
 
-@inline z_f_cross_U(i, j, k, grid, coriolis::SphericalCoriolis, formulation::HydrostaticFormulationU) = zero(grid)
+@inline z_f_cross_U(i, j, k, grid, ::HydrostaticSphericalCoriolis, U) = zero(grid)
 
 #####
 ##### Active Point Enstrophy-conserving scheme

From c7e4b3197810c39776945ab44c2ea7ae23efd727 Mon Sep 17 00:00:00 2001
From: "Francis J. Poulin" <fpoulin@uwaterloo.ca>
Date: Wed, 20 Aug 2025 16:30:59 -0400
Subject: [PATCH 15/17] Update src/Coriolis/spherical_coriolis.jl

Co-authored-by: Simone Silvestri <silvestri.simone0@gmail.com>
---
 src/Coriolis/spherical_coriolis.jl | 9 ---------
 1 file changed, 9 deletions(-)

diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index 500ae78320..f5c01408bc 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -6,15 +6,6 @@ using Oceananigans.Fields
 using Oceananigans.AbstractOperations: KernelFunctionOperation
 using Oceananigans.ImmersedBoundaries
 
-
-"""
-Will use the following syntax:
-
-hydrostatic_coriolis = HydrostaticSphericalCoriolis()
-full_coriolis = SphericalCoriolis()
-
-"""
-
 struct HydrostaticFormulation end
 struct NonhydrostaticFormulation end
 

From e96172e1f445f08cd6e1433c74a9199e398ac0e2 Mon Sep 17 00:00:00 2001
From: "Francis J. Poulin" <fpoulin@uwaterloo.ca>
Date: Wed, 20 Aug 2025 16:31:19 -0400
Subject: [PATCH 16/17] Update src/Coriolis/spherical_coriolis.jl

Co-authored-by: Simone Silvestri <silvestri.simone0@gmail.com>
---
 src/Coriolis/spherical_coriolis.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index f5c01408bc..020c0d2849 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -9,7 +9,6 @@ using Oceananigans.ImmersedBoundaries
 struct HydrostaticFormulation end
 struct NonhydrostaticFormulation end
 
-
 """
     struct SphericalCoriolis{S, FT, F} 
 

From 1b65af7ddcc5670e17ec4fbff3dd6b82e7400b1a Mon Sep 17 00:00:00 2001
From: "Francis J. Poulin" <fpoulin@uwaterloo.ca>
Date: Wed, 20 Aug 2025 16:31:36 -0400
Subject: [PATCH 17/17] Update src/Coriolis/spherical_coriolis.jl

Co-authored-by: Simone Silvestri <silvestri.simone0@gmail.com>
---
 src/Coriolis/spherical_coriolis.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Coriolis/spherical_coriolis.jl b/src/Coriolis/spherical_coriolis.jl
index 020c0d2849..e9ad5ad3a9 100644
--- a/src/Coriolis/spherical_coriolis.jl
+++ b/src/Coriolis/spherical_coriolis.jl
@@ -25,7 +25,7 @@ end
     SphericalCoriolis([FT=Float64;]
                                  rotation_rate = Ω_Earth,
                                  scheme = EnstrophyConserving())
-                                 formulation = HydrostaticFormulation/NonhydrostaticFormulation 
+                                 formulation = HydrostaticFormulation()
 
 Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`.