diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index 4d243b46..ef7b558a 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -1,74 +1,108 @@
 using ClimaOcean
-using ClimaOcean.ECCO: ECCO4Monthly
+using ClimaOcean.ECCO: ECCO4Monthly, NearestNeighborInpainting
 using OrthogonalSphericalShellGrids
 using Oceananigans
 using Oceananigans.Units
 using CFTime
 using Dates
 using Printf
+using CUDA: @allowscalar, device!
 
-arch = CPU()
-z = exponential_z_faces(Nz=30, depth=6000)
-Nx = 120
-Ny = 60
-Nz = length(z) - 1
+using Oceananigans.Grids: znode
 
-grid = TripolarGrid(arch; z, size = (Nx, Ny, Nz), north_poles_latitude=55, first_pole_longitude=70)
+device!(3)
+arch = GPU()
 
-@show grid
+#####
+##### Grid and Bathymetry
+#####
 
-bottom_height = regrid_bathymetry(grid;
+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 = 5,
-                                  major_basins = 3)
+                                  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))
+
+#####
+##### Closures
+#####
 
-# Closure
 gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000)
-catke = Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivity()
-closure = (gm, catke)
-
-# Polar restoring
-@inline function restoring_mask(λ, φ, z, t=0)
-    ϵN = (φ - 75) / 5
-    ϵN = clamp(ϵN, zero(ϵN), one(ϵN))
-    ϵS = - (φ + 75) / 5
-    ϵS = clamp(ϵS, zero(ϵS), one(ϵS))
-    return ϵN + ϵS
-end
+catke = ClimaOcean.OceanSimulations.default_ocean_closure()
+viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=10000)
+
+closure = (gm, catke, viscous_closure)
 
-restoring_rate = 1 / 1days
+#####
+##### Restoring
+#####
 
-restoring_mask_field = CenterField(grid)
-set!(restoring_mask_field, restoring_mask)
+restoring_rate  = 1 / 2days
+z_below_surface = @allowscalar znode(1, 1, grid.Nz, grid, Center(), Center(), Face())
 
-@inline sponge_layer(λ, φ, z, t, c, ω) = - restoring_mask(λ, φ, z, t) * ω * c
-Fu = Forcing(sponge_layer, field_dependencies=:u, parameters=restoring_rate)
-Fv = Forcing(sponge_layer, field_dependencies=:v, parameters=restoring_rate)
+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, ECCO4Monthly())
-salinity = ECCOMetadata(:salinity, dates, ECCO4Monthly())
+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="./")
 
-FT = ECCORestoring(arch, temperature; grid, mask=restoring_mask_field, rate=restoring_rate)
-FS = ECCORestoring(arch, salinity;    grid, mask=restoring_mask_field, rate=restoring_rate)
-forcing = (T=FT, S=FS, u=Fu, v=Fv)
+# inpainting = NearestNeighborInpainting(30) should be enough to fill the gaps near bathymetry
+FT = ECCORestoring(arch, temperature; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(30))
+FS = ECCORestoring(arch, salinity;    grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(30))
+forcing = (T=FT, S=FS)
+
+#####
+##### Ocean simulation
+##### 
 
 momentum_advection = VectorInvariant()
-tracer_advection = Centered(order=2)
+tracer_advection   = Centered(order=2)
+
+# Should we add a side drag since this is at a coarser resolution?
 ocean = ocean_simulation(grid; momentum_advection, tracer_advection,
                          closure, forcing,
                          tracers = (:T, :S, :e))
 
-set!(ocean.model,
-     T = ECCOMetadata(:temperature; dates=first(dates)),
-     S = ECCOMetadata(:salinity; dates=first(dates)))
+set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)),
+                  S=ECCOMetadata(:salinity;    dates=first(dates)))
+
+#####
+##### Atmospheric forcing
+#####
+
+radiation  = Radiation(arch)
+atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20))
+
+#####
+##### Coupled simulation
+#####
 
-radiation = Radiation(arch)
-atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(41))
 sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce()
 coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
 
-simulation = Simulation(coupled_model; Δt=1, stop_iteration=10)
+simulation = Simulation(coupled_model; Δt=15minutes, stop_time=2*365days)
+
+#####
+##### Run it!
+##### 
 
 wall_time = Ref(time_ns())
 
@@ -76,9 +110,12 @@ function progress(sim)
     ocean = sim.model.ocean
     u, v, w = ocean.model.velocities
     T = ocean.model.tracers.T
-    Tmax = maximum(T)
-    Tmin = minimum(T)
-    umax = maximum(abs, u), maximum(abs, v), maximum(abs, w)
+    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[])
 
     @info @sprintf("Time: %s, n: %d, Δt: %s, max|u|: (%.2e, %.2e, %.2e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, wall time: %s \n",
@@ -90,6 +127,6 @@ function progress(sim)
      return nothing
 end
 
-add_callback!(simulation, progress, IterationInterval(1))
+add_callback!(simulation, progress, IterationInterval(10))
 
 run!(simulation)
diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl
index 35100b0e..fc22cfad 100644
--- a/src/DataWrangling/ECCO/ECCO.jl
+++ b/src/DataWrangling/ECCO/ECCO.jl
@@ -18,12 +18,14 @@ using Oceananigans.Utils
 
 using KernelAbstractions: @kernel, @index
 using NCDatasets
+using JLD2
 using Downloads: download
 using Dates
 using Adapt
 using Scratch
 
 download_ECCO_cache::String = ""
+
 function __init__()
     global download_ECCO_cache = @get_scratch!("ECCO")
 end
@@ -123,43 +125,68 @@ end
 """
     ECCO_field(metadata::ECCOMetadata;
                architecture = CPU(),
-               horizontal_halo = (3, 3))
-
-Retrieve the ecco field corresponding to `metadata`. 
-The data is loaded from `filename` on `architecture` with `horizontal_halo`
-in the x and y direction. The halo in the z-direction is one.
+               inpainting = nothing,
+               mask = nothing,
+               horizontal_halo = (7, 7),
+               cache_inpainted_data = false)
+
+Return a `Field` on `architecture` described by `ECCOMetadata` with
+`horizontal_halo` size.
+If not `nothing`, the `inpainting` method is used to fill the cells
+within the specified `mask`. `mask` is set to `ECCO_mask` for non-nothing
+`inpainting`.
 """
 function ECCO_field(metadata::ECCOMetadata;
                     architecture = CPU(),
-                    horizontal_halo = (7, 7))
+                    inpainting = nothing,
+                    mask = nothing,
+                    horizontal_halo = (7, 7),
+                    cache_inpainted_data = false)
+
+    # Respect user-supplied mask, but otherwise build default ECCO mask.
+    if !isnothing(inpainting) && isnothing(mask)
+        mask = ECCO_mask(metadata, architecture)
+    end
+
+    field = empty_ECCO_field(metadata; architecture, horizontal_halo)
+    inpainted_path = inpainted_metadata_path(metadata)
+
+    if !isnothing(inpainting) && isfile(inpainted_path)
+        file = jldopen(inpainted_path, "r")
+        maxiter = file["inpainting_maxiter"]
+
+        # read data if generated with the same inpainting
+        if maxiter == inpainting.maxiter
+            data = file["data"]
+            close(file)
+            copyto!(parent(field), data)
+            return field
+        end
+    end
 
     download_dataset(metadata)
     path = metadata_path(metadata)
     ds = Dataset(path)
-
     shortname = short_name(metadata)
 
     if variable_is_three_dimensional(metadata)
         data = ds[shortname][:, :, :, 1]
-
-        # The surface layer in three-dimensional ECCO fields is at `k = 1`
         data = reverse(data, dims=3)
     else
         data = ds[shortname][:, :, 1]
     end        
 
     close(ds)
-
-    field = empty_ECCO_field(metadata; architecture, horizontal_halo)
     
+    # Convert data from Union(FT, missing} to FT
     FT = eltype(field)
     data[ismissing.(data)] .= 1e10 # Artificially large number!
-    data = if location(field)[2] == Face
+    data = if location(field)[2] == Face # ?
         new_data = zeros(FT, size(field))
         new_data[:, 1:end-1, :] .= data
         new_data    
     else
-        convert.(FT, data)
+        data = Array{FT}(data)
     end
     
     # ECCO4 data is on a -180, 180 longitude grid as opposed to ECCO2 data that
@@ -173,67 +200,53 @@ function ECCO_field(metadata::ECCOMetadata;
     set!(field, data)
     fill_halo_regions!(field)
 
+    if !isnothing(inpainting)
+        # Make sure all values are extended properly
+        name = string(metadata.name)
+        date = string(metadata.dates)
+        version = summary(metadata.version)
+        @info string("Inpainting ", version, " ", name, " data from ", date, "...")
+        start_time = time_ns()
+        
+        inpaint_mask!(field, mask; inpainting)
+        fill_halo_regions!(field)
+
+        elapsed = 1e-9 * (time_ns() - start_time)
+        @info string(" ... (", prettytime(elapsed), ")")
+    
+        # We cache the inpainted data to avoid recomputing it
+        @root if cache_inpainted_data
+            file = jldopen(inpainted_path, "w+")
+            file["data"] = on_architecture(CPU(), parent(field))
+            file["inpainting_maxiter"] = inpainting.maxiter
+            close(file)
+        end
+    end
+
     return field
 end
 
 # Fallback
 ECCO_field(var_name::Symbol; kw...) = ECCO_field(ECCOMetadata(var_name); kw...)
 
-"""
-    inpainted_ECCO_field(metadata::ECCOMetadata;
-                         architecture = CPU(),
-                         mask = ECCO_mask(metadata, architecture),
-                         inpainting = NearestNeighborInpainting(Inf),
-                         kw...)
-    
-Retrieve the ECCO field corresponding to `metadata` inpainted to fill all the missing
-values in the original dataset.
-
-Arguments
-=========
-
-- `metadata`: the metadata corresponding to the dataset.
-
-Keyword Arguments
-=================
-
-- `architecture`: either `CPU()` or `GPU()`.
-- `mask`: the mask used to inpaint the field, see [`inpaint_mask!`](@ref).
-- `inpainting`: the inpainting algorithm, see [`inpaint_mask!`](@ref). Default: `NearestNeighborInpainting(Inf)`.
-"""
-function inpainted_ECCO_field(metadata::ECCOMetadata; 
-                              architecture = CPU(),
-                              mask = ECCO_mask(metadata, architecture),
-                              inpainting = NearestNeighborInpainting(Inf),
-                              kw...)
-
-    # Make sure all values are extended properly
-    name = string(metadata.name)
-    date = string(metadata.dates)
-    version = summary(metadata.version)
-    @info string("Inpainting ", version, " ", name, " data from ", date, "...")
-    start_time = time_ns()
-    
-    f = ECCO_field(metadata; architecture, kw...)
-    inpaint_mask!(f, mask; inpainting)
-    fill_halo_regions!(f)
-    elapsed = 1e-9 * (time_ns() - start_time)
-    @info string(" ... (", prettytime(elapsed), ")")
-    return f
+function inpainted_metadata_filename(metadata::ECCOMetadata)
+    original_filename = metadata_filename(metadata)
+    without_extension = original_filename[1:end-3]
+    return without_extension * "_inpainted.jld2"
 end
 
-inpainted_ECCO_field(variable_name::Symbol; kw...) = inpainted_ECCO_field(ECCOMetadata(variable_name); kw...)
-    
-function set!(field::Field, ecco_metadata::ECCOMetadata; kw...)
+inpainted_metadata_path(metadata::ECCOMetadata) = joinpath(metadata.dir, inpainted_metadata_filename(metadata))
+
+function set!(field::Field, ECCO_metadata::ECCOMetadata; kw...)
 
     # Fields initialized from ECCO
     grid = field.grid
     arch = child_architecture(grid)
-    mask = ECCO_mask(ecco_metadata, arch)
+    mask = ECCO_mask(ECCO_metadata, arch)
 
-    f = inpainted_ECCO_field(ecco_metadata; mask,
-                             architecture = arch,
-                             kw...)
+    f = ECCO_field(ECCO_metadata; mask,
+                   architecture = arch,
+                   kw...)
 
     interpolate!(field, f)
 
@@ -241,3 +254,4 @@ function set!(field::Field, ecco_metadata::ECCOMetadata; kw...)
 end
 
 end # Module 
+
diff --git a/src/DataWrangling/ECCO/ECCO_metadata.jl b/src/DataWrangling/ECCO/ECCO_metadata.jl
index 27125914..a0f07101 100644
--- a/src/DataWrangling/ECCO/ECCO_metadata.jl
+++ b/src/DataWrangling/ECCO/ECCO_metadata.jl
@@ -23,7 +23,7 @@ Metadata information for an ECCO dataset:
 - `dir`: The directory where the dataset is stored.
 """
 struct ECCOMetadata{D, V}
-    name  :: Symbol
+    name :: Symbol
     dates :: D
     version :: V
     dir :: String
diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index fc08b53f..a19916bd 100644
--- a/src/DataWrangling/ECCO/ECCO_restoring.jl
+++ b/src/DataWrangling/ECCO/ECCO_restoring.jl
@@ -19,20 +19,20 @@ import Oceananigans.OutputReaders: new_backend, update_field_time_series!
 @inline instantiate(T::DataType) = T()
 @inline instantiate(T) = T
 
-struct ECCONetCDFBackend{N, I, M} <: AbstractInMemoryBackend{Int}
+struct ECCONetCDFBackend{N, C, I, M} <: AbstractInMemoryBackend{Int}
     start :: Int
     length :: Int
     inpainting :: I
     metadata :: M
 
-    function ECCONetCDFBackend{N}(start::Int, length::Int, inpainting, metadata) where N
+    function ECCONetCDFBackend{N, C}(start::Int, length::Int, inpainting, metadata) where {N, C}
         M = typeof(metadata)
         I = typeof(inpainting)
-        return new{N, I, M}(start, length, inpainting, metadata)
+        return new{N, C, I, M}(start, length, inpainting, metadata)
     end
 end
 
-Adapt.adapt_structure(to, b::ECCONetCDFBackend{N}) where N = ECCONetCDFBackend{N}(b.start, b.length, nothing, nothing)
+Adapt.adapt_structure(to, b::ECCONetCDFBackend{N, C}) where {N, C} = ECCONetCDFBackend{N, C}(b.start, b.length, nothing, nothing)
 
 """
     ECCONetCDFBackend(length; on_native_grid = false, inpainting = NearestNeighborInpainting(Inf))
@@ -42,9 +42,10 @@ Each time instance is stored in an individual file.
 """
 function ECCONetCDFBackend(length, metadata;
                            on_native_grid = false, 
+                           cache_inpainted_data = false,
                            inpainting = NearestNeighborInpainting(Inf))
 
-    return ECCONetCDFBackend{on_native_grid}(1, length, inpainting, metadata)
+    return ECCONetCDFBackend{on_native_grid, cache_inpainted_data}(1, length, inpainting, metadata)
 end
 
 Base.length(backend::ECCONetCDFBackend)  = backend.length
@@ -52,21 +53,23 @@ Base.summary(backend::ECCONetCDFBackend) = string("ECCONetCDFBackend(", backend.
 
 const ECCONetCDFFTS{N} = FlavorOfFTS{<:Any, <:Any, <:Any, <:Any, <:ECCONetCDFBackend{N}} where N
 
-new_backend(b::ECCONetCDFBackend{native}, start, length) where native =
+new_backend(b::ECCONetCDFBackend{native, cache_data}, start, length) where {native, cache_data} =
     ECCONetCDFBackend{native}(start, length, b.inpainting, b.metadata)
 
 on_native_grid(::ECCONetCDFBackend{native}) where native = native
+cache_inpainted_data(::ECCONetCDFBackend{native, cache_data}) where {native, cache_data} = cache_data
 
 function set!(fts::ECCONetCDFFTS) 
     backend = fts.backend
     start   = backend.start
     inpainting = backend.inpainting
+    cache_data = cache_inpainted_data(backend)
     len = backend.length
 
     for t in start:start+len-1
         # Set each element of the time-series to the associated file
         metadatum = @inbounds backend.metadata[t] 
-        set!(fts[t], metadatum; inpainting)
+        set!(fts[t], metadatum; inpainting, cache_inpainted_data=cache_data)
     end
 
     fill_halo_regions!(fts)
@@ -130,12 +133,17 @@ Keyword Arguments
 - `inpainting`: The inpainting algorithm to use for ECCO interpolation.
                 The only option is `NearestNeighborInpainting(maxiter)`, 
                 where an average of the valid surrounding values is used `maxiter` times.
+
+- `cache_inpainted_data`: If `true`, the data is cached to disk after inpainting for later retrieving. 
+                          Default: `true`.
+
 """
 function ECCO_field_time_series(metadata::ECCOMetadata;	
                                 architecture = CPU(),	
                                 time_indices_in_memory = 2,	
                                 time_indexing = Cyclical(),
                                 inpainting = NearestNeighborInpainting(prod(size(metadata))),
+                                cache_inpainted_data = true,
                                 grid = nothing)	
 
     # Make sure all the required individual files are downloaded
@@ -151,7 +159,7 @@ function ECCO_field_time_series(metadata::ECCOMetadata;
     loc = LX, LY, LZ = location(metadata)
     boundary_conditions = FieldBoundaryConditions(grid, loc)
 
-    backend = ECCONetCDFBackend(time_indices_in_memory, metadata; on_native_grid, inpainting)
+    backend = ECCONetCDFBackend(time_indices_in_memory, metadata; on_native_grid, inpainting, cache_inpainted_data)
 
     fts = FieldTimeSeries{LX, LY, LZ}(grid, times; backend, time_indexing, boundary_conditions)
     set!(fts)	
@@ -242,7 +250,8 @@ end
                   mask = 1,
                   rate = 1,
                   grid = nothing,
-                  inpainting = NearestNeighborInpainting(prod(size(metadata))))
+                  inpainting = NearestNeighborInpainting(prod(size(metadata))),
+                  cache_inpainted_data = true)
 
 Create a forcing term that restores to values stored in an ECCO field time series.
 The restoring is applied as a forcing on the right hand side of the evolution equations calculated as
@@ -284,14 +293,15 @@ Keyword Arguments
 
 - `time_indices_in_memory:` how many time instances are loaded in memory; the remaining are loaded lazily.
 
-- `grid`: if not a `nothing`, the ECCO data is directly interpolated on the `grid`.
-
 - `inpainting`: inpainting algorithm, see [`inpaint_mask!`](@ref). Default: `NearestNeighborInpainting(Inf)`.
 
 - `grid`: If `isnothing(grid)`, ECCO data is interpolated on-the-fly to the simulation grid.
           If `!isnothing(grid)`, ECCO data is pre-interpolated to `grid`.
           Default: nothing.
 
+- `cache_inpainted_data`: If `true`, the data is cached to disk after inpainting for later retrieving. 
+                          Default: `true`.
+                          
 It is possible to also pass an `ECCOMetadata` type as the first argument without the need for the 
 `variable_name` argument and the `version` and `dates` keyword arguments.
 """
@@ -305,22 +315,38 @@ function ECCORestoring(arch::AbstractArchitecture,
     return ECCORestoring(arch, metadata; kw...)
 end
 
-function ECCORestoring(architecture::AbstractArchitecture,
+function ECCORestoring(arch::AbstractArchitecture,
                        metadata::ECCOMetadata;
                        rate,
                        mask = 1,
                        grid = nothing,
                        time_indices_in_memory = 2, # Not more than this if we want to use GPU!
                        time_indexing = Cyclical(),
-                       inpainting = NearestNeighborInpainting(Inf))
+                       inpainting = NearestNeighborInpainting(Inf),
+                       cache_inpainted_data = true)
 
-    fts = ECCO_field_time_series(metadata; grid, architecture, time_indices_in_memory, time_indexing, inpainting)
+    # Validate architecture
+    if !isnothing(grid) && architecture(grid) != arch
+        throw(ArgumentError("The architecture of ECCORestoring must match the architecture of the grid."))
+    end
+
+    fts = ECCO_field_time_series(metadata; 
+                                 grid, 
+                                 architecture = arch, 
+                                 time_indices_in_memory, 
+                                 time_indexing, 
+                                 inpainting,
+                                 cache_inpainted_data)
 
     # Grab the correct Oceananigans field to restore
     variable_name = metadata.name
     field_name = oceananigans_fieldname[variable_name]
 
-    return ECCORestoring(fts, fts.grid, mask, field_name, rate)
+    # If we pass the grid we do not need to interpolate
+    # so we can save parameter space by setting the native grid to nothing
+    native_grid = isnothing(grid) ? fts.grid : nothing 
+
+    return ECCORestoring(fts, native_grid, mask, field_name, rate)
 end
 
 # Make sure we can call ECCORestoring with architecture as the first positional argument
diff --git a/src/DataWrangling/JRA55.jl b/src/DataWrangling/JRA55.jl
index d4493bca..ce8daa0f 100644
--- a/src/DataWrangling/JRA55.jl
+++ b/src/DataWrangling/JRA55.jl
@@ -31,6 +31,7 @@ import Oceananigans.OutputReaders: new_backend, update_field_time_series!
 using Downloads: download
 
 download_jra55_cache::String = ""
+
 function __init__()
     global download_jra55_cache = @get_scratch!("JRA55")
 end
diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl
index dfee7ad7..9d583e6e 100644
--- a/src/OceanSimulations/OceanSimulations.jl
+++ b/src/OceanSimulations/OceanSimulations.jl
@@ -44,7 +44,8 @@ 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?
-default_free_surface(grid::TripolarGrid) = SplitExplicitFreeSurface(grid; substeps = 70)
+const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}}
+default_free_surface(grid::TripolarOfSomeKind) = SplitExplicitFreeSurface(grid; substeps=70)
 
 function default_ocean_closure()
     mixing_length = CATKEMixingLength(Cᵇ=0.01)
diff --git a/test/test_ecco.jl b/test/test_ecco.jl
index 061657d5..9a3bb9c3 100644
--- a/test/test_ecco.jl
+++ b/test/test_ecco.jl
@@ -21,37 +21,42 @@ inpainting = NearestNeighborInpainting(2)
 @testset "ECCO fields utilities" begin
     for arch in test_architectures
         A = typeof(arch)
-        @info "Testing ECCO_field on $A..."
-
-        temperature = ECCOMetadata(:temperature, dates, ECCO4Monthly())
-        t_restoring = ECCORestoring(temperature; rate = 1 / 1000.0, inpainting)
-
-        ECCO_fts = t_restoring.field_time_series
-
-        for metadatum in temperature
-            @test isfile(metadata_path(metadatum))
+        for name in (:temperature, :salinity)
+            @info "Testing ECCO_field on $A..."
+            metadata = ECCOMetadata(name, dates, ECCO4Monthly())
+            restoring = ECCORestoring(metadata ; rate = 1 / 1000.0, inpainting)
+
+            for datum in metadata 
+                @test isfile(metadata_path(datum))
+            end
+
+            fts = restoring.field_time_series
+            @test fts isa FieldTimeSeries
+            @test fts.grid isa LatitudeLongitudeGrid
+            @test topology(fts.grid) == (Periodic, Bounded, Bounded)
+
+            Nx, Ny, Nz = size(interior(fts))
+            Nt = length(fts.times)
+
+            @test Nx == size(metadata)[1]
+            @test Ny == size(metadata)[2]
+            @test Nz == size(metadata)[3]
+            @test Nt == size(metadata)[4]
+
+            @test fts.times[1] == ECCO_times(metadata)[1]
+            @test fts.times[end] == ECCO_times(metadata)[end]
+
+            datum = first(metadata)
+            ψ = ECCO_field(datum, architecture=arch, inpainting=NearestNeighborInpainting(2))
+            datapath = ClimaOcean.DataWrangling.ECCO.inpainted_metadata_path(datum)
+            @test isfile(datapath)
         end
-
-        @test ECCO_fts isa FieldTimeSeries
-        @test ECCO_fts.grid isa LatitudeLongitudeGrid
-        @test topology(ECCO_fts.grid) == (Periodic, Bounded, Bounded)
-
-        Nx, Ny, Nz = size(interior(ECCO_fts))
-        Nt = length(ECCO_fts.times)
-
-        @test Nx == size(temperature)[1]
-        @test Ny == size(temperature)[2]
-        @test Nz == size(temperature)[3]
-        @test Nt == size(temperature)[4]
-
-        @test ECCO_fts.times[1]   == ECCO_times(temperature)[1]
-        @test ECCO_fts.times[end] == ECCO_times(temperature)[end]
     end
 end
 
 @testset "Inpainting algorithm" begin
     for arch in test_architectures
-        temperature = ECCOMetadata(:temperature, dates[1], ECCO4Monthly())
+        T_metadata = ECCOMetadata(:temperature, dates[1], ECCO4Monthly())
 
         grid = LatitudeLongitudeGrid(arch,
                                      size = (100, 100, 10),
@@ -63,8 +68,8 @@ end
         fully_inpainted_field = CenterField(grid)
         partially_inpainted_field = CenterField(grid)
 
-        set!(fully_inpainted_field,     temperature; inpainting = NearestNeighborInpainting(Inf))
-        set!(partially_inpainted_field, temperature; inpainting = NearestNeighborInpainting(1))
+        set!(fully_inpainted_field,     T_metadata; inpainting = NearestNeighborInpainting(Inf))
+        set!(partially_inpainted_field, T_metadata; inpainting = NearestNeighborInpainting(1))
 
         fully_inpainted_interior = on_architecture(CPU(), interior(fully_inpainted_field))
         partially_inpainted_interior = on_architecture(CPU(), interior(partially_inpainted_field))