diff --git a/examples/freely_decaying_mediterraneum.jl b/examples/freely_decaying_mediterraneum.jl
index ac068612..27fcf15f 100644
--- a/examples/freely_decaying_mediterraneum.jl
+++ b/examples/freely_decaying_mediterraneum.jl
@@ -8,6 +8,7 @@ using ClimaOcean.VerticalGrids: stretched_vertical_faces, PowerLawStretching
 using ClimaOcean.InitialConditions: three_dimensional_regrid!, adjust_tracers!
 using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity
 using Oceananigans.Coriolis: ActiveCellEnstrophyConserving
+using Oceananigans.Units
 
 #####
 ##### Regional Mediterranean grid 
@@ -42,7 +43,7 @@ Nx = 20 * 42 # 1 / 20th of a degree
 Ny = 20 * 15 # 1 / 20th of a degree
 Nz = length(z) - 1
 
-grid = LatitudeLongitudeGrid(GPU();
+grid = LatitudeLongitudeGrid(CPU();
                              size = (Nx, Ny, Nz),
                              latitude = (30, 45),
                              longitude = (0, 42),
diff --git a/examples/generate_bathymetry.jl b/examples/generate_bathymetry.jl
index 64ebc579..a989069f 100644
--- a/examples/generate_bathymetry.jl
+++ b/examples/generate_bathymetry.jl
@@ -35,8 +35,12 @@ display(fig)
 ##### Regional Mediterranean grid 
 #####
 
+# 1/25th degree resolution
+Nλ = 25 * 55
+Nφ = 25 * 25
+
 grid = LatitudeLongitudeGrid(CPU();
-                             size = (25 * 50, 55 * 50, 1),
+                             size = (Nλ, Nφ, 1),
                              latitude = (25, 50),
                              longitude = (-10, 45),
                              z = (0, 1),
diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl
index 37ae85e5..c4232195 100644
--- a/src/Bathymetry.jl
+++ b/src/Bathymetry.jl
@@ -4,8 +4,13 @@ using ..DataWrangling: download_progress
 
 using Oceananigans
 using Oceananigans.Architectures: architecture
-using Oceananigans.Grids: halo_size, λnodes
-using Oceananigans.Utils: pretty_filesize
+using Oceananigans.Grids: halo_size, λnodes, φnodes
+using Oceananigans.Grids: x_domain, y_domain
+using Oceananigans.Grids: topology
+using Oceananigans.Utils: pretty_filesize, launch!
+using Oceananigans.Fields: interpolate!
+using Oceananigans.BoundaryConditions
+using KernelAbstractions: @kernel, @index
 
 using NCDatasets
 using Downloads
@@ -22,14 +27,38 @@ using Printf
 Regrid bathymetry associated with the NetCDF file at `path = joinpath(dir, filename)` to `target_grid`.
 If `path` does not exist, then a download is attempted from `joinpath(url, filename)`.
 
-TODO: describe keyword arguments.
+Arguments:
+==========
+
+- target_grid: grid to interpolate onto
+
+Keyword Arguments:
+==================
+
+- height_above_water: limits the maximum height of above-water topography (where h > 0). If
+                      `nothing` the original topography is retained
+
+- minimum_depth: minimum depth for the shallow regions. `h > minimum_depth` will be considered land
+
+- dir: directory of the bathymetry-containing file
+
+- filename: file containing bathymetric data. Must be netcdf with fields:
+            (1) `lat` vector of latitude nodes
+            (2) `lon` vector of longitude nodes
+            (3) `z` matrix of depth values
+
+- interpolation_passes: regridding/interpolation passes. The bathymetry is interpolated in
+                        `interpolation_passes - 1` intermediate steps. With more steps the 
+                        final bathymetry will be smoother.
+
 """
 function regrid_bathymetry(target_grid;
                            height_above_water = nothing,
                            minimum_depth = 0,
                            dir = joinpath(@__DIR__, "..", "data"),
                            url = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf", 
-                           filename = "ETOPO_2022_v1_60s_N90W180_surface.nc")
+                           filename = "ETOPO_2022_v1_60s_N90W180_surface.nc",
+                           interpolation_passes = 10)
 
     filepath = joinpath(dir, filename)
 
@@ -106,8 +135,8 @@ function regrid_bathymetry(target_grid;
     end
 
     if minimum_depth > 0
-        ocean = h_data .<= 0
-        h_data[ocean] .= min.(-minimum_depth, h_data[ocean])
+        shallow_ocean = h_data .> minimum_depth
+        h_data[shallow_ocean] .= height_above_water
     end
 
     # Build the "native" grid of the bathymetry and make a bathymetry field.
@@ -125,28 +154,57 @@ function regrid_bathymetry(target_grid;
     native_h = Field{Center, Center, Nothing}(native_grid)
     set!(native_h, h_data)
 
-    Nxi = Nxt
-    Nyi = Nyn
-    Nzi = Nzn
+    # 
+    target_h = interpolate_bathymetry_in_passes(native_h, target_grid; passes = interpolation_passes)
 
-    if parent(parent(λt)) isa StepRangeLen # longitude is equispaced
-        longitude = (λ₁, λ₂)
-    else
-        longitude = λnodes(target_grid, Face(), Center(), Center())
+    return target_h
+end
+
+# Here we can either use `regrid!` (three dimensional version) or `interpolate`
+function interpolate_bathymetry_in_passes(native_h, target_grid; passes = 10)
+
+    Nλt, Nφt = Nt = size(target_grid)
+    Nλn, Nφn = Nn = size(native_h)
+    
+    if any(Nt .> Nn) # We are refining the grid (at least in one direction), more passes will not help!
+        new_h = Field{Center, Center, Nothing}(target_grid)
+        interpolate!(new_h, native_h)
+        return new_h
     end
+ 
+    latitude  = y_domain(target_grid)
+    longitude = x_domain(target_grid)
+
+    ΔNλ = floor((Nλn - Nλt) / passes)
+    ΔNφ = floor((Nφn - Nφt) / passes)
+
+    Nλ = [Nλn - ΔNλ * pass for pass in 1:passes-1]
+    Nφ = [Nφn - ΔNφ * pass for pass in 1:passes-1]
 
-    intermediate_grid = LatitudeLongitudeGrid(arch;
-                                              size = (Nxi, Nyi, Nzi),
-                                              latitude = (φ₁, φ₂),
-                                              longitude,
-                                              z = (0, 1),
-                                              halo = halo_size(target_grid))
+    Nλ = Int[Nλ..., Nλt]
+    Nφ = Int[Nφ..., Nφt]
 
-    intermediate_h =  Field{Center, Center, Nothing}(intermediate_grid)
-    regrid!(intermediate_h, native_h)
+    old_h     = native_h
+    TX, TY, _ = topology(target_grid)
+
+    for pass = 1:passes - 1
+        new_size = (Nλ[pass], Nφ[pass], 1)
+
+        @info "pass number $pass with size $new_size"
+        new_grid = LatitudeLongitudeGrid(size = new_size, 
+                                     latitude = (latitude[1],  latitude[2]), 
+                                    longitude = (longitude[1], longitude[2]), 
+                                            z = (0, 1),
+                                     topology = (TX, TY, Bounded))
+
+        new_h = Field{Center, Center, Nothing}(new_grid)
+
+        interpolate!(new_h, old_h)
+        old_h = new_h
+    end
 
     target_h = Field{Center, Center, Nothing}(target_grid)
-    regrid!(target_h, intermediate_h)
+    interpolate!(target_h, old_h)
 
     return target_h
 end
diff --git a/src/DataWrangling/ECCO2.jl b/src/DataWrangling/ECCO2.jl
index 186f72d2..71ee81eb 100644
--- a/src/DataWrangling/ECCO2.jl
+++ b/src/DataWrangling/ECCO2.jl
@@ -2,6 +2,7 @@ module ECCO2
 
 using Oceananigans
 using Oceananigans.BoundaryConditions
+using KernelAbstractions: @kernel, @index
 using NCDatasets
 
 temperature_filename = "THETA.1440x720x50.19920102.nc"
diff --git a/src/InitialConditions.jl b/src/InitialConditions.jl
index 53407763..774d38ed 100644
--- a/src/InitialConditions.jl
+++ b/src/InitialConditions.jl
@@ -50,6 +50,11 @@ end
 
 propagate_horizontally!(field, ::Nothing; kw...) = nothing
 
+""" 
+    propagate_horizontally!(field, mask; max_iter = Inf)
+
+propagate horizontally 
+"""
 function propagate_horizontally!(field, mask; max_iter = Inf) 
     iter  = 0
     grid  = field.grid