diff --git a/.buildkite/examples_build.yml b/.buildkite/examples_build.yml index b5be821e..9101af75 100644 --- a/.buildkite/examples_build.yml +++ b/.buildkite/examples_build.yml @@ -1,7 +1,7 @@ agents: queue: new-central slurm_mem: 8G - modules: climacommon/2024_05_27 + modules: climacommon/2024_10_10 slurm_time: 48:00:00 env: diff --git a/docs/Project.toml b/docs/Project.toml index 4a554956..0402899d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -9,9 +9,9 @@ OrthogonalSphericalShellGrids = "c2be9673-fb75-4747-82dc-aa2bb9f4aed0" SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40" [compat] -CairoMakie = "0.10.12, 0.11, 0.12" +CairoMakie = "0.10.12, 0.11, 0.12, 0.13" DataDeps = "0.7" Documenter = "1" -Oceananigans = "0.95.7 - 0.99" +Oceananigans = "0.95.12 - 0.99" OrthogonalSphericalShellGrids = "0.2.1" -SeawaterPolynomials = "0.3.4" +SeawaterPolynomials = "0.3.5" diff --git a/docs/make.jl b/docs/make.jl index 4692edec..f543a89a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -17,6 +17,7 @@ to_be_literated = [ "generate_bathymetry.jl", "generate_surface_fluxes.jl", "single_column_os_papa_simulation.jl", + "one_degree_simulation.jl", # "mediterranean_simulation_with_ecco_restoring.jl", "near_global_ocean_simulation.jl" ] @@ -46,6 +47,7 @@ pages = [ "Surface fluxes" => "literated/generate_surface_fluxes.md", "Single-column simulation" => "literated/single_column_os_papa_simulation.md", # "Mediterranean simulation with ECCO restoring" => "literated/mediterranean_simulation_with_ecco_restoring.md", + "One-degree Ocean simulation" => "literated/one_degree_simulation.md", "Near-global Ocean simulation" => "literated/near_global_ocean_simulation.md", ], diff --git a/examples/ecco_mixed_layer_depth.jl b/examples/ecco_mixed_layer_depth.jl index 7eed00a5..643367a2 100644 --- a/examples/ecco_mixed_layer_depth.jl +++ b/examples/ecco_mixed_layer_depth.jl @@ -11,8 +11,8 @@ using SeawaterPolynomials: TEOS10EquationOfState using Oceananigans.BuoyancyFormulations: buoyancy arch = CPU() -Nx = 360 ÷ 1 -Ny = 160 ÷ 1 +Nx = 360 +Ny = 160 z = ClimaOcean.DataWrangling.ECCO.ECCO_z z = z[20:end] diff --git a/examples/generate_surface_fluxes.jl b/examples/generate_surface_fluxes.jl index 4cfd35dd..e97f6fad 100644 --- a/examples/generate_surface_fluxes.jl +++ b/examples/generate_surface_fluxes.jl @@ -22,6 +22,8 @@ using CairoMakie grid = ECCO_immersed_grid() +# We visualize the bottom height of the ECCO grid using CairoMakie. + fig = Figure() ax = Axis(fig[1, 1]) heatmap!(ax, interior(grid.immersed_boundary.bottom_height, :, :, 1)) @@ -45,7 +47,7 @@ save("ECCO_continents.png", fig) #hide # January 1st (at 00:00 AM and 03:00 AM). atmosphere = JRA55PrescribedAtmosphere(1:2; backend = InMemory()) -ocean = ocean_simulation(grid) +ocean = ocean_simulation(grid, closure=nothing) # Now that we have an atmosphere and ocean, we `set!` the ocean temperature and salinity # to the ECCO2 data by first creating T, S metadata objects, @@ -68,7 +70,7 @@ coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation=Radiation()) # Now that the surface fluxes are computed, we can extract and visualize them. # The turbulent fluxes are stored in `coupled_model.fluxes.turbulent`. -fluxes = coupled_model.fluxes.turbulent.fields +fluxes = coupled_model.interfaces.atmosphere_ocean_interface.fluxes λ, φ, z = nodes(fluxes.sensible_heat) fig = Figure(size = (800, 800), fontsize = 15) diff --git a/examples/one_degree_simulation.jl b/examples/one_degree_simulation.jl new file mode 100644 index 00000000..977c659a --- /dev/null +++ b/examples/one_degree_simulation.jl @@ -0,0 +1,260 @@ +# # One-degree global ocean simulation +# +# This example configures a global ocean--sea ice simulation at 1ᵒ horizontal resolution with +# realistic bathymetry and some closures. +# +# For this example, we need Oceananigans, ClimaOcean, OrthogonalSphericalShellGrids, and +# CairoMakie to visualize the simulation. Also we need CFTime and Dates for date handling. + +using ClimaOcean +using ClimaOcean.ECCO +using Oceananigans +using Oceananigans.Units +using OrthogonalSphericalShellGrids +using CFTime +using Dates +using Printf + +arch = GPU() + +# ### Grid and Bathymetry + +Nx = 360 +Ny = 180 +Nz = 100 + +r_faces = exponential_z_faces(; Nz, depth=5000, h=34) +z_faces = Oceananigans.MutableVerticalDiscretization(r_faces) + +underlying_grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (5, 5, 4), + first_pole_longitude = 70, + north_poles_latitude = 55) + +bottom_height = regrid_bathymetry(underlying_grid; + minimum_depth = 10, + interpolation_passes = 75, + major_basins = 2) + +# For this bathymetry at this horizontal resolution we need to manually open the Gibraltar strait. +view(bottom_height, 102:103, 124, 1) .= -400 +grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); active_cells_map=true) + +# ### Restoring +# +# We include temperature and salinity surface restoring to ECCO data. + +restoring_rate = 1 / 10days +z_below_surface = r_faces[end-1] + +mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(z_below_surface, 0)) + +dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1994, 1, 1) +temperature = ECCOMetadata(:temperature; dates, version=ECCO4Monthly(), dir="./") +salinity = ECCOMetadata(:salinity; dates, version=ECCO4Monthly(), dir="./") + +FT = ECCORestoring(temperature, grid; mask, rate=restoring_rate) +FS = ECCORestoring(salinity, grid; mask, rate=restoring_rate) +forcing = (T=FT, S=FS) + +# ### Closures +# +# We include a Gent-McWilliam isopycnal diffusivity as a parameterization for the mesoscale +# eddy fluxes. For vertical mixing at the upper-ocean boundary layer we include the CATKE +# parameterization. We also include some explicit horizontal diffusivity. + +using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity, + DiffusiveFormulation + +eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3, + skew_flux_formulation=DiffusiveFormulation()) +vertical_mixing = ClimaOcean.OceanSimulations.default_ocean_closure() + +closure = (eddy_closure, vertical_mixing) + +# ### Ocean simulation +# +# Now we bring everything together to construct the ocean simulation. +# We use a split-explicit timestepping with 30 substeps for the barotropic +# mode. + +free_surface = SplitExplicitFreeSurface(grid; substeps=30) + +momentum_advection = WENOVectorInvariant(vorticity_order=3) +tracer_advection = Centered() + +ocean = ocean_simulation(grid; + momentum_advection, + tracer_advection, + closure, + forcing, + free_surface) + +# ### Initial condition + +# We initialize the ocean from the ECCO state estimate. + +set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)), + S=ECCOMetadata(:salinity; dates=first(dates))) + +# ### Atmospheric forcing + +# We force the simulation with an JRA55-do atmospheric reanalysis. + +radiation = Radiation(arch) +atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(20)) + +# ### Coupled simulation + +# Now we are ready to build the coupled ocean--sea ice model and bring everything +# together into a `simulation`. + +# We use a relatively short time step initially and only run for a few days to +# avoid numerical instabilities from the initial "shock" of the adjustment of the +# flow fields. + +coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) +simulation = Simulation(coupled_model; Δt=1minutes, stop_time=10days) + +# ### A progress messenger +# +# We write a function that prints out a helpful progress message while the simulation runs. + +wall_time = Ref(time_ns()) + +function progress(sim) + ocean = sim.model.ocean + u, v, w = ocean.model.velocities + T = ocean.model.tracers.T + Tmax = maximum(interior(T)) + Tmin = minimum(interior(T)) + umax = (maximum(abs, interior(u)), + maximum(abs, interior(v)), + maximum(abs, interior(w))) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max|u|: (%.2e, %.2e, %.2e) m s⁻¹, ", umax...) + msg3 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg4 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg3 * msg4 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. + +add_callback!(simulation, progress, IterationInterval(10)) + +# ### Output +# +# We are almost there! We need to save some output. Below we choose to save _only surface_ +# fields using the `indices` keyword argument. We save all velocity and tracer components. +# Note, that besides temperature and salinity, the CATKE vertical mixing parameterization +# also uses a prognostic turbulent kinetic energy, `e`, to diagnose the vertical mixing length. + +outputs = merge(ocean.model.tracers, ocean.model.velocities) +ocean.output_writers[:surface] = JLD2OutputWriter(ocean.model, outputs; + schedule = TimeInterval(5days), + filename = "global_surface_fields", + indices = (:, :, grid.Nz), + overwrite_existing = true) + +# ### Ready to run + +# We are ready to press the big red button and run the simulation. + +# After we run for a short time (here we set up the simulation with `stop_time = 10days`), +# we increase the timestep and run for longer. + +run!(simulation) + +simulation.Δt = 20minutes +simulation.stop_time = 360days + +run!(simulation) + +# ### A pretty movie +# +# We load the saved output and make a pretty movie of the simulation. First we plot a snapshot: + +using CairoMakie + +u = FieldTimeSeries("global_surface_fields.jld2", "u"; backend = OnDisk()) +v = FieldTimeSeries("global_surface_fields.jld2", "v"; backend = OnDisk()) +T = FieldTimeSeries("global_surface_fields.jld2", "T"; backend = OnDisk()) +e = FieldTimeSeries("global_surface_fields.jld2", "e"; backend = OnDisk()) + +times = u.times +Nt = length(times) + +n = Observable(Nt) + +# We create a land mask and use it to fill land points with `NaN`s. + +land = interior(T.grid.immersed_boundary.bottom_height) .≥ 0 + +Tn = @lift begin + Tn = interior(T[$n]) + Tn[land] .= NaN + view(Tn, :, :, 1) +end + +en = @lift begin + en = interior(e[$n]) + en[land] .= NaN + view(en, :, :, 1) +end + +# We compute the surface speed. + +un = Field{Face, Center, Nothing}(u.grid) +vn = Field{Center, Face, Nothing}(v.grid) +s = Field(sqrt(un^2 + vn^2)) + +sn = @lift begin + parent(un) .= parent(u[$n]) + parent(vn) .= parent(v[$n]) + compute!(s) + sn = interior(s) + sn[land] .= NaN + view(sn, :, :, 1) +end + +# Finally, we plot a snapshot of the surface speed, temperature, and the turbulent +# eddy kinetic energy from the CATKE vertical mixing parameterization. + +fig = Figure(size = (800, 1200)) + +axs = Axis(fig[1, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)") +axT = Axis(fig[2, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)") +axe = Axis(fig[3, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)") + +hm = heatmap!(axs, sn, colorrange = (0, 0.5), colormap = :deep, nan_color=:lightgray) +Colorbar(fig[1, 2], hm, label = "Surface speed (m s⁻¹)") + +hm = heatmap!(axT, Tn, colorrange = (-1, 30), colormap = :magma, nan_color=:lightgray) +Colorbar(fig[2, 2], hm, label = "Surface Temperature (ᵒC)") + +hm = heatmap!(axe, en, colorrange = (0, 1e-3), colormap = :solar, nan_color=:lightgray) +Colorbar(fig[3, 2], hm, label = "Turbulent Kinetic Energy (m² s⁻²)") + +save("global_snapshot.png", fig) +nothing #hide + +# ![](global_snapshot.png) + +# And now a movie: + +record(fig, "one_degree_global_ocean_surface.mp4", 1:Nt, framerate = 8) do nn + n[] = nn +end +nothing #hide + +# ![](one_degree_global_ocean_surface.mp4) diff --git a/examples/single_column_os_papa_simulation.jl b/examples/single_column_os_papa_simulation.jl index 7c421efc..3dc7f43c 100644 --- a/examples/single_column_os_papa_simulation.jl +++ b/examples/single_column_os_papa_simulation.jl @@ -101,8 +101,7 @@ current_figure() # the skin temperature from a balance between internal and external heat fluxes. radiation = Radiation() -similarity_theory = SimilarityTheoryTurbulentFluxes(grid; surface_temperature_type=SkinTemperature()) -coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) +coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) simulation = Simulation(coupled_model, Δt=ocean.Δt, stop_time=30days) wall_clock = Ref(time_ns()) @@ -122,9 +121,9 @@ function progress(sim) S = sim.model.ocean.model.tracers.S e = sim.model.ocean.model.tracers.e - τx = first(sim.model.fluxes.total.ocean.momentum.u) - τy = first(sim.model.fluxes.total.ocean.momentum.v) - Q = first(sim.model.fluxes.total.ocean.heat) + τx = first(sim.model.interfaces.net_fluxes.ocean_surface.u) + τy = first(sim.model.interfaces.net_fluxes.ocean_surface.v) + Q = first(sim.model.interfaces.net_fluxes.ocean_surface.Q) u★ = sqrt(sqrt(τx^2 + τy^2)) @@ -142,15 +141,15 @@ end simulation.callbacks[:progress] = Callback(progress, IterationInterval(100)) # Build flux outputs -τx = coupled_model.fluxes.total.ocean.momentum.u -τy = coupled_model.fluxes.total.ocean.momentum.v -JT = coupled_model.fluxes.total.ocean.tracers.T -Js = coupled_model.fluxes.total.ocean.tracers.S -E = coupled_model.fluxes.turbulent.fields.water_vapor -Qc = coupled_model.fluxes.turbulent.fields.sensible_heat -Qv = coupled_model.fluxes.turbulent.fields.latent_heat -ρₒ = coupled_model.fluxes.ocean_reference_density -cₚ = coupled_model.fluxes.ocean_heat_capacity +τx = coupled_model.interfaces.atmosphere_ocean_interface.fluxes.x_momentum +τy = coupled_model.interfaces.atmosphere_ocean_interface.fluxes.y_momentum +JT = coupled_model.interfaces.net_fluxes.ocean_surface.T +Js = coupled_model.interfaces.net_fluxes.ocean_surface.S +E = coupled_model.interfaces.atmosphere_ocean_interface.fluxes.water_vapor +Qc = coupled_model.interfaces.atmosphere_ocean_interface.fluxes.sensible_heat +Qv = coupled_model.interfaces.atmosphere_ocean_interface.fluxes.latent_heat +ρₒ = coupled_model.interfaces.ocean_properties.reference_density +cₚ = coupled_model.interfaces.ocean_properties.heat_capacity Q = ρₒ * cₚ * JT ρτx = ρₒ * τx @@ -249,7 +248,7 @@ tn = @lift times[$n] colors = Makie.wong_colors() -ρₒ = coupled_model.fluxes.ocean_reference_density +ρₒ = coupled_model.interfaces.ocean_properties.reference_density τx = interior(ρτx, 1, 1, 1, :) ./ ρₒ τy = interior(ρτy, 1, 1, 1, :) ./ ρₒ u★ = @. (τx^2 + τy^2)^(1/4) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl deleted file mode 100644 index 03db70ca..00000000 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ /dev/null @@ -1,142 +0,0 @@ -using ClimaOcean -using ClimaOcean.ECCO: ECCO4Monthly, NearestNeighborInpainting -using OrthogonalSphericalShellGrids -using Oceananigans -using Oceananigans.Units -using CFTime -using Dates -using Printf -using CUDA: @allowscalar, device! - -using Oceananigans.Grids: znode - -arch = GPU() - -##### -##### Grid and Bathymetry -##### - -Nx = 360 -Ny = 180 -Nz = 100 - -z_faces = exponential_z_faces(; Nz, depth=5000, h=34) - -underlying_grid = TripolarGrid(arch; - size = (Nx, Ny, Nz), - z = z_faces, - first_pole_longitude = 70, - north_poles_latitude = 55) - -bottom_height = regrid_bathymetry(underlying_grid; - minimum_depth = 10, - interpolation_passes = 75, - major_basins = 2) - -# Open Gibraltar strait -# TODO: find a better way to do this -tampered_bottom_height = deepcopy(bottom_height) -view(tampered_bottom_height, 102:103, 124, 1) .= -400 - -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_height); active_cells_map=true) - -##### -##### Closures -##### - -gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000) -catke = ClimaOcean.OceanSimulations.default_ocean_closure() -viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=2000) - -closure = (gm, catke, viscous_closure) - -##### -##### Restoring -##### - -restoring_rate = 1 / 2days -z_below_surface = @allowscalar znode(1, 1, grid.Nz, grid, Center(), Center(), Face()) - -mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(z_below_surface, 0)) - -dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1993, 11, 1) -temperature = ECCOMetadata(:temperature; dates, version=ECCO4Monthly(), dir="./") -salinity = ECCOMetadata(:salinity; dates, version=ECCO4Monthly(), dir="./") - -# inpainting = NearestNeighborInpainting(30) should be enough to fill the gaps near bathymetry -FT = ECCORestoring(arch, temperature; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(50)) -FS = ECCORestoring(arch, salinity; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(50)) -forcing = (T=FT, S=FS) - -##### -##### Ocean simulation -##### - -momentum_advection = VectorInvariant() -tracer_advection = Centered(order=2) - -free_surface = SplitExplicitFreeSurface(grid; substeps=30) - -# Should we add a side drag since this is at a coarser resolution? -ocean = ocean_simulation(grid; momentum_advection, tracer_advection, - free_surface, forcing) - -set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)), - S=ECCOMetadata(:salinity; dates=first(dates))) - -##### -##### Sea ice simulation -##### - -sea_ice = sea_ice_simulation(grid) - -set!(sea_ice.model.ice_thickness, ECCOMetadata(:sea_ice_thickness; dates=first(dates)); inpainting=NearestNeighborInpainting(0)) -set!(sea_ice.model.ice_concentration, ECCOMetadata(:sea_ice_concentration; dates=first(dates)); inpainting=NearestNeighborInpainting(0)) - -##### -##### Atmospheric forcing -##### - -radiation = Radiation(arch) -atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(20)) - -##### -##### Coupled simulation -##### - -coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -simulation = Simulation(coupled_model; Δt=1minutes, stop_time=2*365days) - -##### -##### Run it! -##### - -wall_time = Ref(time_ns()) - -function progress(sim) - ocean = sim.model.ocean - sea_ice = sim.model.sea_ice - u, v, w = ocean.model.velocities - T = ocean.model.tracers.T - h = sea_ice.model.ice_thickness - Tmax = maximum(interior(T)) - Tmin = minimum(interior(T)) - hmax = maximum(interior(h)) - umax = (maximum(abs, interior(u)), - maximum(abs, interior(v)), - maximum(abs, interior(w))) - - step_time = 1e-9 * (time_ns() - wall_time[]) - - @info @sprintf("Time: %s, n: %d, Δt: %s, max|u|: (%.2e, %.2e, %.2e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, max(h): %.2f m, wall time: %s \n", - prettytime(sim), iteration(sim), prettytime(sim.Δt), - umax..., Tmax, Tmin, hmax, prettytime(step_time)) - - wall_time[] = time_ns() - - return nothing -end - -add_callback!(simulation, progress, IterationInterval(10)) - -run!(simulation) diff --git a/src/OceanSimulations.jl b/src/OceanSimulations.jl index 5b698c10..ade751a6 100644 --- a/src/OceanSimulations.jl +++ b/src/OceanSimulations.jl @@ -5,18 +5,23 @@ export ocean_simulation using Oceananigans using Oceananigans.Units using Oceananigans.Utils: with_tracers +using Oceananigans.Grids: architecture using Oceananigans.Advection: FluxFormAdvection +using Oceananigans.DistributedComputations: DistributedGrid, all_reduce using Oceananigans.BoundaryConditions: DefaultBoundaryCondition using Oceananigans.Coriolis: ActiveCellEnstrophyConserving using Oceananigans.ImmersedBoundaries: immersed_peripheral_node, inactive_node, MutableGridOfSomeKind using OrthogonalSphericalShellGrids +using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization + using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, CATKEMixingLength, CATKEEquation using SeawaterPolynomials.TEOS10: TEOS10EquationOfState +using Statistics: mean using Oceananigans.BuoyancyFormulations: g_Earth using Oceananigans.Coriolis: Ω_Earth @@ -43,18 +48,53 @@ default_or_override(override, alternative_default=nothing) = override # Some defaults default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7) -# 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? +function estimate_maximum_Δt(grid) + arch = architecture(grid) + Δx = mean(xspacings(grid)) + Δy = mean(yspacings(grid)) + Δθ = rad2deg(mean([Δx, Δy])) / grid.radius + + # The maximum Δt is roughly 30minutes / Δθ, giving: + # - 30 minutes for a 1 degree ocean + # - 15 minutes for a 1/4 degree ocean + # - 7.5 minutes for a 1/8 degree ocean + # - 3.75 minutes for a 1/16 degree ocean + # - 1.875 minutes for a 1/32 degree ocean + + Δt = 30minutes / Δθ + + return all_reduce(min, Δt, arch) +end + const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}} -default_free_surface(grid::TripolarOfSomeKind) = SplitExplicitFreeSurface(grid; substeps=70) + +function default_free_surface(grid::TripolarOfSomeKind; + fixed_Δt = estimate_maximum_Δt(grid), + cfl = 0.7) + free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt) + return free_surface +end + +function default_free_surface(grid::DistributedGrid; + fixed_Δt = compute_maximum_Δt(grid), + cfl = 0.7) + + free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt) + substeps = length(free_surface.substepping.averaging_weights) + substeps = all_reduce(max, substeps, architecture(grid)) + free_surface = SplitExplicitFreeSurface(grid; substeps) + @info "Using a $(free_surface)" + return free_surface +end default_vertical_coordinate(grid) = Oceananigans.Models.ZCoordinate() default_vertical_coordinate(::MutableGridOfSomeKind) = Oceananigans.Models.ZStar() -function default_ocean_closure() +function default_ocean_closure(FT=Oceananigans.defaults.FloatType) mixing_length = CATKEMixingLength(Cᵇ=0.01) turbulent_kinetic_energy_equation = CATKEEquation(Cᵂϵ=1.0) - return CATKEVerticalDiffusivity(; mixing_length, turbulent_kinetic_energy_equation) + + return CATKEVerticalDiffusivity(VerticallyImplicitTimeDiscretization(), FT; mixing_length, turbulent_kinetic_energy_equation) end default_momentum_advection() = VectorInvariant(; vorticity_scheme = WENO(order=9), @@ -79,7 +119,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 = 5minutes, + Δt = estimate_maximum_Δt(grid), closure = default_ocean_closure(), tracers = (:T, :S), free_surface = default_free_surface(grid),