do not change OceanSimulations.jl here

src/OceanSimulations/OceanSimulations.jl

@@ -16,7 +16,6 @@ using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities:
     CATKEEquation
 
 using SeawaterPolynomials.TEOS10: TEOS10EquationOfState
-using Statistics: mean
 
 using Oceananigans.BuoyancyModels: g_Earth
 using Oceananigans.Coriolis: Ω_Earth
@@ -43,30 +42,10 @@ default_or_override(override, alternative_default=nothing) = override
 # Some defaults
 default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7)
 
-function infer_maximum_Δt(grid)
-    Δx = mean(xspacings(grid))
-    Δy = mean(yspacings(grid))
-    Δθ = rad2deg(mean([Δx, Δy])) / grid.radius
-
-    # The maximum Δt is roughly 40minutes / Δθ, giving:
-    # - 40 minutes for a 1 degree ocean
-    # - 20 minutes for a 1/4 degree ocean
-    # - 10 minutes for a 1/8 degree ocean
-    # - 5 minutes for a 1/16 degree ocean
-    # - 2.5 minutes for a 1/32 degree ocean
-
-    return 40minutes / Δθ
-end
-
+# 70 substeps is a safe rule of thumb for an ocean at 1/4 - 1/10th of a degree
+# TODO: pass the cfl and a given Δt to calculate the number of substeps?
 const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}}
-
-function default_free_surface(grid::TripolarOfSomeKind; 
-                              fixed_Δt = compute_maximum_Δt(grid),
-                              cfl = 0.7) 
-    free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt)
-    @info "Using a $(free_surface)"
-    return free_surface
-end
+default_free_surface(grid::TripolarOfSomeKind) = SplitExplicitFreeSurface(grid; substeps=70)
 
 function default_ocean_closure()
     mixing_length = CATKEMixingLength(Cᵇ=0.01)
@@ -96,7 +75,7 @@ default_tracer_advection() = FluxFormAdvection(WENO(order=7),
 # TODO: Specify the grid to a grid on the sphere; otherwise we can provide a different
 # function that requires latitude and longitude etc for computing coriolis=FPlane...
 function ocean_simulation(grid;
-                          Δt = infer_maximum_Δt(grid),
+                          Δt = 5minutes,
                           closure = default_ocean_closure(),
                           tracers = (:T, :S),
                           free_surface = default_free_surface(grid),