diff --git a/.buildkite/examples_build.yml b/.buildkite/examples_build.yml
new file mode 100644
index 00000000..6f06a800
--- /dev/null
+++ b/.buildkite/examples_build.yml
@@ -0,0 +1,57 @@
+agents:
+  queue: new-central
+  slurm_mem: 8G
+  modules: climacommon/2024_05_27
+  slurm_time: 48:00:00
+
+timeout_in_minutes: 1440
+
+env:
+  JULIA_LOAD_PATH: "${JULIA_LOAD_PATH}:${BUILDKITE_BUILD_CHECKOUT_PATH}/.buildkite"
+  OPENBLAS_NUM_THREADS: 1
+  OMPI_MCA_opal_warn_on_missing_libcuda: 0
+
+steps:
+  - label: "initialize"
+    key: "init"
+    command:
+      - "echo '--- Instantiate project'"
+      - "julia --project -e 'using Pkg; Pkg.instantiate(; verbose=true); Pkg.precompile(; strict=true)'"
+      # force the initialization of the CUDA runtime as it is lazily loaded by default
+      - "julia --project -e 'using CUDA; CUDA.precompile_runtime()'"
+
+    agents:
+      slurm_mem: 120G
+      slurm_gpus: 1
+      slurm_cpus_per_task: 8
+    env:
+      JULIA_NUM_PRECOMPILE_TASKS: 8
+
+  - wait
+
+  - label: "Run documentation"
+    key: "build_documentation"
+    commands:
+      - "julia --color=yes --project=docs/ -e 'using Pkg; Pkg.instantiate(); Pkg.develop(PackageSpec(path=pwd()))'"
+      - "julia --color=yes --project=docs/ docs/make.jl"
+
+    agents:
+      slurm_mem: 120G
+      slurm_gpus: 1
+      slurm_cpus_per_task: 8
+      slurm_ntasks: 1
+      slurm_gpus_per_task: 1
+      slurm_time: 48:00:00
+
+    env:
+      JULIA_NUM_PRECOMPILE_TASKS: 8
+      JULIA_DEBUG: "Documenter"
+      # This environment variable is needed to avoid SSL verification errors when Downloads.jl 
+      # tries to download the bathymetry data. It should not be required so we need to fix our certificates.
+      # and remove this environment variable. ref: https://github.com/JuliaLang/Downloads.jl/issues/97
+      JULIA_SSL_NO_VERIFY: "**" 
+
+    timeout_in_minutes: 1440
+
+  - wait: ~
+    continue_on_failure: true
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index db573661..3c27f771 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -20,6 +20,7 @@ steps:
 
     agents:
       slurm_mem: 32G
+      slurm_gpus: 1
       slurm_cpus_per_task: 8
       slurm_gpus: 1
       slurm_ntasks: 1
diff --git a/Manifest.toml b/Manifest.toml
index aba8f91e..65a8c4e0 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -885,7 +885,9 @@ version = "1.2.0"
 
 [[deps.Oceananigans]]
 deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"]
-git-tree-sha1 = "942c78ad3c95e47bcbab3eeacb0b8a846bcfb33b"
+git-tree-sha1 = "937a5ba25152cbcd0b74e66d21ce931a620a6cc4"
+repo-rev = "main"
+repo-url = "https://github.com/CliMA/Oceananigans.jl.git"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 version = "0.91.4"
 
diff --git a/docs/Project.toml b/docs/Project.toml
index 7331d2d9..c5baf5a7 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,10 +3,9 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
+OrthogonalSphericalShellGrids = "c2be9673-fb75-4747-82dc-aa2bb9f4aed0"
 
 [compat]
 CairoMakie = "0.10.12"
 DataDeps = "0.7"
 Documenter = "1"
-Oceananigans = "0.91.4"
diff --git a/docs/make.jl b/docs/make.jl
index 3b401b64..e551fb70 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -15,15 +15,19 @@ const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
 const OUTPUT_DIR   = joinpath(@__DIR__, "src/literated")
 
 to_be_literated = [
-    "inspect_ecco_data.jl",
+    # "inspect_ecco_data.jl",
+    "generate_bathymetry.jl",
     "generate_surface_fluxes.jl",
-    "single_column_simulation.jl",
-    # "near_global_omip_simulation.jl"
+    # "single_column_simulation.jl",
+    # "mediterranean_simulation_with_ecco_restoring.jl",
+    "near_global_ocean_simulation.jl"
 ]
 
 for file in to_be_literated
     filepath = joinpath(EXAMPLES_DIR, file)
-    Literate.markdown(filepath, OUTPUT_DIR; flavor = Literate.DocumenterFlavor())
+    withenv("JULIA_DEBUG" => "Literate") do
+        Literate.markdown(filepath, OUTPUT_DIR; flavor = Literate.DocumenterFlavor(), execute = true)
+    end
 end
 
 #####
@@ -31,10 +35,10 @@ end
 #####
 
 format = Documenter.HTML(
-       collapselevel = 2,
-      size_threshold = nothing,
-          prettyurls = get(ENV, "CI", nothing) == "true",
-           canonical = "https://clima.github.io/ClimaOceanDocumentation/dev/",
+   collapselevel = 2,
+  size_threshold = nothing,
+      prettyurls = get(ENV, "CI", nothing) == "true",
+       canonical = "https://clima.github.io/ClimaOceanDocumentation/dev/",
 )
 
 pages = [
@@ -48,9 +52,12 @@ pages = [
         ],
 
     "Examples" => [
-        "Inspect ECCO2 data" => "literated/inspect_ecco_data.md",
+        # "Inspect ECCO2 data" => "literated/inspect_ecco_data.md",
+        "Generate bathymetry" => "literated/generate_bathymetry.md",
         "Surface fluxes" => "literated/generate_surface_fluxes.md",
-        "Single column simulation" => "literated/single_column_simulation.md",
+        # "Single column simulation" => "literated/single_column_simulation.md",
+        # "Mediterranean simulation with ECCO restoring" => "literated/mediterranean_simulation_with_ecco_restoring.md",
+        "Near-global Ocean simulation" => "literated/near_global_ocean_simulation.md",
         ]
 ]
 
@@ -86,10 +93,8 @@ for file in files
     rm(file)
 end
 
-withenv("GITHUB_REPOSITORY" => "CliMA/ClimaOceanDocumentation") do
-    deploydocs(        repo = "github.com/CliMA/ClimaOceanDocumentation.git",
-                   versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
-                  forcepush = true,
-                  devbranch = "main",
-               push_preview = true)
-end
+deploydocs(repo = "github.com/CliMA/ClimaOceanDocumentation.git",
+       versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
+      forcepush = true,
+      devbranch = "main",
+   push_preview = true)
diff --git a/examples/freely_decaying_mediterranean.jl b/examples/freely_decaying_mediterranean.jl
deleted file mode 100644
index 2a05e79b..00000000
--- a/examples/freely_decaying_mediterranean.jl
+++ /dev/null
@@ -1,189 +0,0 @@
-using GLMakie
-using Oceananigans
-using Oceananigans: architecture
-using ClimaOcean
-using ClimaOcean.ECCO
-using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity
-using Oceananigans.Coriolis: ActiveCellEnstrophyConserving
-using Oceananigans.Units
-using Printf
-
-#####
-##### Regional Mediterranean grid 
-#####
-
-# Domain and Grid size
-#
-# We construct a grid that represents the Mediterranean sea, 
-# with a resolution of 1/10th of a degree (roughly 10 km resolution)
-λ₁, λ₂  = ( 0, 42) # domain in longitude
-φ₁, φ₂  = (30, 45) # domain in latitude
-# A stretched vertical grid with a Δz of 1.5 meters in the first 50 meters
-z_faces = stretched_vertical_faces(depth = 5000, 
-                             surface_layer_Δz = 2.5, 
-                             stretching = PowerLawStretching(1.070), 
-                             surface_layer_height = 50)
-
-Nx = 15 * 42 # 1 / 15th of a degree resolution
-Ny = 15 * 15 # 1 / 15th of a degree resolution
-Nz = length(z_faces) - 1
-
-grid = LatitudeLongitudeGrid(CPU();
-                             size = (Nx, Ny, Nz),
-                             latitude  = (φ₁, φ₂),
-                             longitude = (λ₁, λ₂),
-                             z = z_faces,
-                             halo = (7, 7, 7))
-
-# Interpolating the bathymetry onto the grid
-#
-# We regrid the bathymetry onto the grid.
-# we allow a minimum depth of 10 meters (all shallower regions are 
-# considered land) and we use 25 intermediate grids (interpolation_passes = 25)
-# Note that more interpolation passes will smooth the bathymetry
-bottom_height = regrid_bathymetry(grid, 
-                                  height_above_water = 1,
-                                  minimum_depth = 10,
-                                  interpolation_passes = 25)
-
-# Let's use an active cell map to elide computation in inactive cells
-grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map = true)
-
-# Correct oceananigans
-import Oceananigans.Advection: nothing_to_default
-
-nothing_to_default(user_value; default) = isnothing(user_value) ? default : user_value
-
-# Constructing the model
-#
-# We construct a model that evolves two tracers, temperature (:T), salinity (:S)
-# We do not have convection since the simulation is just slumping to equilibrium so we do not need a turbulence closure
-# We select a linear equation of state for buoyancy calculation, and WENO schemes for both tracer and momentum advection.
-# The free surface utilizes a split-explicit formulation with a barotropic CFL of 0.75 based on wave speed.
-model = HydrostaticFreeSurfaceModel(; grid,
-                            momentum_advection = WENOVectorInvariant(),
-                              tracer_advection = WENO(grid; order = 7),
-                                  free_surface = SplitExplicitFreeSurface(; cfl = 0.75, grid),
-                                      buoyancy = SeawaterBuoyancy(),
-                                      tracers  = (:T, :S, :c),
-                                      coriolis = HydrostaticSphericalCoriolis(scheme = ActiveCellEnstrophyConserving()))
-
-# Initializing the model
-#
-# the model can be initialized with custom values or with ecco fields.
-# In this case, our ECCO dataset has access to a temperature and a salinity
-# field, so we initialize T and S from ECCO. 
-# We initialize our passive tracer with a surface blob near to the coasts of Libia
-@info "initializing model"
-libia_blob(x, y, z) = z > -20 || (x - 15)^2 + (y - 34)^2 < 1.5 ? 1 : 0
-
-set!(model, T = ECCOMetadata(:temperature), S = ECCOMetadata(:salinity), c = libia_blob)
-
-fig = Figure()
-ax  = Axis(fig[1, 1])
-heatmap!(ax, interior(model.tracers.T, :, :, Nz), colorrange = (10, 20), colormap = :thermal)
-ax  = Axis(fig[1, 2])
-heatmap!(ax, interior(model.tracers.S, :, :, Nz), colorrange = (35, 40), colormap = :haline)
-
-simulation = Simulation(model, Δt = 10minutes, stop_time = 10*365days)
-
-function progress(sim) 
-    u, v, w = sim.model.velocities  
-    T, S = sim.model.tracers
-
-    @info @sprintf("Time: %s, Iteration %d, Δt %s, max(vel): (%.2e, %.2e, %.2e), max(T, S): %.2f, %.2f\n",
-                   prettytime(sim.model.clock.time),
-                   sim.model.clock.iteration,
-                   prettytime(sim.Δt),
-                   maximum(abs, u), maximum(abs, v), maximum(abs, w),
-                   maximum(abs, T), maximum(abs, S))
-end
-
-simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))
-
-# Simulation warm up!
-#
-# We have regridded from a coarse solution (1/4er of a degree) to a
-# fine grid (1/15th of a degree). Also, the bathymetry has little mismatches 
-# that might crash our simulation. We warm up the simulation with a little 
-# time step for few iterations to allow the solution to adjust to the new_grid
-# bathymetry
-simulation.Δt = 10
-simulation.stop_iteration = 1000
-run!(simulation)
-
-# Run the real simulation
-#
-# Now that the solution has adjusted to the bathymetry we can ramp up the time
-# step size. We use a `TimeStepWizard` to automatically adapt to a cfl of 0.2
-wizard = TimeStepWizard(; cfl = 0.2, max_Δt = 10minutes, max_change = 1.1)
-
-simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))
-
-# Let's reset the maximum number of iterations
-simulation.stop_iteration = Inf
-
-simulation.output_writers[:surface_fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers);
-                                                              indices = (:, :, Nz),
-                                                              schedule = TimeInterval(1day),
-                                                              overwrite_existing = true,
-                                                              filename = "med_surface_field")
-
-run!(simulation)
-
-# Record a video
-#
-# Let's read the data and record a video of the Mediterranean Sea's surface
-# (1) Zonal velocity (u)
-# (2) Meridional velocity (v)
-# (3) Temperature (T)
-# (4) Salinity (S)
-u_series = FieldTimeSeries("med_surface_field.jld2", "u")
-v_series = FieldTimeSeries("med_surface_field.jld2", "v")
-T_series = FieldTimeSeries("med_surface_field.jld2", "T")
-S_series = FieldTimeSeries("med_surface_field.jld2", "S")
-c_series = FieldTimeSeries("med_surface_field.jld2", "c")
-iter = Observable(1)
-
-u = @lift begin
-    f = interior(u_series[$iter], :, :, 1)
-    f[f .== 0] .= NaN
-    f
-end
-v = @lift begin
-    f = interior(v_series[$iter], :, :, 1)
-    f[f .== 0] .= NaN
-    f
-end
-T = @lift begin
-    f = interior(T_series[$iter], :, :, 1)
-    f[f .== 0] .= NaN
-    f
-end
-S = @lift begin
-    f = interior(S_series[$iter], :, :, 1)
-    f[f .== 0] .= NaN
-    f
-end
-c = @lift begin
-    f = interior(c_series[$iter], :, :, 1)
-    f[f .== 0] .= NaN
-    f
-end
-
-fig = Figure()
-ax  = Axis(fig[1, 1], title = "surface zonal velocity ms⁻¹")
-heatmap!(u)
-ax  = Axis(fig[1, 2], title = "surface meridional velocity ms⁻¹")
-heatmap!(v)
-ax  = Axis(fig[2, 1], title = "surface temperature ᵒC")
-heatmap!(T)
-ax  = Axis(fig[2, 2], title = "surface salinity psu")
-heatmap!(S)
-ax  = Axis(fig[2, 3], title = "passive tracer -")
-heatmap!(c)
-
-GLMakie.record(fig, "mediterranean_video.mp4", 1:length(u_series.times); framerate = 5) do i
-    @info "recording iteration $i"
-    iter[] = i    
-end
\ No newline at end of file
diff --git a/examples/generate_bathymetry.jl b/examples/generate_bathymetry.jl
index 4db9c54b..fcf3652f 100644
--- a/examples/generate_bathymetry.jl
+++ b/examples/generate_bathymetry.jl
@@ -1,70 +1,78 @@
-using GLMakie
+# # Generate bathymetry data for the Mediterranean Sea
+#
+# This script shows how to configure an Immersed boundary grid with realistic bathymetry using ClimaOcean.jl
+# by generating the bathymetry data for the Mediterranean Sea.
+#
+# For this example, we need Oceananigans for the LatitudeLongitudeGrid and Field utilities, 
+# ClimaOcean to donwload and regrid the bathymetry, and CairoMakie to visualize the grid.
+
+using ClimaOcean
 using Oceananigans
-using ClimaOcean.Bathymetry: regrid_bathymetry
+using CairoMakie
+
+# We start by defining a gridded domain for the Mediterranean Sea using the `LatitudeLongitudeGrid` from Oceananigans.
+#
+# The Mediterranean sea is positioned roughly between 28ᵒ and 48ᵒ latitude and 0ᵒ and 42ᵒ longitude.
+# We define a grid in this region and to have a reasonable resolution, we set a grid resolution to 1/25ᵒ
+# in both latitude and longitude directions.
+
+latitude_range = (28, 48)
+longitude_range = (0, 42)
+
+Nφ = 25 * (latitude_range[2] - latitude_range[1])
+Nλ = 25 * (longitude_range[2] - longitude_range[1])
+
+grid = LatitudeLongitudeGrid(size = (Nλ, Nφ, 1),
+                             latitude = latitude_range,
+                             longitude = longitude_range,
+                             z = (0, 1))
+
+# Next, we generate the bathymetry data for the Mediterranean Sea using the
+# `regrid_bathymetry` function from ClimaOcean. The function downloads the bathymetry
+# data from the ETOPO1 dataset, regrids it to the provided grid, and returns the
+# bathymetry field. The three different regidding procedures below demonstrate the effect
+# of different parameters on the generated bathymetry:
+#
+# - `h_rough`  shows the output of the function with default parameters, which means only
+#   one interpolation passes and no restrictions on connected regions.
+# - `h_smooth` shows the output of the function with 40 interpolation passes, which results
+#    in a smoother bathymetry.
+# - `h_no_connected_regions` shows the output of the function with `connected_regions_allowed = 0`, which
+#    means that the function does not allow connected regions in the bathymetry  (e.g., lakes)
+#    and fills them with land.
+
+h_rough = regrid_bathymetry(grid)
+h_smooth = regrid_bathymetry(grid; interpolation_passes = 40)
+h_no_connected_regions = regrid_bathymetry(grid; connected_regions_allowed = 0)
+nothing # hide
+
+# Finally, we visualize the generated bathymetry data for the Mediterranean Sea using CairoMakie.
 
-#####
-##### Quarter degree near-global grid
-#####
-
-grid = LatitudeLongitudeGrid(CPU();
-                             size = (4 * 360, 4 * 160, 1),
-                             latitude = (-80, 80),
-                             longitude = (-180, 180),
-                             z = (0, 1),
-                             halo = (4, 4, 4))
-
-h = regrid_bathymetry(grid, height_above_water=1, minimum_depth=5)
-
-λ, φ, z = nodes(h)
+λ, φ, z = nodes(h_smooth)
 
-land = interior(h) .>= 0
-interior(h)[land] .= NaN
+land_smooth = interior(h_smooth) .>= 0
+interior(h_smooth)[land_smooth] .= NaN
 
-fig = Figure(size=(2400, 1200))
-ax = Axis(fig[1, 1])
-heatmap!(ax, λ, φ, interior(h, :, :, 1), nan_color=:white, colorrange=(-5000, 0))
+land_rough = interior(h_rough) .>= 0
+interior(h_rough)[land_rough] .= NaN
 
-λp = -112.45
-φp = 42.86
-text = "😻"
-text!(ax, λp, φp; text, fontsize=30)
+land_no_connected_regions = interior(h_no_connected_regions) .>= 0
+interior(h_no_connected_regions)[land_no_connected_regions] .= NaN
 
-display(fig)
+fig = Figure(resolution=(850, 1150))
 
-#####
-##### Regional Mediterranean grid 
-#####
+ax = Axis(fig[1, 1], title = "Rough bathymetry", xlabel = "Longitude", ylabel = "Latitude")
+hm = heatmap!(ax, λ, φ, - interior(h_rough, :, :, 1), nan_color=:white, colormap = Reverse(:deep))
 
-# 1/25th degree resolution
-Nλ = 25 * 55
-Nφ = 25 * 25
+ax = Axis(fig[2, 1], title = "Smooth bathymetry", xlabel = "Longitude", ylabel = "Latitude")
+hm = heatmap!(ax, λ, φ, - interior(h_smooth, :, :, 1), nan_color=:white, colormap = Reverse(:deep))
 
-grid = LatitudeLongitudeGrid(CPU();
-                             size = (Nλ, Nφ, 1),
-                             latitude = (25, 50),
-                             longitude = (-10, 45),
-                             z = (0, 1),
-                             halo = (4, 4, 4))
+ax = Axis(fig[3, 1], title = "Bathymetry without connected regions}", xlabel = "Longitude", ylabel = "Latitude")
+hm = heatmap!(ax, λ, φ, - interior(h_no_connected_regions, :, :, 1), nan_color=:white, colormap = Reverse(:deep))
 
-h_smooth  = regrid_bathymetry(grid, height_above_water=1, minimum_depth=10, interpolation_passes = 40)
-h_rough   = regrid_bathymetry(grid, height_above_water=1, minimum_depth=10, interpolation_passes = 1)
-h_nolakes = regrid_bathymetry(grid, height_above_water=1, minimum_depth=10, connected_regions_allowed = 0)
+cb = Colorbar(fig[1:3, 2], hm, height = Relative(3/4), label = "Depth (m)")
 
-λ, φ, z = nodes(h_smooth)
+save("different_bottom_heights.png", fig)
+nothing #hide
 
-land_smooth = interior(h_smooth) .>= 0
-interior(h_smooth)[land_smooth] .= NaN
-land_rough = interior(h_rough) .>= 0
-interior(h_rough)[land_rough] .= NaN
-land_nolakes = interior(h_nolakes) .>= 0
-interior(h_nolakes)[land_nolakes] .= NaN
-
-fig = Figure(resolution=(2400, 800))
-ax = Axis(fig[1, 1])
-heatmap!(ax, λ, φ, interior(h_smooth,  :, :, 1), nan_color=:white) #, colorrange=(-5000, 0))
-ax = Axis(fig[1, 2])
-heatmap!(ax, λ, φ, interior(h_rough,   :, :, 1), nan_color=:white) #, colorrange=(-5000, 0))
-ax = Axis(fig[1, 3])
-heatmap!(ax, λ, φ, interior(h_nolakes, :, :, 1), nan_color=:white) #, colorrange=(-5000, 0))
-
-display(fig)
+# ![](different_bottom_heights.png)
diff --git a/examples/generate_surface_fluxes.jl b/examples/generate_surface_fluxes.jl
index 1f54ba95..078dbf58 100644
--- a/examples/generate_surface_fluxes.jl
+++ b/examples/generate_surface_fluxes.jl
@@ -1,4 +1,4 @@
-# Calculating surface fluxes with an ocean and an atmosphere
+# # Surface fluxes from prescribed ocean and atmosphere
 #
 # ClimaOcean uses bulk formulae to estimate the surface exchange of momentum,
 # heat, and water vapor between the atmosphere and the ocean.
@@ -37,6 +37,7 @@ nothing #hide
 # ![](ecco_continents.png)
 
 # Next, we construct our atmosphere and ocean.
+#
 # The atmosphere is prescribed, downloaded from the JRA55 dataset.
 # It contains:
 # - zonal wind `u`
@@ -67,7 +68,7 @@ ocean = ocean_simulation(grid; momentum_advection = nothing,
 # our ocean with initial conditions. To do this, we can use the ECCO2 dataset by
 # `set!`ting the model with the `ECCOMetadata`. If no date is specified,
 # the fields corresponding to January 1st, 1992 (the first available date in
-# ECCO2) are used.
+# ECCO2 dataset) are used.
 # This command will download the fields to the local machine.
 
 set!(ocean.model;
@@ -91,13 +92,13 @@ coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
 # Now that the surface fluxes are computed, we can extract and visualize them.
 # The turbulent fluxes are stored in `coupled_model.fluxes.turbulent`.
 # 
-# Qs = coupled_model.fluxes.turbulent.fields.sensible_heat : the sensible heat flux 
-# Ql = coupled_model.fluxes.turbulent.fields.latent_heat   : the latent heat flux 
-# τx = coupled_model.fluxes.turbulent.fields.x_momentum    : the zonal wind stress
-# τy = coupled_model.fluxes.turbulent.fields.y_momentum    : the meridional wind stress
-# Mv = coupled_model.fluxes.turbulent.fields.water_vapor   : evaporation
+# Qs = coupled_model.fluxes.turbulent.fields.sensible_heat : the sensible heat flux (in Wm⁻²)
+# Ql = coupled_model.fluxes.turbulent.fields.latent_heat   : the latent heat flux  (in Wm⁻²)
+# τx = coupled_model.fluxes.turbulent.fields.x_momentum    : the zonal wind stress (in Nm)
+# τy = coupled_model.fluxes.turbulent.fields.y_momentum    : the meridional wind stress (in Nm)
+# Mv = coupled_model.fluxes.turbulent.fields.water_vapor   : evaporation (in kg m⁻²s⁻¹)
 #
-# They are 3D fields with one point in the vertical. To extract the data, we use the 
+# They are 2D fields (3D data structures with one point in the vertical). To extract the data, we use the 
 # `interior` functionality from Oceananigans.
 
 turbulent_fluxes = coupled_model.fluxes.turbulent.fields
@@ -107,26 +108,45 @@ Ql = interior(turbulent_fluxes.latent_heat,   :, :, 1)
 τx = interior(turbulent_fluxes.x_momentum,    :, :, 1)
 τy = interior(turbulent_fluxes.y_momentum,    :, :, 1)
 Mv = interior(turbulent_fluxes.water_vapor,   :, :, 1)
-nothing
-
-fig = Figure()
+nothing #hide
 
-ax = Axis(fig[1, 1], title = "Sensible heat flux")
-heatmap!(ax, Qs; colormap = :bwr)
+fig = Figure(size = (800, 400))
+ax = Axis(fig[1, 1], title = "Sensible heat flux (Wm⁻²)")
+hm = heatmap!(ax, Qs; colormap = :bwr)
+hidedecorations!(ax)
+save("sensible_heat_flux.png", fig)
+nothing #hide
+# ![](sensible_heat_flux.png)
 
-ax = Axis(fig[1, 2], title = "Latent heat flux")
+fig = Figure(size = (800, 400))
+ax = Axis(fig[1, 2], title = "Latent heat flux (Wm⁻²)")
 heatmap!(ax, Ql; colormap = :bwr)
+hidedecorations!(ax)
+save("latent_heat_flux.png", fig)
+nothing #hide
+# ![](latent_heat_flux.png)
 
-ax = Axis(fig[2, 1], title = "Zonal wind stress")
+fig = Figure(size = (800, 400))
+ax = Axis(fig[2, 1], title = "Zonal wind stress (Nm)")
 heatmap!(ax, τx; colormap = :bwr)
+hidedecorations!(ax)
+save("zonal_wind_stress.png", fig)
+nothing #hide
+# ![](zonal_wind_stress.png)
 
-ax = Axis(fig[2, 2], title = "Meridional wind stress")
+fig = Figure(size = (800, 400))
+ax = Axis(fig[2, 2], title = "Meridional wind stress (Nm)")
 heatmap!(ax, τy; colormap = :bwr)
+hidedecorations!(ax)
+save("meridional_wind_stress.png", fig)
+nothing #hide
+# ![](meridional_wind_stress.png)
 
-ax = Axis(fig[3, 1], title = "Evaporation")
+fig = Figure(size = (800, 400))
+ax = Axis(fig[3, 1], title = "Water vapor flux (kg m⁻²s⁻¹)")
 heatmap!(ax, Mv; colormap = :bwr)
-
-save("turbulent_fluxes.png", fig)
+hidedecorations!(ax)
+save("water_vapor_flux.png", fig)
 nothing #hide
+# ![](water_vapor_flux.png)
 
-# ![](turbulent_fluxes.png)
diff --git a/examples/inspect_JRA55_data.jl b/examples/inspect_JRA55_data.jl
index 37403da3..23445799 100644
--- a/examples/inspect_JRA55_data.jl
+++ b/examples/inspect_JRA55_data.jl
@@ -1,5 +1,5 @@
 using ClimaOcean
-using GLMakie
+using CairoMakie
 using Oceananigans
 using Oceananigans.Units
 using Printf
diff --git a/examples/inspect_ecco_data.jl b/examples/inspect_ecco_data.jl
index 71ec6609..9ee3ab69 100644
--- a/examples/inspect_ecco_data.jl
+++ b/examples/inspect_ecco_data.jl
@@ -76,7 +76,7 @@ text!(axS, 50, 50, text=depth_str, justification=:center, fontsize=24)
 stillframes = 10
 movingframes = Nz
 
-record(fig, "ECCO_temperature_salinity.gif", framerate=4) do io
+record(fig, "ECCO_temperature_salinity.mp4", framerate=4) do io
 
     [recordframe!(io) for _ = 1:stillframes]
 
@@ -94,4 +94,6 @@ record(fig, "ECCO_temperature_salinity.gif", framerate=4) do io
 
     [recordframe!(io) for _ = 1:stillframes]
 end
+nothing #hide
 
+# ![](ECCO_temperature_salinity.mp4)
diff --git a/examples/mediterranean_simulation_with_ecco_restoring.jl b/examples/mediterranean_simulation_with_ecco_restoring.jl
index c5e37945..6a33fb3c 100644
--- a/examples/mediterranean_simulation_with_ecco_restoring.jl
+++ b/examples/mediterranean_simulation_with_ecco_restoring.jl
@@ -9,12 +9,12 @@
 
 # ## Initial Setup with Package Imports
 #
-# The script begins by importing necessary Julia packages for visualization (GLMakie), 
+# The script begins by importing necessary Julia packages for visualization (CairoMakie), 
 # ocean modeling (Oceananigans, ClimaOcean), and handling of dates and times (CFTime, Dates). 
 # These packages provide the foundational tools for creating the simulation environment, 
 # including grid setup, physical processes modeling, and data visualization.
 
-using GLMakie
+using CairoMakie
 using Oceananigans
 using Oceananigans: architecture
 using ClimaOcean
@@ -45,7 +45,7 @@ Nx = 15 * 42 # 1 / 15th of a degree resolution
 Ny = 15 * 15 # 1 / 15th of a degree resolution
 Nz = length(z_faces) - 1
 
-grid = LatitudeLongitudeGrid(CPU();
+grid = LatitudeLongitudeGrid(GPU();
                              size = (Nx, Ny, Nz),
                              latitude  = (φ₁, φ₂),
                              longitude = (λ₁, λ₂),
@@ -64,7 +64,7 @@ bottom_height = regrid_bathymetry(grid,
                                   interpolation_passes = 25,
                                   connected_regions_allowed = 1)
 
-grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map = true)
+grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
 
 # ## Downloading ECCO data
 #
@@ -77,8 +77,8 @@ dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGre
 temperature = ECCOMetadata(:temperature, dates, ECCO4Monthly())
 salinity    = ECCOMetadata(:salinity,    dates, ECCO4Monthly())
 
-FT = ECCO_restoring_forcing(temperature; timescale = 2days)
-FS = ECCO_restoring_forcing(salinity;    timescale = 2days)
+FT = ECCO_restoring_forcing(temperature; architecture = GPU(), timescale = 2days)
+FS = ECCO_restoring_forcing(salinity;    architecture = GPU(), timescale = 2days)
 
 # Constructing the Simulation
 #
@@ -134,10 +134,11 @@ run!(ocean)
 
 wizard = TimeStepWizard(; cfl = 0.2, max_Δt = 10minutes, max_change = 1.1)
 
-coean.callbacks[:wizard] = Callback(wizard, IterationInterval(10))
+ocean.callbacks[:wizard] = Callback(wizard, IterationInterval(10))
 
 # Let's reset the maximum number of iterations
 ocean.stop_iteration = Inf
+ocean.stop_time = 200days
 
 ocean.output_writers[:surface_fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers);
                                                          indices = (:, :, Nz),
@@ -200,7 +201,7 @@ heatmap!(S)
 ax  = Axis(fig[2, 3], title = "passive tracer -")
 heatmap!(c)
 
-GLMakie.record(fig, "mediterranean_video.mp4", 1:length(u_series.times); framerate = 5) do i
+CairoMakie.record(fig, "mediterranean_video.mp4", 1:length(u_series.times); framerate = 5) do i
     @info "recording iteration $i"
     iter[] = i    
 end
\ No newline at end of file
diff --git a/examples/near_global_ocean_simulation.jl b/examples/near_global_ocean_simulation.jl
new file mode 100644
index 00000000..12df9c99
--- /dev/null
+++ b/examples/near_global_ocean_simulation.jl
@@ -0,0 +1,290 @@
+# # Near-global ocean simulation
+#
+# This Julia script sets up and runs a near-global ocean simulation using the Oceananigans.jl and ClimaOcean.jl packages. 
+# The simulation covers latitudes from 75°S to 75°N with a horizontal resolution of 1/4 degree and 40 vertical levels. 
+#
+# The simulation runs for one year, and the results are visualized using the CairoMakie.jl package.
+#
+# ## Initial setup with package imports
+#
+# The script begins by importing the necessary Julia packages for visualization (CairoMakie), 
+# ocean modeling (Oceananigans, ClimaOcean), and handling dates and times (CFTime, Dates). 
+# These packages provide the foundational tools for setting up the simulation environment, 
+# including grid setup, physical processes modeling, and data visualization.
+
+using Printf
+using Oceananigans
+using Oceananigans.Units
+using ClimaOcean
+using CairoMakie
+
+using CFTime
+using Dates
+
+# ### Grid configuration 
+#
+# We define a global grid with a horizontal resolution of 1/4 degree and 40 vertical levels. 
+# The grid is a `LatitudeLongitudeGrid` capped at 75°S to 75°N.
+# We use an exponential vertical spacing to better resolve the upper ocean layers. The total depth of the domain is set to 6000 meters.
+# Finally, we specify the architecture for the simulation, which in this case is a GPU.
+
+arch = GPU() 
+
+z_faces = exponential_z_faces(Nz=40, depth=6000)
+
+Nx = 1440
+Ny = 600
+Nz = length(z_faces) - 1
+
+grid = LatitudeLongitudeGrid(arch;
+                             size = (Nx, Ny, Nz),
+                             halo = (7, 7, 7),
+                             z = z_faces, 
+                             latitude  = (-75, 75),
+                             longitude = (0, 360))
+
+# ### Bathymetry and immersed boundary
+#
+# We retrieve the bathymetry from the ETOPO1 data, ensuring a minimum depth of 10 meters
+# (depths shallower than this are considered land). The `interpolation_passes` parameter
+# specifies the number of passes to interpolate the bathymetry data. A larger number
+# results in a smoother bathymetry. We also remove all connected regions (such as inland
+# lakes) from the bathymetry data by specifying `connected_regions_allowed = 2` (Mediterrean
+# sea an North sea in addition to the ocean).
+
+bottom_height = regrid_bathymetry(grid; 
+                                  minimum_depth = 10,
+                                  interpolation_passes = 5,
+                                  connected_regions_allowed = 2)
+ 
+grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) 
+
+bathymetry = deepcopy(Array(interior(bottom_height, :, :, 1)))
+bathymetry[bathymetry .>= 0] .= NaN
+
+fig = Figure(size = (1200, 400))
+ax  = Axis(fig[1, 1])
+hm = heatmap!(ax, bathymetry, colormap = :deep, colorrange = (-6000, 0))
+cb = Colorbar(fig[0, 1], hm, label = "Bottom height (m)", vertical = false)
+hidedecorations!(ax)
+
+save("bathymetry.png", fig)
+nothing #hide
+
+# ![](bathymetry.png)
+
+# ### Ocean model configuration
+#
+# To configure the ocean simulation, we use the `ocean_simulation` function from ClimaOcean.jl. This function allows us to build
+# an ocean simulation with default parameters and numerics. The defaults include:
+# - CATKE turbulence closure for vertical mixing, see [`CATKEVerticalDiffusivity`](@ref)
+# - WENO-based advection scheme for momentum in the vector invariant form, see [`WENOVectorInvariant`](@ref)
+# - WENO-based advection scheme for tracers, see [`WENO`](@ref)
+# - `SplitExplicitFreeSurfaceSolver` with 75 substeps, see [`SplitExplicitFreeSurface`](@ref)
+# - TEOS-10 equation of state, see [`TEOS10EquationOfState`](@ref)
+# - Quadratic bottom drag with a drag coefficient of 0.003
+#
+# The ocean model is then initialized with the ECCO2 temperature and salinity fields for January 1, 1993.
+
+ocean = ocean_simulation(grid)
+model = ocean.model
+
+date  = DateTimeProlepticGregorian(1993, 1, 1)
+
+set!(model, 
+     T = ECCOMetadata(:temperature; date),
+     S = ECCOMetadata(:salinity;    date))
+nothing #hide
+
+# ### Prescribed atmosphere and radiation
+#
+# The atmospheric data is prescribed using the JRA55 dataset, which is loaded
+# into memory in 4 snapshots at a time. The JRA55 dataset provides atmospheric
+# data such as temperature, humidity, and wind fields to calculate turbulent fluxes
+# using bulk formulae, see [`CrossRealmFluxes`](@ref).
+#
+# The radiation model specifies an ocean albedo emissivity to compute the net radiative
+# fluxes. The default ocean albedo is based on Payne (1982) and depends on cloud cover
+# (calculated from the ratio of maximum possible incident solar radiation to actual
+# incident solar radiation) and latitude. The ocean emissivity is set to 0.97.
+
+backend    = JRA55NetCDFBackend(41) 
+atmosphere = JRA55_prescribed_atmosphere(arch; backend)
+radiation  = Radiation(arch)
+nothing #hide
+
+# ### Sea ice model 
+#
+# This simulation includes a simplified representation of ice cover where the
+# air-sea fluxes are shut down whenever the sea surface temperature is below
+# the freezing point. Only heating fluxes are allowed. This is not a full
+# sea ice model, but it prevents the temperature from dropping excessively
+# low by including atmosphere-ocean fluxes.
+
+sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce()
+nothing #hide
+
+# ## The coupled simulation
+#
+# Finally, we define the coupled model, which includes the ocean, atmosphere,
+# and radiation parameters. The model is constructed using the `OceanSeaIceModel`
+# constructor.
+#
+# We then create a coupled simulation, starting with a time step of 10 seconds
+# and running the simulation for 10 days.
+# We will eventually increase the time step size and end time as the simulation
+# progresses and initialization shocks dissipate.
+#
+# We also define a callback function to monitor the simulation's progress.
+# This function prints the current time, iteration, time step,
+# as well as the maximum velocities and tracers in the domain. The wall time
+# is also printed to monitor the time taken for each iteration.
+
+coupled_model      = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
+coupled_simulation = Simulation(coupled_model; Δt=10, stop_time = 10days)
+
+wall_time = [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[1])
+
+    @info @sprintf("Time: %s, Iteration %d, Δt %s, max(vel): (%.2e, %.2e, %.2e), max(T): %.2f, min(T): %.2f, wtime: %s \n",
+                   prettytime(ocean.model.clock.time),
+                   ocean.model.clock.iteration,
+                   prettytime(ocean.Δt),
+                   umax..., Tmax, Tmin, prettytime(step_time))
+
+     wall_time[1] = time_ns()
+end
+
+coupled_simulation.callbacks[:progress] = Callback(progress, IterationInterval(1000))
+nothing #hide
+
+# ### Set up output writers
+#
+# We define output writers to save the simulation data at regular intervals.
+# In this case, we save the surface fluxes and surface fields at a relatively high frequency (every day).
+# The `indices` keyword argument allows us to save down a slice at the surface, which is located at `k = grid.Nz`
+
+ocean.output_writers[:surface] = JLD2OutputWriter(model, merge(model.tracers, model.velocities);
+                                                  schedule = TimeInterval(1days),
+                                                  filename = "surface",
+                                                  indices = (:, :, grid.Nz),
+                                                  overwrite_existing = true,
+                                                  array_type = Array{Float32})
+nothing #hide
+
+# ### Spinning up the simulation
+#
+# As an initial condition, we have interpolated ECCO tracer fields onto our custom grid.
+# The bathymetry of the original ECCO data may differ from our grid, so the initialization of the velocity
+# field might cause shocks if a large time step is used.
+#
+# Therefore, we spin up the simulation with a small time step to ensure that the interpolated initial
+# conditions adapt to the model numerics and parameterization without causing instability. A 10-day
+# integration with a maximum time step of 1.5 minutes should be sufficient to dissipate spurious
+# initialization shocks.
+# We use an adaptive time step that maintains the [CFL condition](https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition) equal to 0.1.
+# For this scope, we use the Oceananigans utility `conjure_time_step_wizard!` (see Oceanigans's documentation).
+
+ocean.stop_time = 10days
+conjure_time_step_wizard!(ocean; cfl = 0.1, max_Δt = 90, max_change = 1.1)
+run!(coupled_simulation)
+nothing #hide
+
+# ### Running the simulation
+#
+# Now that the simulation has spun up, we can run it for the full 100 days.
+# We increase the maximum time step size to 10 minutes and let the simulation run for 100 days.
+# This time, we set the CFL in the time_step_wizard to be 0.25 as this is the maximum recommended CFL to be
+# used in conjunction with Oceananigans' hydrostatic time-stepping algorithm ([two step Adams-Bashfort](https://en.wikipedia.org/wiki/Linear_multistep_method))
+
+ocean.stop_time = 100days
+coupled_simulation.stop_time = 100days
+conjure_time_step_wizard!(ocean; cfl = 0.25, max_Δt = 10minutes, max_change = 1.1)
+run!(coupled_simulation)
+nothing #hide
+
+# ## Visualizing the results
+# 
+# The simulation has finished, let's visualize the results.
+# In this section we pull up the saved data and create visualizations using the CairoMakie.jl package.
+# In particular, we generate an animation of the evolution of surface fields:
+# surface speed (s), surface temperature (T), and turbulent kinetic energy (e).
+
+u = FieldTimeSeries("surface.jld2", "u"; backend = OnDisk())
+v = FieldTimeSeries("surface.jld2", "v"; backend = OnDisk())
+T = FieldTimeSeries("surface.jld2", "T"; backend = OnDisk())
+e = FieldTimeSeries("surface.jld2", "e"; backend = OnDisk())
+
+times = u.times
+Nt = length(times)
+
+iter = Observable(Nt)
+
+Ti = @lift begin
+     Ti = interior(T[$iter], :, :, 1)
+     Ti[Ti .== 0] .= NaN
+     Ti
+end
+
+ei = @lift begin
+     ei = interior(e[$iter], :, :, 1)
+     ei[ei .== 0] .= NaN
+     ei
+end
+
+si = @lift begin
+     s = Field(sqrt(u[$iter]^2 + v[$iter]^2))
+     compute!(s)
+     s = interior(s, :, :, 1)
+     s[s .== 0] .= NaN
+     s
+end
+
+
+fig = Figure(size = (800, 400))
+ax = Axis(fig[1, 1])
+hm = heatmap!(ax, si, colorrange = (0, 0.5), colormap = :deep)
+cb = Colorbar(fig[0, 1], hm, vertical = false, label = "Surface speed (ms⁻¹)")
+hidedecorations!(ax)
+
+CairoMakie.record(fig, "near_global_ocean_surface_s.mp4", 1:Nt, framerate = 8) do i
+    iter[] = i
+end
+nothing #hide
+ 
+ # ![](near_global_ocean_surface_s.mp4)
+ 
+fig = Figure(size = (800, 400))
+ax = Axis(fig[1, 1])
+hm = heatmap!(ax, Ti, colorrange = (-1, 30), colormap = :magma)
+cb = Colorbar(fig[0, 1], hm, vertical = false, label = "Surface Temperature (Cᵒ)")
+hidedecorations!(ax)
+
+CairoMakie.record(fig, "near_global_ocean_surface_T.mp4", 1:Nt, framerate = 8) do i
+    iter[] = i
+end
+nothing #hide
+ 
+# ![](near_global_ocean_surface_T.mp4)
+
+fig = Figure(size = (800, 400))
+ax = Axis(fig[1, 1])
+hm = heatmap!(ax, ei, colorrange = (0, 1e-3), colormap = :solar)
+cb = Colorbar(fig[0, 1], hm, vertical = false, label = "Turbulent Kinetic Energy (m²s⁻²)")
+hidedecorations!(ax)
+
+CairoMakie.record(fig, "near_global_ocean_surface_e.mp4", 1:Nt, framerate = 8) do i
+    iter[] = i
+end
+nothing #hide
+
+# ![](near_global_ocean_surface_e.mp4)
diff --git a/examples/near_global_omip_simulation.jl b/examples/near_global_omip_simulation.jl
deleted file mode 100644
index eca17321..00000000
--- a/examples/near_global_omip_simulation.jl
+++ /dev/null
@@ -1,195 +0,0 @@
-using Printf
-using Oceananigans
-using Oceananigans.Units
-using Oceananigans: architecture
-using ClimaOcean
-using ClimaOcean.ECCO
-using ClimaOcean.OceanSimulations
-using ClimaOcean.OceanSeaIceModels
-
-using CFTime
-using Dates
-
-#####
-##### Near - Global Ocean at 1/4th of a degree
-#####
-
-# 40 vertical levels
-z_faces = exponential_z_faces(Nz=40, depth=6000)
-
-Nx = 1440
-Ny = 600
-Nz = length(z_faces) - 1
-
-# Running on a GPU
-arch = GPU() 
-
-# A near-global grid from 75ᵒ S to 75ᵒ N
-grid = LatitudeLongitudeGrid(arch; 
-                             size = (Nx, Ny, Nz), 
-                             halo = (7, 7, 7), 
-                             z = z_faces, 
-                             longitude = (0, 360),
-                             latitude = (-75, 75))
-
-# We retrieve the bathymetry from the ETOPO1 data by ensuring a 
-# minimum depth of 10 meters (everyhting shallower is considered land)
-# and removing all connected regions (inland lakes)
-bottom_height = retrieve_bathymetry(grid; 
-                                    minimum_depth = 10,
-                                    dir = "./",
-                                    interpolation_passes = 20,
-                                    connected_regions_allowed = 0)
- 
-# An immersed boundary using a staircase representation of bathymetry
-grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map = true) 
-
-#####
-##### The Ocean component
-#####                             
-
-# We retain all the defaults for the ocean model
-ocean = ocean_simulation(grid) 
-model = ocean.model
-
-date  = DateTimeProlepticGregorian(1993, 1, 1)
-
-# We interpolate the initial conditions from the ECCO2 dataset 
-# (for the moment these are both 1st January 1993)
-set!(model, 
-     T = ECCO2Metadata(:temperature, date, ECCO2Daily()),
-     S = ECCO2Metadata(:salinity,    date, ECCO2Daily()))
-
-#####
-##### The atmosphere
-#####
-
-# The whole prescribed atmosphere is loaded in memory 
-# 4 snapshots at the time.
-backend    = JRA55NetCDFBackend(4) 
-atmosphere = JRA55_prescribed_atmosphere(arch; backend)
-
-# Tabulated ocean albedo from Payne (1982)
-# ocean emissivity is the default 0.97
-radiation  = Radiation(arch)
-
-#####
-##### The atmospheric-forced coupled ocean-seaice model
-#####
-
-# Simplistic sea ice that ensure "no-cooling-fluxes" where `T < T_minimum`
-# to change with a thermodynamic sea-ice model
-sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce()
-
-# The complete coupled model
-coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
-
-wall_time = [time_ns()]
-
-# We define a progress function that
-# shows the maximum values of velocity and temperature
-# to make sure everything proceedes as planned
-function progress(sim) 
-    u, v, w = sim.model.velocities  
-    T = sim.model.tracers.T
-
-    Tmax = maximum(interior(T))
-    Tmin = minimum(interior(T))
-    umax = maximum(interior(u)), maximum(interior(v)), maximum(interior(w))
-    step_time = 1e-9 * (time_ns() - wall_time[1])
-
-    @info @sprintf("Time: %s, Iteration %d, Δt %s, max(vel): (%.2e, %.2e, %.2e), max(trac): %.2f, %.2f, wtime: %s \n",
-                   prettytime(sim.model.clock.time),
-                   sim.model.clock.iteration,
-                   prettytime(sim.Δt),
-                   umax..., Tmax, Tmin, prettytime(step_time))
-
-     wall_time[1] = time_ns()
-end
-
-ocean.callbacks[:progress] = Callback(progress, IterationInterval(10))
-
-fluxes = (u = model.velocities.u.boundary_conditions.top.condition,
-          v = model.velocities.v.boundary_conditions.top.condition,
-          T = model.tracers.T.boundary_conditions.top.condition,
-          S = model.tracers.S.boundary_conditions.top.condition)
-
-
-output_kwargs = (; overwrite_existing = true, array_type = Array{Float32})
-
-# We add a couple of outputs: the surface state every 12 hours, the surface fluxes
-# every 12 hours, and the whole state every ten days
-ocean.output_writers[:fluxes] = JLD2OutputWriter(model, fluxes;
-                                                  schedule = TimeInterval(0.5days),
-                                                  overwrite_existing = true,
-                                                  filename = "surface_fluxes",
-                                                  output_kwargs...)
-
-ocean.output_writers[:surface] = JLD2OutputWriter(model, merge(model.tracers, model.velocities);
-                                                  schedule = TimeInterval(0.5days),
-                                                  filename = "surface",
-                                                  indices = (:, :, grid.Nz),
-                                                  output_kwargs...)
-
-ocean.output_writers[:snapshots] = JLD2OutputWriter(model, merge(model.tracers, model.velocities);
-                                                    schedule = TimeInterval(10days),
-                                                    filename = "snapshots",
-                                                    output_kwargs...)
-
-# Checkpointer for restarting purposes
-ocean.output_writers[:checkpoint] = Checkpointer(model;
-                                                 schedule = TimeInterval(60days),
-                                                 overwrite_existing = true,
-                                                 prefix = "checkpoint")
-
-#####
-##### Run the simulation!
-#####
-
-# warm up the simulation to ensure that the model adapts 
-# to the interpolated initial conditions without crashing
-ocean.Δt = 10
-ocean.stop_iteration = 1
-wizard = TimeStepWizard(; cfl = 0.1, max_Δt = 90, max_change = 1.1)
-ocean.callbacks[:wizard] = Callback(wizard, IterationInterval(1))
-
-# Finally, the coupled simulation!
-coupled_simulation = Simulation(coupled_model; Δt=1, stop_time = 20days)
-
-run!(coupled_simulation)
-
-#####
-##### Visualization
-#####
-
-using CairoMakie
-
-u, v, w = model.velocities
-T, S, e = model.tracers
-
-using Oceananigans.Models.HydrostaticFreeSurfaceModel: VerticalVorticityField
-
-ζ = VerticalVorticityField(model)
-s = Field(sqrt(u^2 + v^2))
-
-compute!(ζ)
-compute!(s)
-
-ζ = on_architecture(CPU(), ζ)
-s = on_architecture(CPU(), s)
-T = on_architecture(CPU(), T)
-e = on_architecture(CPU(), e)
-
-fig = Figure(size = (1000, 800))
-
-ax = Axis(fig[1, 1], title = "Vertical vorticity [s⁻¹]")
-heatmap!(ax, interior(ζ, :, :, grid.Nz), colorrange = (-4e-5, 4e-5), colormap = :bwr)
-
-ax = Axis(fig[1, 2], title = "Surface speed [ms⁻¹]")
-heatmap!(ax, interior(s, :, :, grid.Nz), colorrange = (0, 0.5), colormap = :deep)
-
-ax = Axis(fig[2, 1], title = "Surface Temperature [Cᵒ]")
-heatmap!(ax, interior(T, :, :, grid.Nz), colorrange = (-1, 30), colormap = :magma)
-
-ax = Axis(fig[2, 1], title = "Turbulent Kinetic Energy [m²s⁻²]")
-heatmap!(ax, interior(e, :, :, grid.Nz), colorrange = (0, 1e-3), colormap = :solar)
diff --git a/examples/single_column_simulation.jl b/examples/single_column_simulation.jl
index e8ef7370..904c3bee 100644
--- a/examples/single_column_simulation.jl
+++ b/examples/single_column_simulation.jl
@@ -267,7 +267,7 @@ lines!(axτ, times, interior(ρτy, 1, 1, 1, :), label="Meridional")
 vlines!(axτ, tn, linewidth=4, color=(:black, 0.5))
 axislegend(axτ)
 
-lines!(axT, times, Tat[1:Nt] .- 273.15,             color=colors[1], linewidth=2, linestyle=:dash, label="Atmosphere temperature")
+lines!(axT, times, Tat[1:Nt] .- 273.15,      color=colors[1], linewidth=2, linestyle=:dash, label="Atmosphere temperature")
 lines!(axT, times, interior(T, 1, 1, Nz, :), color=colors[2], linewidth=4, label="Ocean surface temperature")
 vlines!(axT, tn, linewidth=4, color=(:black, 0.5))
 axislegend(axT)
@@ -279,7 +279,6 @@ lines!(axQ, times, - interior(Qlw, 1, 1, 1, 1:Nt), color=colors[5], label="Longw
 vlines!(axQ, tn, linewidth=4, color=(:black, 0.5))
 axislegend(axQ)
 
-#lines!(axF, times, interior(Jˢt, 1, 1, 1, :), label="Net freshwater flux")
 lines!(axF, times, Pt[1:Nt], label="Prescribed freshwater flux")
 lines!(axF, times, - interior(Ev, 1, 1, 1, 1:Nt), label="Evaporation")
 vlines!(axF, tn, linewidth=4, color=(:black, 0.5))
@@ -324,13 +323,10 @@ Smax = maximum(interior(S))
 Smin = minimum(interior(S))
 xlims!(axSz, Smin - 0.2, Smax + 0.2)
 
-save("single_column_profiles.png", fig)
-nothing # hide
-
-# ![](single_column_profiles.png)
-
-# record(fig, "$(location_name)_single_column_simulation.mp4", 1:8:Nt, framerate=24) do nn
-#     @info "Drawing frame $nn of $Nt..."
-#     n[] = nn
-# end
+record(fig, "single_column_profiles.mp4", 1:8:Nt, framerate=24) do nn
+    @info "Drawing frame $nn of $Nt..."
+    n[] = nn
+end
+nothing #hide
 
+# ![](single_column_profiles.mp4)
\ No newline at end of file
diff --git a/experiments/freely_decaying_regional_simulation.jl b/experiments/freely_decaying_regional_simulation.jl
index 6130064e..fb22ead0 100644
--- a/experiments/freely_decaying_regional_simulation.jl
+++ b/experiments/freely_decaying_regional_simulation.jl
@@ -9,7 +9,7 @@ using ClimaOcean.OceanSeaIceModels: Radiation
 using ClimaOcean.DataWrangling.JRA55: JRA55_prescribed_atmosphere
 using ClimaOcean.DataWrangling.ECCO: ecco_field
 
-# using GLMakie
+# using CairoMakie
 using Printf
 using Dates
 
diff --git a/experiments/one_degree_simulations/inspect_one_degree_near_global_simulation.jl b/experiments/one_degree_simulations/inspect_one_degree_near_global_simulation.jl
index 13654f4e..092d5399 100644
--- a/experiments/one_degree_simulations/inspect_one_degree_near_global_simulation.jl
+++ b/experiments/one_degree_simulations/inspect_one_degree_near_global_simulation.jl
@@ -2,7 +2,7 @@ using Oceananigans
 using Oceananigans.BuoyancyModels: ∂z_b, buoyancy_perturbation
 using SeawaterPolynomials.TEOS10: TEOS10EquationOfState
 using JLD2
-using GLMakie
+using CairoMakie
 using Printf
 
 dir = "../data" #/storage1/greg"
diff --git a/experiments/plot_freely_decaying_simulation.jl b/experiments/plot_freely_decaying_simulation.jl
index 9238a70f..935035c8 100644
--- a/experiments/plot_freely_decaying_simulation.jl
+++ b/experiments/plot_freely_decaying_simulation.jl
@@ -1,4 +1,4 @@
-using GLMakie
+using CairoMakie
 using Oceananigans
 using JLD2
 
diff --git a/prototype_omip_simulation/Manifest.toml b/prototype_omip_simulation/Manifest.toml
index 9ec3daf5..3d4f43c0 100644
--- a/prototype_omip_simulation/Manifest.toml
+++ b/prototype_omip_simulation/Manifest.toml
@@ -180,12 +180,6 @@ git-tree-sha1 = "bcba305388e16aa5c879e896726db9e71b4942c6"
 uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
 version = "0.14.0+1"
 
-[[deps.Cairo_jll]]
-deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "a2f1c8c668c8e3cb4cca4e57a8efdb09067bb3fd"
-uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
-version = "1.18.0+2"
-
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra"]
 git-tree-sha1 = "71acdbf594aab5bbb2cec89b208c41b4c411e49f"
@@ -438,18 +432,6 @@ git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec"
 uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
 version = "0.1.10"
 
-[[deps.FFMPEG]]
-deps = ["FFMPEG_jll"]
-git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
-uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
-version = "0.4.1"
-
-[[deps.FFMPEG_jll]]
-deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
-git-tree-sha1 = "466d45dc38e15794ec7d5d63ec03d776a9aff36e"
-uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
-version = "4.4.4+1"
-
 [[deps.FFTW]]
 deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
 git-tree-sha1 = "4820348781ae578893311153d69049a93d05f39d"
@@ -477,12 +459,6 @@ git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 version = "0.8.5"
 
-[[deps.Fontconfig_jll]]
-deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"]
-git-tree-sha1 = "db16beca600632c95fc8aca29890d83788dd8b23"
-uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
-version = "2.13.96+0"
-
 [[deps.ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
 git-tree-sha1 = "cf0fe81336da9fb90944683b8c41984b08793dad"
@@ -493,18 +469,6 @@ weakdeps = ["StaticArrays"]
     [deps.ForwardDiff.extensions]
     ForwardDiffStaticArraysExt = "StaticArrays"
 
-[[deps.FreeType2_jll]]
-deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "5c1d8ae0efc6c2e7b1fc502cbe25def8f661b7bc"
-uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
-version = "2.13.2+0"
-
-[[deps.FriBidi_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "1ed150b39aebcc805c26b93a8d0122c940f64ce2"
-uuid = "559328eb-81f9-559d-9380-de523a88c83c"
-version = "1.0.14+0"
-
 [[deps.Future]]
 deps = ["Random"]
 uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
@@ -532,12 +496,6 @@ git-tree-sha1 = "518ebd058c9895de468a8c255797b0c53fdb44dd"
 uuid = "61eb1bfa-7361-4325-ad38-22787b887f55"
 version = "0.26.5"
 
-[[deps.Gettext_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
-git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
-uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
-version = "0.21.0+0"
-
 [[deps.Git]]
 deps = ["Git_jll"]
 git-tree-sha1 = "04eff47b1354d702c3a85e8ab23d539bb7d5957e"
@@ -550,12 +508,6 @@ git-tree-sha1 = "d18fb8a1f3609361ebda9bf029b60fd0f120c809"
 uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
 version = "2.44.0+2"
 
-[[deps.Glib_jll]]
-deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "7c82e6a6cd34e9d935e9aa4051b66c6ff3af59ba"
-uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
-version = "2.80.2+0"
-
 [[deps.Glob]]
 git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496"
 uuid = "c27321d9-0574-5035-807b-f59d2c89b15c"
@@ -567,12 +519,6 @@ git-tree-sha1 = "383db7d3f900f4c1f47a8a04115b053c095e48d3"
 uuid = "0951126a-58fd-58f1-b5b3-b08c7c4a876d"
 version = "3.8.4+0"
 
-[[deps.Graphite2_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
-uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
-version = "1.3.14+0"
-
 [[deps.HDF5_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"]
 git-tree-sha1 = "38c8874692d48d5440d5752d6c74b0c6b0b60739"
@@ -585,12 +531,6 @@ git-tree-sha1 = "d1d712be3164d61d1fb98e7ce9bcbc6cc06b45ed"
 uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 version = "1.10.8"
 
-[[deps.HarfBuzz_jll]]
-deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
-git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
-uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
-version = "2.8.1+1"
-
 [[deps.HostCPUFeatures]]
 deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
 git-tree-sha1 = "eb8fed28f4994600e29beef49744639d985a04b2"
@@ -727,12 +667,6 @@ version = "0.9.21"
     [deps.KernelAbstractions.weakdeps]
     EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
 
-[[deps.LAME_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "170b660facf5df5de098d866564877e119141cbd"
-uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
-version = "3.100.2+0"
-
 [[deps.LLVM]]
 deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"]
 git-tree-sha1 = "389aea28d882a40b5e1747069af71bdbd47a1cae"
@@ -769,12 +703,6 @@ weakdeps = ["Serialization"]
     [deps.LRUCache.extensions]
     SerializationExt = ["Serialization"]
 
-[[deps.LZO_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "70c5da094887fd2cae843b8db33920bac4b6f07d"
-uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
-version = "2.10.2+0"
-
 [[deps.LaTeXStrings]]
 git-tree-sha1 = "50901ebc375ed41dbf8058da26f9de442febbbec"
 uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
@@ -822,42 +750,12 @@ version = "1.11.0+1"
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
-[[deps.Libffi_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
-uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
-version = "3.2.2+1"
-
-[[deps.Libgcrypt_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll"]
-git-tree-sha1 = "9fd170c4bbfd8b935fdc5f8b7aa33532c991a673"
-uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
-version = "1.8.11+0"
-
-[[deps.Libgpg_error_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "fbb1f2bef882392312feb1ede3615ddc1e9b99ed"
-uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
-version = "1.49.0+0"
-
 [[deps.Libiconv_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175"
 uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
 version = "1.17.0+0"
 
-[[deps.Libmount_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "0c4f9c4f1a50d8f35048fa0532dabbadf702f81e"
-uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
-version = "2.40.1+0"
-
-[[deps.Libuuid_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "5ee6203157c120d79034c748a2acba45b82b8807"
-uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
-version = "2.40.1+0"
-
 [[deps.LinearAlgebra]]
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -1067,12 +965,6 @@ weakdeps = ["Adapt"]
     [deps.OffsetArrays.extensions]
     OffsetArraysAdaptExt = "Adapt"
 
-[[deps.Ogg_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f"
-uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
-version = "1.3.5+1"
-
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
@@ -1107,24 +999,18 @@ git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
 uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 version = "0.5.5+0"
 
-[[deps.Opus_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
-uuid = "91d4177d-7536-5919-b921-800302f37372"
-version = "1.3.2+0"
-
 [[deps.OrderedCollections]]
 git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 version = "1.6.3"
 
 [[deps.OrthogonalSphericalShellGrids]]
-deps = ["Adapt", "CUDA", "Documenter", "DocumenterTools", "FFMPEG", "JLD2", "JSON3", "KernelAbstractions", "Oceananigans", "OffsetArrays", "Printf"]
-git-tree-sha1 = "1a31c10d9eeec98bfdbb66282b39be643db31fea"
+deps = ["Adapt", "CUDA", "Documenter", "DocumenterTools", "JLD2", "KernelAbstractions", "MPI", "Oceananigans", "OffsetArrays"]
+git-tree-sha1 = "f9614694b607bd2168aef61a819f98a43b9b81c5"
 repo-rev = "main"
-repo-url = "https://github.com/simone-silvestri/OrthogonalSphericalShellGrids.jl.git"
+repo-url = "https://github.com/CliMA/OrthogonalSphericalShellGrids.jl"
 uuid = "c2be9673-fb75-4747-82dc-aa2bb9f4aed0"
-version = "0.1.0"
+version = "0.1.2"
 
 [[deps.P11Kit_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -1177,12 +1063,6 @@ git-tree-sha1 = "bd69f3f0ee248cfb4241800aefb705b5ded592ff"
 uuid = "4a48f351-57a6-4416-9ec4-c37015456aae"
 version = "0.15.1"
 
-[[deps.Pixman_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"]
-git-tree-sha1 = "35621f10a7531bc8fa58f74610b1bfb70a3cfc6b"
-uuid = "30392449-352a-5448-841d-b1acce4e97dc"
-version = "0.43.4+0"
-
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -1637,66 +1517,12 @@ git-tree-sha1 = "52ff2af32e591541550bd753c0da8b9bc92bb9d9"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
 version = "2.12.7+0"
 
-[[deps.XSLT_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
-git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
-uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
-version = "1.1.34+0"
-
 [[deps.XZ_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "ac88fb95ae6447c8dda6a5503f3bafd496ae8632"
 uuid = "ffd25f8a-64ca-5728-b0f7-c24cf3aae800"
 version = "5.4.6+0"
 
-[[deps.Xorg_libX11_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
-git-tree-sha1 = "afead5aba5aa507ad5a3bf01f58f82c8d1403495"
-uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
-version = "1.8.6+0"
-
-[[deps.Xorg_libXau_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "6035850dcc70518ca32f012e46015b9beeda49d8"
-uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
-version = "1.0.11+0"
-
-[[deps.Xorg_libXdmcp_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "34d526d318358a859d7de23da945578e8e8727b7"
-uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
-version = "1.1.4+0"
-
-[[deps.Xorg_libXext_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
-git-tree-sha1 = "d2d1a5c49fae4ba39983f63de6afcbea47194e85"
-uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
-version = "1.3.6+0"
-
-[[deps.Xorg_libXrender_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
-git-tree-sha1 = "47e45cd78224c53109495b3e324df0c37bb61fbe"
-uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
-version = "0.9.11+0"
-
-[[deps.Xorg_libpthread_stubs_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "8fdda4c692503d44d04a0603d9ac0982054635f9"
-uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
-version = "0.1.1+0"
-
-[[deps.Xorg_libxcb_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
-git-tree-sha1 = "b4bfde5d5b652e22b9c790ad00af08b6d042b97d"
-uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
-version = "1.15.0+0"
-
-[[deps.Xorg_xtrans_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "e92a1a012a10506618f10b7047e478403a046c77"
-uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
-version = "1.5.0+0"
-
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
@@ -1714,41 +1540,11 @@ git-tree-sha1 = "46bf7be2917b59b761247be3f317ddf75e50e997"
 uuid = "477f73a3-ac25-53e9-8cc3-50b2fa2566f0"
 version = "1.1.2+0"
 
-[[deps.libaom_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "1827acba325fdcdf1d2647fc8d5301dd9ba43a9d"
-uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b"
-version = "3.9.0+0"
-
-[[deps.libass_jll]]
-deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
-uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
-version = "0.15.1+0"
-
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
 version = "5.8.0+1"
 
-[[deps.libfdk_aac_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
-uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
-version = "2.0.2+0"
-
-[[deps.libpng_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "d7015d2e18a5fd9a4f47de711837e980519781a4"
-uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
-version = "1.6.43+1"
-
-[[deps.libvorbis_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
-git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c"
-uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
-version = "1.3.7+1"
-
 [[deps.libzip_jll]]
 deps = ["Artifacts", "Bzip2_jll", "GnuTLS_jll", "JLLWrappers", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
 git-tree-sha1 = "3282b7d16ae7ac3e57ec2f3fa8fafb564d8f9f7f"
@@ -1770,15 +1566,3 @@ version = "2021.12.0+0"
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
 version = "17.4.0+2"
-
-[[deps.x264_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
-uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
-version = "2021.5.5+0"
-
-[[deps.x265_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
-uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
-version = "3.5.0+0"
diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl
index f3f479ca..3d0993d2 100644
--- a/src/Bathymetry.jl
+++ b/src/Bathymetry.jl
@@ -98,7 +98,13 @@ function regrid_bathymetry(target_grid;
         end
 
         fileurl = joinpath(url, filename)
-        Downloads.download(fileurl, filepath; progress=download_progress, verbose=true)
+
+        try 
+            Downloads.download(fileurl, filepath; progress=download_progress, verbose=true)
+        catch 
+            cmd = `wget --no-check-certificate -O $filepath $fileurl`
+            run(cmd)
+        end
     end
 
     dataset = Dataset(filepath)
diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl
index 86bd2339..9c22df5e 100644
--- a/src/ClimaOcean.jl
+++ b/src/ClimaOcean.jl
@@ -7,6 +7,7 @@ export
     SimilarityTheoryTurbulentFluxes,
     JRA55_prescribed_atmosphere,
     JRA55NetCDFBackend,
+    ECCOMetadata,
     ecco2_field,
     regrid_bathymetry,
     retrieve_bathymetry,
@@ -71,7 +72,7 @@ using .OceanSeaIceModels
 using .OceanSimulations
 using .DataWrangling: JRA55, ECCO
 using ClimaOcean.DataWrangling.JRA55: JRA55_prescribed_atmosphere, JRA55NetCDFBackend
-using ClimaOcean.DataWrangling.ECCO: ecco_field
+using ClimaOcean.DataWrangling.ECCO: ecco_field, ECCOMetadata
 
 end # module
 
diff --git a/src/DataWrangling/JRA55.jl b/src/DataWrangling/JRA55.jl
index a7427ec7..3cfd6dea 100644
--- a/src/DataWrangling/JRA55.jl
+++ b/src/DataWrangling/JRA55.jl
@@ -455,7 +455,7 @@ function JRA55_field_time_series(variable_name;
 
     boundary_conditions = FieldBoundaryConditions(JRA55_native_grid, (Center, Center, Nothing))
     times = jra55_times(native_times)
-    
+
     if backend isa JRA55NetCDFBackend
         fts = FieldTimeSeries{Center, Center, Nothing}(JRA55_native_grid, times;
                                                        backend,
@@ -534,7 +534,7 @@ function JRA55_field_time_series(variable_name;
     else
         copyto!(interior(fts, :, :, 1, :), new_data[:, :, :])
     end
-                     
+
     while n <= all_Nt
         print("        ... processing time index $n of $all_Nt \r")
 
@@ -599,7 +599,7 @@ function JRA55_prescribed_atmosphere(architecture::AA, time_indices=Colon();
                                      backend = nothing,
                                      time_indexing = Cyclical(),
                                      reference_height = 10,  # meters
-                                     include_rivers_and_icebergs = true, # rivers and icebergs are not needed in single column simulations
+                                     include_rivers_and_icebergs = false, # rivers and icebergs are not needed in single column simulations
                                      other_kw...)
 
     if isnothing(backend) # apply a default
@@ -621,7 +621,6 @@ function JRA55_prescribed_atmosphere(architecture::AA, time_indices=Colon();
     va  = JRA55_field_time_series(:northward_velocity;              kw...)
     Ta  = JRA55_field_time_series(:temperature;                     kw...)
     qa  = JRA55_field_time_series(:specific_humidity;               kw...)
-    ra  = JRA55_field_time_series(:relative_humidity;               kw...)
     pa  = JRA55_field_time_series(:sea_level_pressure;              kw...)
     Fra = JRA55_field_time_series(:rain_freshwater_flux;            kw...)
     Fsn = JRA55_field_time_series(:snow_freshwater_flux;            kw...)
@@ -648,10 +647,8 @@ function JRA55_prescribed_atmosphere(architecture::AA, time_indices=Colon();
                   v = va)
 
     tracers = (T = Ta,
-               q = qa,
-               r = ra)
+               q = qa)
 
-                       
     pressure = pa
 
     downwelling_radiation = TwoBandDownwellingRadiation(shortwave=Qs, longwave=Ql)
@@ -672,4 +669,3 @@ function JRA55_prescribed_atmosphere(architecture::AA, time_indices=Colon();
 end
 
 end # module
-
diff --git a/src/DataWrangling/ecco_metadata.jl b/src/DataWrangling/ecco_metadata.jl
index 3855bcac..ceb1d30c 100644
--- a/src/DataWrangling/ecco_metadata.jl
+++ b/src/DataWrangling/ecco_metadata.jl
@@ -1,6 +1,6 @@
 using CFTime
 using Dates
-import Dates: year, month, day
+
 import Oceananigans.Fields: set!
 import Base
 
diff --git a/src/DataWrangling/ecco_restoring.jl b/src/DataWrangling/ecco_restoring.jl
index b7b24b02..fc617961 100644
--- a/src/DataWrangling/ecco_restoring.jl
+++ b/src/DataWrangling/ecco_restoring.jl
@@ -8,7 +8,7 @@ using Base
 
 using NCDatasets
 using JLD2 
-using Dates
+using Dates: Second
 
 using ClimaOcean: stateindex
 
@@ -177,11 +177,11 @@ A struct representing ECCO restoring.
 - `λ⁻¹`: The reciprocal of the restoring timescale.
 """
 struct ECCORestoring{FTS, G, M, V, N} <: Function
-    ecco_fts  :: FTS
-    ecco_grid :: G
-    mask      :: M
-    variable_name   :: V
-    λ⁻¹       :: N
+    ecco_fts      :: FTS
+    ecco_grid     :: G
+    mask          :: M
+    variable_name :: V
+    λ⁻¹           :: N
 end
 
 Adapt.adapt_structure(to, p::ECCORestoring) = 
diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl
index 885898e5..6c479bae 100644
--- a/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl
+++ b/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl
@@ -10,6 +10,8 @@ export Radiation,
 using ..OceanSeaIceModels: default_gravitational_acceleration
 
 import ..OceanSeaIceModels: surface_velocities,
+                            surface_horizontal_velocities,
+                            surface_active_tracers,
                             surface_tracers
 
 import ClimaOcean: stateindex
@@ -27,6 +29,14 @@ function surface_flux(f::Field)
     end
 end
 
+function surface_horizontal_velocities(ocean::Simulation{<:HydrostaticFreeSurfaceModel})
+    grid = ocean.model.grid
+    Nz = size(grid, 3)
+    u = view(ocean.model.velocities.u.data, :, :, Nz)
+    v = view(ocean.model.velocities.v.data, :, :, Nz)
+    return (; u, v)
+end
+
 function surface_velocities(ocean::Simulation{<:HydrostaticFreeSurfaceModel})
     grid = ocean.model.grid
     Nz = size(grid, 3)
@@ -36,6 +46,14 @@ function surface_velocities(ocean::Simulation{<:HydrostaticFreeSurfaceModel})
     return (; u, v, w)
 end
 
+function surface_active_tracers(ocean::Simulation{<:HydrostaticFreeSurfaceModel})
+    grid = ocean.model.grid
+    Nz = size(grid, 3)
+    T = view(ocean.model.tracers.T.data, :, :, Nz)
+    S = view(ocean.model.tracers.S.data, :, :, Nz)
+    return (; T, S)
+end
+
 function surface_tracers(ocean::Simulation{<:HydrostaticFreeSurfaceModel})
     grid = ocean.model.grid
     Nz = size(grid, 3)
diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl
index 682b5f90..36f51782 100644
--- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl
+++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl
@@ -18,8 +18,8 @@ function compute_atmosphere_ocean_fluxes!(coupled_model)
     clock = coupled_model.clock
 
     # Ocean, atmosphere, and sea ice state
-    ocean_velocities  = surface_velocities(ocean)
-    ocean_tracers     = surface_tracers(ocean)
+    ocean_velocities  = surface_horizontal_velocities(ocean)
+    ocean_tracers     = surface_active_tracers(ocean)
 
     # Fluxes, and flux contributors
     centered_velocity_fluxes = (u = coupled_model.fluxes.total.ocean.momentum.uᶜᶜᶜ,
diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl
index 2b6b896b..84d3e4b8 100644
--- a/src/OceanSeaIceModels/OceanSeaIceModels.jl
+++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl
@@ -24,6 +24,8 @@ using KernelAbstractions.Extras.LoopInfo: @unroll
 
 function surface_velocities end
 function surface_tracers end
+function surface_horizontal_velocities end
+function surface_active_tracers end
 function downwelling_radiation end
 function freshwater_flux end
 function reference_density end
@@ -102,4 +104,3 @@ end
 end
 
 end # module
-
diff --git a/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl b/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl
index 5fcf0375..6fab535e 100644
--- a/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl
+++ b/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl
@@ -28,7 +28,7 @@ function limit_fluxes_over_sea_ice!(grid, kernel_parameters, sea_ice::MinimumTem
     launch!(architecture(grid), grid, kernel_parameters, _cap_fluxes_on_sea_ice!,
             centered_velocity_fluxes,
             net_tracer_fluxes,
-            grid, 
+            grid,
             sea_ice.minimum_temperature,
             ocean_temperature)
 
@@ -37,7 +37,7 @@ end
 
 @kernel function _cap_fluxes_on_sea_ice!(centered_velocity_fluxes,
                                          net_tracer_fluxes,
-                                         grid, 
+                                         grid,
                                          minimum_temperature,
                                          ocean_temperature)    
 
@@ -55,8 +55,8 @@ end
         cooling_sea_ice = sea_ice & (Jᵀ[i, j, 1] > 0)
 
         # Don't allow the ocean to cool below the minimum temperature! (make sure it heats up though!)
-	Jᵀ[i, j, 1] = ifelse(cooling_sea_ice, zero(grid), Jᵀ[i, j, 1]) 
-        
+        Jᵀ[i, j, 1] = ifelse(cooling_sea_ice, zero(grid), Jᵀ[i, j, 1]) 
+
         # If we are in a "sea ice" region we remove all fluxes
         Jˢ[i, j, 1] = ifelse(sea_ice, zero(grid), Jˢ[i, j, 1])
         Jᵘ[i, j, 1] = ifelse(sea_ice, zero(grid), Jᵘ[i, j, 1]) 
diff --git a/src/OceanSeaIceModels/ocean_only_model.jl b/src/OceanSeaIceModels/ocean_only_model.jl
index 8ab28eb8..0b32a1a2 100644
--- a/src/OceanSeaIceModels/ocean_only_model.jl
+++ b/src/OceanSeaIceModels/ocean_only_model.jl
@@ -23,14 +23,14 @@ function time_step!(coupled_model::NoSeaIceModel, Δt; callbacks=[], compute_ten
 
     tick!(coupled_model.clock, ocean.Δt) # An Ocean-only model advances with the ocean time-step!
     update_state!(coupled_model, callbacks; compute_tendencies)
-    
+
     return nothing
 end
 
 function update_state!(coupled_model::NoSeaIceModel, callbacks=[]; compute_tendencies=false)
     time = Time(coupled_model.clock.time)
     update_model_field_time_series!(coupled_model.atmosphere, time)
-    
+
     ocean_model = coupled_model.ocean.model
 
     # Do we really have to do this?
@@ -42,4 +42,3 @@ function update_state!(coupled_model::NoSeaIceModel, callbacks=[]; compute_tende
     compute_atmosphere_ocean_fluxes!(coupled_model) 
     return nothing
 end
-