From 6adbbbfddabfac1a0ebef516f076a360c4087d50 Mon Sep 17 00:00:00 2001
From: Gregory Wagner <wagner.greg@gmail.com>
Date: Sat, 9 Nov 2024 23:26:49 -0700
Subject: [PATCH 01/30] Update compat

---
 Project.toml                                          |  4 ++--
 .../one_degree_simulation/one_degree_simulation.jl    | 11 +++++------
 src/DataWrangling/JRA55.jl                            |  1 +
 3 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/Project.toml b/Project.toml
index ff81ec6c..33cb4757 100644
--- a/Project.toml
+++ b/Project.toml
@@ -38,9 +38,9 @@ Downloads = "1.6"
 ImageMorphology = "0.4"
 JLD2 = "0.4, 0.5"
 KernelAbstractions = "0.9"
+Oceananigans = "0.93 - 0.94"
+OrthogonalSphericalShellGrids = "0.1.7"
 NCDatasets = "0.12, 0.13, 0.14"
-Oceananigans = "0.93.1"
-OrthogonalSphericalShellGrids = "0.1.2"
 Scratch = "1"
 SeawaterPolynomials = "0.3.4"
 StaticArrays = "1"
diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index 4d243b46..f43f4c19 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -12,16 +12,15 @@ z = exponential_z_faces(Nz=30, depth=6000)
 Nx = 120
 Ny = 60
 Nz = length(z) - 1
-
-grid = TripolarGrid(arch; z, size = (Nx, Ny, Nz), north_poles_latitude=55, first_pole_longitude=70)
-
-@show grid
+grid = TripolarGrid(arch; z, size=(Nx, Ny, Nz))
 
 bottom_height = regrid_bathymetry(grid;
                                   minimum_depth = 10,
                                   interpolation_passes = 5,
                                   major_basins = 3)
 
+grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
+
 # Closure
 gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000)
 catke = Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivity()
@@ -45,7 +44,7 @@ set!(restoring_mask_field, restoring_mask)
 Fu = Forcing(sponge_layer, field_dependencies=:u, parameters=restoring_rate)
 Fv = Forcing(sponge_layer, field_dependencies=:v, parameters=restoring_rate)
 
-dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1994, 1, 1)
+dates = DateTimeProlepticGregorian(1993, 11, 1) : Month(1) : DateTimeProlepticGregorian(1994, 11, 1)
 temperature = ECCOMetadata(:temperature, dates, ECCO4Monthly())
 salinity = ECCOMetadata(:salinity, dates, ECCO4Monthly())
 
@@ -68,7 +67,7 @@ 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=10minutes, stop_iteration=10)
 
 wall_time = Ref(time_ns())
 
diff --git a/src/DataWrangling/JRA55.jl b/src/DataWrangling/JRA55.jl
index df2ffce8..c02a3628 100644
--- a/src/DataWrangling/JRA55.jl
+++ b/src/DataWrangling/JRA55.jl
@@ -27,6 +27,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

From 746b5992e3dfe12b334e9e8cc1cff64104e3ad57 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Sun, 10 Nov 2024 17:30:12 +1100
Subject: [PATCH 02/30] don't relax compass

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 33cb4757..d9818525 100644
--- a/Project.toml
+++ b/Project.toml
@@ -38,7 +38,7 @@ Downloads = "1.6"
 ImageMorphology = "0.4"
 JLD2 = "0.4, 0.5"
 KernelAbstractions = "0.9"
-Oceananigans = "0.93 - 0.94"
+Oceananigans = "0.93.1 - 0.94"
 OrthogonalSphericalShellGrids = "0.1.7"
 NCDatasets = "0.12, 0.13, 0.14"
 Scratch = "1"

From 4ac5f21607e746ca24247afd66fad9b9f55df15e Mon Sep 17 00:00:00 2001
From: "Gregory L. Wagner" <wagner.greg@gmail.com>
Date: Sat, 9 Nov 2024 23:39:39 -0700
Subject: [PATCH 03/30] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index d9818525..cf3c2985 100644
--- a/Project.toml
+++ b/Project.toml
@@ -39,7 +39,7 @@ ImageMorphology = "0.4"
 JLD2 = "0.4, 0.5"
 KernelAbstractions = "0.9"
 Oceananigans = "0.93.1 - 0.94"
-OrthogonalSphericalShellGrids = "0.1.7"
+OrthogonalSphericalShellGrids = "0.1.8"
 NCDatasets = "0.12, 0.13, 0.14"
 Scratch = "1"
 SeawaterPolynomials = "0.3.4"

From 2736766699b94aff3ed12ead15a9a21b840f0560 Mon Sep 17 00:00:00 2001
From: "Gregory L. Wagner" <wagner.greg@gmail.com>
Date: Sat, 9 Nov 2024 23:41:52 -0700
Subject: [PATCH 04/30] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index cf3c2985..935ddead 100644
--- a/Project.toml
+++ b/Project.toml
@@ -38,7 +38,7 @@ Downloads = "1.6"
 ImageMorphology = "0.4"
 JLD2 = "0.4, 0.5"
 KernelAbstractions = "0.9"
-Oceananigans = "0.93.1 - 0.94"
+Oceananigans = "0.93.3 - 0.94"
 OrthogonalSphericalShellGrids = "0.1.8"
 NCDatasets = "0.12, 0.13, 0.14"
 Scratch = "1"

From 8545662f19456fae11cdd7979f21755f1efca270 Mon Sep 17 00:00:00 2001
From: "Gregory L. Wagner" <wagner.greg@gmail.com>
Date: Sat, 9 Nov 2024 23:45:34 -0700
Subject: [PATCH 05/30] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 935ddead..342c978f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -39,7 +39,7 @@ ImageMorphology = "0.4"
 JLD2 = "0.4, 0.5"
 KernelAbstractions = "0.9"
 Oceananigans = "0.93.3 - 0.94"
-OrthogonalSphericalShellGrids = "0.1.8"
+OrthogonalSphericalShellGrids = "0.1.7"
 NCDatasets = "0.12, 0.13, 0.14"
 Scratch = "1"
 SeawaterPolynomials = "0.3.4"

From e1481942be32eb29cf0cdd7e5cb1200460cfdc2d Mon Sep 17 00:00:00 2001
From: Gregory Wagner <wagner.greg@gmail.com>
Date: Sat, 9 Nov 2024 23:51:35 -0700
Subject: [PATCH 06/30] Generalize default free surface for IBG

---
 src/OceanSimulations/OceanSimulations.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl
index ea8f8296..a711d16d 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)

From 9c3aa7eb8d6d96ef42f3d0665b47bfcdf3d8f978 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Sun, 10 Nov 2024 19:44:28 +1100
Subject: [PATCH 07/30] Update Project.toml

---
 docs/Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/Project.toml b/docs/Project.toml
index ec9e24ea..64370ad1 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,3 +12,4 @@ CairoMakie = "0.10.12, 0.11, 0.12"
 DataDeps = "0.7"
 Documenter = "1"
 Oceananigans = "0.93.1"
+OrthogonalSphericalShellGrids = "0.1.7"

From ce7c11636a997b3ee0a9e5def4b527e1cbe3ee2a Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Sun, 10 Nov 2024 04:26:31 -0500
Subject: [PATCH 08/30] use interior for max of field

---
 .../one_degree_simulation/one_degree_simulation.jl  | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index f43f4c19..ec0eef52 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -67,7 +67,7 @@ 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=10minutes, stop_iteration=10)
+simulation = Simulation(coupled_model; Δt=30minutes, stop_time=2*365days)
 
 wall_time = Ref(time_ns())
 
@@ -75,9 +75,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",
@@ -89,6 +92,6 @@ function progress(sim)
      return nothing
 end
 
-add_callback!(simulation, progress, IterationInterval(1))
+add_callback!(simulation, progress, IterationInterval(10))
 
 run!(simulation)

From 75098198660b5a273d2e8cd54b8379416d32869e Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Sun, 10 Nov 2024 05:09:17 -0500
Subject: [PATCH 09/30] =?UTF-8?q?smaller=20=CE=94t?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 experiments/one_degree_simulation/one_degree_simulation.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index ec0eef52..601a06fb 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -67,7 +67,7 @@ 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=30minutes, stop_time=2*365days)
+simulation = Simulation(coupled_model; Δt=15minutes, stop_time=2*365days)
 
 wall_time = Ref(time_ns())
 

From 870e57a605e7bd9beb8eb68357422e7f90b05d72 Mon Sep 17 00:00:00 2001
From: Gregory Wagner <wagner.greg@gmail.com>
Date: Mon, 11 Nov 2024 00:53:26 -0500
Subject: [PATCH 10/30] Implement inpainting saving

---
 .../one_degree_simulation.jl                  |  4 +-
 src/DataWrangling/ECCO/ECCO.jl                | 85 +++++++++++++++----
 src/DataWrangling/ECCO/ECCO_metadata.jl       |  2 +-
 src/DataWrangling/JRA55.jl                    |  1 +
 4 files changed, 74 insertions(+), 18 deletions(-)

diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index 601a06fb..de5e30ca 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -8,10 +8,10 @@ using Dates
 using Printf
 
 arch = CPU()
-z = exponential_z_faces(Nz=30, depth=6000)
 Nx = 120
 Ny = 60
-Nz = length(z) - 1
+Nz = 20
+z = exponential_z_faces(; Nz, depth=6000)
 grid = TripolarGrid(arch; z, size=(Nx, Ny, Nz))
 
 bottom_height = regrid_bathymetry(grid;
diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl
index c5b1f743..2e2d7c63 100644
--- a/src/DataWrangling/ECCO/ECCO.jl
+++ b/src/DataWrangling/ECCO/ECCO.jl
@@ -15,12 +15,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
@@ -120,43 +122,61 @@ end
 """
     ECCO_field(metadata::ECCOMetadata;
                architecture = CPU(),
+               inpainting = nothing, #NearestNeighborInpainting(Inf),
+               mask = nothing, #ECCO_mask(metadata, architecture),
                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.
+
+K
+- `mask`: the mask used to inpaint the field, see [`inpaint_mask!`](@ref).
+- `inpainting`: the inpainting algorithm, see [`inpaint_mask!`](@ref). Default: `NearestNeighborInpainting(Inf)`.
 """
 function ECCO_field(metadata::ECCOMetadata;
                     architecture = CPU(),
+                    inpainting = nothing, #NearestNeighborInpainting(Inf),
+                    mask = nothing, #ECCO_mask(metadata, architecture),
                     horizontal_halo = (7, 7))
 
+    !isnothing(inpainting) && isnothing(mask) &&
+        throw(ArgumentError("Must provide a mask to use inpainting=$inpainting."))
+
+    field = empty_ECCO_field(metadata; architecture, horizontal_halo)
+    inpainted_path = inpainted_metadata_path(metadata)
+
+    if !isnothing(inpainting) && isfile(inpainted_path)
+        file = jldopen(inpainted_path)
+        data = file["data"]
+        close(file)
+        copyto!(parent(field), data)
+        return field
+    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
@@ -170,12 +190,40 @@ 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), ")")
+
+        file = jldopen(inpainted_path, "w+")
+        file["data"] = on_architecture(CPU(), parent(field))
+        close(file)
+    end
+
     return field
 end
 
 # Fallback
 ECCO_field(var_name::Symbol; kw...) = ECCO_field(ECCOMetadata(var_name); kw...)
 
+function inpainted_metadata_filename(metadata::ECCOMetadata)
+    original_filename = metadata_filename(metadata)
+    without_extension = original_filename[1:end-3]
+    return "inpainted_" * without_extension * ".jld2"
+end
+
+inpainted_metadata_path(metadata::ECCOMetadata) = joinpath(metadata.dir, inpainted_metadata_filename(metadata))
+
+#=
 """
     inpainted_ECCO_field(metadata::ECCOMetadata;
                          architecture = CPU(),
@@ -222,8 +270,12 @@ function inpainted_ECCO_field(metadata::ECCOMetadata;
 end
 
 inpainted_ECCO_field(variable_name::Symbol; kw...) = inpainted_ECCO_field(ECCOMetadata(variable_name); kw...)
+=#
     
-function set!(field::DistributedField, ECCO_metadata::ECCOMetadata; kw...)
+function set!(field::DistributedField, ECCO_metadata::ECCOMetadata;
+              inpainting = NearestNeighborInpainting(Inf),
+              kw...)
+
     # Fields initialized from ECCO
     grid = field.grid
     arch = architecture(grid)
@@ -231,9 +283,9 @@ function set!(field::DistributedField, ECCO_metadata::ECCOMetadata; kw...)
 
     f_ECCO = if arch.local_rank == 0 # Make sure we read/write the file using only one core
         mask = ECCO_mask(ECCO_metadata, child_arch)
-        inpainted_ECCO_field(ECCO_metadata; mask, architecture = child_arch, kw...)
+        ECCO_field(ECCO_metadata; inpainting, mask, architecture=child_arch, kw...)
     else
-        empty_ECCO_field(ECCO_metadata; architecture = child_arch)
+        empty_ECCO_field(ECCO_metadata; architecture=child_arch)
     end
 
     barrier!(arch)
@@ -248,15 +300,18 @@ function set!(field::DistributedField, ECCO_metadata::ECCOMetadata; kw...)
     return field
 end
 
-function set!(field::Field, ECCO_metadata::ECCOMetadata; kw...)
+function set!(field::Field, ECCO_metadata::ECCOMetadata;
+              inpainting = NearestNeighborInpainting(Inf), kw...)
+
     grid = field.grid
     arch = architecture(grid)
     mask = ECCO_mask(ECCO_metadata, arch)
     
-    f = inpainted_ECCO_field(ECCO_metadata; mask, architecture=arch, kw...)
-    f_grid = Field(location(ECCO_metadata), grid)   
-    interpolate!(f_grid, f)
-    set!(field, f_grid)
+    native_field = ECCO_field(ECCO_metadata; inpainting, mask, architecture=arch, kw...)
+    on_grid = Field(location(ECCO_metadata), grid)   
+    #interpolate!(on_grid, native_field)
+    interpolate!(field, native_field)
+    #set!(field, on_grid)
 
     return field
 end
diff --git a/src/DataWrangling/ECCO/ECCO_metadata.jl b/src/DataWrangling/ECCO/ECCO_metadata.jl
index ffd2c87c..4ff9b8ef 100644
--- a/src/DataWrangling/ECCO/ECCO_metadata.jl
+++ b/src/DataWrangling/ECCO/ECCO_metadata.jl
@@ -20,7 +20,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/JRA55.jl b/src/DataWrangling/JRA55.jl
index c02a3628..151ab36a 100644
--- a/src/DataWrangling/JRA55.jl
+++ b/src/DataWrangling/JRA55.jl
@@ -3,6 +3,7 @@ module JRA55
 using Oceananigans
 using Oceananigans.Units
  
+using Oceananigans: location
 using Oceananigans.Architectures: arch_array
 using Oceananigans.DistributedComputations
 using Oceananigans.DistributedComputations: child_architecture

From 854f5a80199199e5e3ddef59be09bbd4b64bacde Mon Sep 17 00:00:00 2001
From: Gregory Wagner <wagner.greg@gmail.com>
Date: Mon, 11 Nov 2024 00:58:12 -0500
Subject: [PATCH 11/30] Clean up

---
 src/DataWrangling/ECCO/ECCO.jl | 79 +++++++++-------------------------
 1 file changed, 21 insertions(+), 58 deletions(-)

diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl
index 2e2d7c63..002bfc98 100644
--- a/src/DataWrangling/ECCO/ECCO.jl
+++ b/src/DataWrangling/ECCO/ECCO.jl
@@ -122,22 +122,33 @@ end
 """
     ECCO_field(metadata::ECCOMetadata;
                architecture = CPU(),
-               inpainting = nothing, #NearestNeighborInpainting(Inf),
-               mask = nothing, #ECCO_mask(metadata, architecture),
-               horizontal_halo = (3, 3))
+               inpainting = nothing,
+               mask = nothing,
+               horizontal_halo = (7, 7))
 
-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.
+Return a `Field` on `architecture` described by `ECCOMetadata`.
+If not `nothing`, the `inpainting` method is used to fill the cells
+within the specified `mask`.
+
+Arguments
+=========
+
+- `metadata`
+
+Keyword arguments
+==================
+- `architecture`:
 
-K
 - `mask`: the mask used to inpaint the field, see [`inpaint_mask!`](@ref).
+
 - `inpainting`: the inpainting algorithm, see [`inpaint_mask!`](@ref). Default: `NearestNeighborInpainting(Inf)`.
+
+- `horizontal_halo`:
 """
 function ECCO_field(metadata::ECCOMetadata;
                     architecture = CPU(),
-                    inpainting = nothing, #NearestNeighborInpainting(Inf),
-                    mask = nothing, #ECCO_mask(metadata, architecture),
+                    inpainting = nothing,
+                    mask = nothing,
                     horizontal_halo = (7, 7))
 
     !isnothing(inpainting) && isnothing(mask) &&
@@ -223,55 +234,6 @@ end
 
 inpainted_metadata_path(metadata::ECCOMetadata) = joinpath(metadata.dir, inpainted_metadata_filename(metadata))
 
-#=
-"""
-    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
-end
-
-inpainted_ECCO_field(variable_name::Symbol; kw...) = inpainted_ECCO_field(ECCOMetadata(variable_name); kw...)
-=#
-    
 function set!(field::DistributedField, ECCO_metadata::ECCOMetadata;
               inpainting = NearestNeighborInpainting(Inf),
               kw...)
@@ -317,3 +279,4 @@ function set!(field::Field, ECCO_metadata::ECCOMetadata;
 end
 
 end # Module 
+

From f6d4b480225d2b0882ecade96e337bd8060f925a Mon Sep 17 00:00:00 2001
From: Gregory Wagner <wagner.greg@gmail.com>
Date: Mon, 11 Nov 2024 01:15:10 -0500
Subject: [PATCH 12/30] Test inpainting caching

---
 src/DataWrangling/ECCO/ECCO.jl |  8 +++--
 test/test_ecco.jl              | 59 ++++++++++++++++++----------------
 2 files changed, 37 insertions(+), 30 deletions(-)

diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl
index 002bfc98..ff0d28f0 100644
--- a/src/DataWrangling/ECCO/ECCO.jl
+++ b/src/DataWrangling/ECCO/ECCO.jl
@@ -151,8 +151,10 @@ function ECCO_field(metadata::ECCOMetadata;
                     mask = nothing,
                     horizontal_halo = (7, 7))
 
-    !isnothing(inpainting) && isnothing(mask) &&
-        throw(ArgumentError("Must provide a mask to use inpainting=$inpainting."))
+    # 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)
@@ -229,7 +231,7 @@ ECCO_field(var_name::Symbol; kw...) = ECCO_field(ECCOMetadata(var_name); kw...)
 function inpainted_metadata_filename(metadata::ECCOMetadata)
     original_filename = metadata_filename(metadata)
     without_extension = original_filename[1:end-3]
-    return "inpainted_" * without_extension * ".jld2"
+    return without_extension * "_inpainted.jld2"
 end
 
 inpainted_metadata_path(metadata::ECCOMetadata) = joinpath(metadata.dir, inpainted_metadata_filename(metadata))
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))

From 6d547ce509702c27f317e66df303357ee8200112 Mon Sep 17 00:00:00 2001
From: Gregory Wagner <wagner.greg@gmail.com>
Date: Mon, 11 Nov 2024 01:18:56 -0500
Subject: [PATCH 13/30] Update ECCO_field docstring

---
 src/DataWrangling/ECCO/ECCO.jl | 21 ++++-----------------
 1 file changed, 4 insertions(+), 17 deletions(-)

diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl
index ff0d28f0..319c8d4a 100644
--- a/src/DataWrangling/ECCO/ECCO.jl
+++ b/src/DataWrangling/ECCO/ECCO.jl
@@ -126,24 +126,11 @@ end
                mask = nothing,
                horizontal_halo = (7, 7))
 
-Return a `Field` on `architecture` described by `ECCOMetadata`.
+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`.
-
-Arguments
-=========
-
-- `metadata`
-
-Keyword arguments
-==================
-- `architecture`:
-
-- `mask`: the mask used to inpaint the field, see [`inpaint_mask!`](@ref).
-
-- `inpainting`: the inpainting algorithm, see [`inpaint_mask!`](@ref). Default: `NearestNeighborInpainting(Inf)`.
-
-- `horizontal_halo`:
+within the specified `mask`. `mask` is set to `ECCO_mask` for non-nothing
+`inpainting`.
 """
 function ECCO_field(metadata::ECCOMetadata;
                     architecture = CPU(),

From 70084390905093bf6d91c7e8d39ef5512c3c6025 Mon Sep 17 00:00:00 2001
From: sbishnu <sbishnu@eofe7.mit.edu>
Date: Wed, 13 Nov 2024 10:26:31 -0500
Subject: [PATCH 14/30] Modify one degree simulation

---
 .../one_degree_simulation.jl                  | 38 +++++++------------
 1 file changed, 13 insertions(+), 25 deletions(-)

diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index de5e30ca..3242c67a 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -7,12 +7,13 @@ using CFTime
 using Dates
 using Printf
 
-arch = CPU()
-Nx = 120
-Ny = 60
-Nz = 20
-z = exponential_z_faces(; Nz, depth=6000)
+arch = GPU()
+Nx = 360
+Ny = 180
+Nz = 60
+z = exponential_z_faces(; Nz, depth=5000)
 grid = TripolarGrid(arch; z, size=(Nx, Ny, Nz))
+@info grid
 
 bottom_height = regrid_bathymetry(grid;
                                   minimum_depth = 10,
@@ -24,33 +25,20 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
 # 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
+viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=1000)
+closure = (gm, catke, viscous_closure)
 
 restoring_rate = 1 / 1days
 
-restoring_mask_field = CenterField(grid)
-set!(restoring_mask_field, restoring_mask)
-
-@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, 80), z=(-10, 0))
 
 dates = DateTimeProlepticGregorian(1993, 11, 1) : Month(1) : DateTimeProlepticGregorian(1994, 11, 1)
 temperature = ECCOMetadata(:temperature, dates, ECCO4Monthly())
 salinity = ECCOMetadata(:salinity, dates, ECCO4Monthly())
 
-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)
+FT = ECCORestoring(arch, temperature; grid, mask, rate=restoring_rate)
+FS = ECCORestoring(arch, salinity;    grid, mask, rate=restoring_rate)
+forcing = (T=FT, S=FS)
 
 momentum_advection = VectorInvariant()
 tracer_advection = Centered(order=2)
@@ -63,7 +51,7 @@ set!(ocean.model,
      S = ECCOMetadata(:salinity; dates=first(dates)))
 
 radiation = Radiation(arch)
-atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(41))
+atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20))
 sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce()
 coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
 

From 13091c999a5e100bcffab40d6335b8eb453117f6 Mon Sep 17 00:00:00 2001
From: Simone Silvestri <silvestri.simone0@gmail.com>
Date: Wed, 13 Nov 2024 18:06:07 +0100
Subject: [PATCH 15/30] some more updates to the one degree simulation

---
 .../one_degree_simulation.jl                  | 94 ++++++++++++++-----
 1 file changed, 69 insertions(+), 25 deletions(-)

diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index 3242c67a..39084342 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -1,62 +1,106 @@
 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
+
+arch = CPU()
+
+#####
+##### Grid and Bathymetry
+#####
 
-arch = GPU()
 Nx = 360
 Ny = 180
-Nz = 60
-z = exponential_z_faces(; Nz, depth=5000)
-grid = TripolarGrid(arch; z, size=(Nx, Ny, Nz))
-@info grid
+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(grid;
+bottom_height = regrid_bathymetry(underlying_grid;
                                   minimum_depth = 10,
-                                  interpolation_passes = 5,
-                                  major_basins = 3)
+                                  interpolation_passes = 75,
+                                  major_basins = 2)
 
-grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
+# 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()
-viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=1000)
+catke = ClimaOcean.OceanSimulations.default_ocean_closure()
+viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=10000)
+
 closure = (gm, catke, viscous_closure)
 
-restoring_rate = 1 / 1days
+#####
+##### Restoring
+#####
 
-mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 80), z=(-10, 0))
+restoring_rate  = 1 / 2days
+z_below_surface = @allowscalar znode(1, 1, grid.Nz, grid, Center(), Center(), Face())
 
-dates = DateTimeProlepticGregorian(1993, 11, 1) : Month(1) : DateTimeProlepticGregorian(1994, 11, 1)
-temperature = ECCOMetadata(:temperature, dates, ECCO4Monthly())
-salinity = ECCOMetadata(:salinity, dates, ECCO4Monthly())
+mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(z_below_surface, 0))
 
-FT = ECCORestoring(arch, temperature; grid, mask, rate=restoring_rate)
-FS = ECCORestoring(arch, salinity;    grid, mask, rate=restoring_rate)
+dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1993, 11, 1)
+temperature = ECCOMetadata(:temperature; dates, version=ECCO4Monthly(), dir="./")
+salinity    = ECCOMetadata(:salinity;    dates, version=ECCO4Monthly(), dir="./")
+
+# inpainting = NearestNeighborInpainting(30) should be enough to fill the gaps near bathymetry
+FT = ECCORestoring(arch, temperature; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(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)
+radiation  = Radiation(arch)
 atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20))
+
+#####
+##### Coupled simulation
+#####
+
 sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce()
 coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
 
 simulation = Simulation(coupled_model; Δt=15minutes, stop_time=2*365days)
 
+#####
+##### Run it!
+##### 
+
 wall_time = Ref(time_ns())
 
 function progress(sim)

From 60b8a259036d857fb6227685100ff394782da018 Mon Sep 17 00:00:00 2001
From: Simone Silvestri <silvestri.simone0@gmail.com>
Date: Wed, 13 Nov 2024 18:18:39 +0100
Subject: [PATCH 16/30] bugfix

---
 experiments/one_degree_simulation/one_degree_simulation.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index 39084342..f32b9cd0 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -18,7 +18,7 @@ Nx = 360
 Ny = 180
 Nz = 100
 
-z_faces = exponential_z_faces(; Nz, depth=5000, h=34),
+z_faces = exponential_z_faces(; Nz, depth=5000, h=34)
 
 underlying_grid = TripolarGrid(arch;
                                size = (Nx, Ny, Nz),

From eef5a1ba43b31edc73e8070a9c1a4d84d223e1a7 Mon Sep 17 00:00:00 2001
From: Simone Silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 11:05:11 +0100
Subject: [PATCH 17/30] cache inpainted data?

---
 src/DataWrangling/ECCO/ECCO.jl           | 28 ++++++++-------
 src/DataWrangling/ECCO/ECCO_restoring.jl | 43 +++++++++++++++++-------
 2 files changed, 47 insertions(+), 24 deletions(-)

diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl
index c9e15e1e..c6f5c486 100644
--- a/src/DataWrangling/ECCO/ECCO.jl
+++ b/src/DataWrangling/ECCO/ECCO.jl
@@ -126,7 +126,8 @@ end
                architecture = CPU(),
                inpainting = nothing,
                mask = nothing,
-               horizontal_halo = (7, 7))
+               horizontal_halo = (7, 7),
+               cache_inpainted_data = false)
 
 Return a `Field` on `architecture` described by `ECCOMetadata` with
 `horizontal_halo` size.
@@ -138,7 +139,8 @@ function ECCO_field(metadata::ECCOMetadata;
                     architecture = CPU(),
                     inpainting = nothing,
                     mask = nothing,
-                    horizontal_halo = (7, 7))
+                    horizontal_halo = (7, 7),
+                    cache_inpainted_data = false)
 
     # Respect user-supplied mask, but otherwise build default ECCO mask.
     if !isnothing(inpainting) && isnothing(mask)
@@ -150,10 +152,15 @@ function ECCO_field(metadata::ECCOMetadata;
 
     if !isnothing(inpainting) && isfile(inpainted_path)
         file = jldopen(inpainted_path)
-        data = file["data"]
-        close(file)
-        copyto!(parent(field), data)
-        return field
+        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)
@@ -206,10 +213,10 @@ function ECCO_field(metadata::ECCOMetadata;
         elapsed = 1e-9 * (time_ns() - start_time)
         @info string(" ... (", prettytime(elapsed), ")")
     
-        if save_inpainted
+        if cache_inpainted_data
             file = jldopen(inpainted_path, "w+")
             file["data"] = on_architecture(CPU(), parent(field))
-            file["inpainting_steps"] = inpainting.steps
+            file["inpainting_maxiter"] = inpainting.maxiter
             close(file)
         end
     end
@@ -228,9 +235,7 @@ end
 
 inpainted_metadata_path(metadata::ECCOMetadata) = joinpath(metadata.dir, inpainted_metadata_filename(metadata))
 
-function set!(field::DistributedField, ECCO_metadata::ECCOMetadata;
-              inpainting = NearestNeighborInpainting(Inf),
-              kw...)
+function set!(field::Field, ECCO_metadata::ECCOMetadata; kw...)
 
     # Fields initialized from ECCO
     grid = field.grid
@@ -239,7 +244,6 @@ function set!(field::DistributedField, ECCO_metadata::ECCOMetadata;
 
     f = ECCO_field(ecco_metadata; mask,
                    architecture = arch,
-                   inpainting,
                    kw...)
 
     interpolate!(field, f)
diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index c8f52d44..86379ae9 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
         M = typeof(metadata)
         I = typeof(inpainting)
-        return new{N, I, M}(start, length, inpainting, metadata)
+        return new{N, I, M}(start, length, inpainting_settings, 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: `false`.
+
 """
 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 = false,
                                 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 = false)
 
 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
@@ -292,6 +301,9 @@ Keyword Arguments
           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: `false`.
+                          
 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.
 """
@@ -312,9 +324,16 @@ function ECCORestoring(architecture::AbstractArchitecture,
                        grid = nothing,
                        time_indices_in_memory = 2, # Not more than this if we want to use GPU!
                        time_indexing = Cyclical(),
-                       inpainting = NearestNeighborInpainting(Inf))
-
-    fts = ECCO_field_time_series(metadata; grid, architecture, time_indices_in_memory, time_indexing, inpainting)
+                       inpainting = NearestNeighborInpainting(Inf),
+                       cache_inpainted_data = false)
+
+    fts = ECCO_field_time_series(metadata; 
+                                 grid, 
+                                 architecture, 
+                                 time_indices_in_memory, 
+                                 time_indexing, 
+                                 inpainting,
+                                 cache_inpainted_data)
 
     # Grab the correct Oceananigans field to restore
     variable_name = metadata.name

From c6532999f588052f0c2d58f81c831f586ecdee75 Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 05:14:41 -0500
Subject: [PATCH 18/30] ebugfix

---
 src/DataWrangling/ECCO/ECCO_restoring.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index 86379ae9..2642c07a 100644
--- a/src/DataWrangling/ECCO/ECCO_restoring.jl
+++ b/src/DataWrangling/ECCO/ECCO_restoring.jl
@@ -28,7 +28,7 @@ struct ECCONetCDFBackend{N, C, I, M} <: AbstractInMemoryBackend{Int}
     function ECCONetCDFBackend{N, C}(start::Int, length::Int, inpainting, metadata) where N
         M = typeof(metadata)
         I = typeof(inpainting)
-        return new{N, I, M}(start, length, inpainting_settings, metadata)
+        return new{N, C, I, M}(start, length, inpainting_settings, metadata)
     end
 end
 

From d1dd2f3d3a1ab4c64c00c5539b35b88b3d8d44a4 Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 06:29:05 -0500
Subject: [PATCH 19/30] advancing

---
 experiments/one_degree_simulation/one_degree_simulation.jl | 2 ++
 src/DataWrangling/ECCO/ECCO_restoring.jl                   | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index f32b9cd0..77e5b80e 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -8,6 +8,8 @@ using Dates
 using Printf
 using CUDA: @allowscalar
 
+using Oceananigans.Grids: znode
+
 arch = CPU()
 
 #####
diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index 2642c07a..dae1e999 100644
--- a/src/DataWrangling/ECCO/ECCO_restoring.jl
+++ b/src/DataWrangling/ECCO/ECCO_restoring.jl
@@ -25,7 +25,7 @@ struct ECCONetCDFBackend{N, C, I, M} <: AbstractInMemoryBackend{Int}
     inpainting :: I
     metadata :: M
 
-    function ECCONetCDFBackend{N, C}(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, C, I, M}(start, length, inpainting_settings, metadata)

From b1ae2ae1e5c7efeef11023ec6d68c7148b7c6a18 Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 06:30:24 -0500
Subject: [PATCH 20/30] small bug

---
 src/DataWrangling/ECCO/ECCO_restoring.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index dae1e999..4bb689d4 100644
--- a/src/DataWrangling/ECCO/ECCO_restoring.jl
+++ b/src/DataWrangling/ECCO/ECCO_restoring.jl
@@ -28,7 +28,7 @@ struct ECCONetCDFBackend{N, C, I, M} <: AbstractInMemoryBackend{Int}
     function ECCONetCDFBackend{N, C}(start::Int, length::Int, inpainting, metadata) where {N, C}
         M = typeof(metadata)
         I = typeof(inpainting)
-        return new{N, C, I, M}(start, length, inpainting_settings, metadata)
+        return new{N, C, I, M}(start, length, inpainting, metadata)
     end
 end
 

From ccbe516549cbd9a1fd2bc280f4962e880d2e1e3d Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 08:56:03 -0500
Subject: [PATCH 21/30] some bugfixes

---
 experiments/one_degree_simulation/one_degree_simulation.jl | 5 +++--
 src/DataWrangling/ECCO/ECCO.jl                             | 4 ++--
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
index 77e5b80e..ef7b558a 100644
--- a/experiments/one_degree_simulation/one_degree_simulation.jl
+++ b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -6,11 +6,12 @@ using Oceananigans.Units
 using CFTime
 using Dates
 using Printf
-using CUDA: @allowscalar
+using CUDA: @allowscalar, device!
 
 using Oceananigans.Grids: znode
 
-arch = CPU()
+device!(3)
+arch = GPU()
 
 #####
 ##### Grid and Bathymetry
diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl
index c6f5c486..fca755d8 100644
--- a/src/DataWrangling/ECCO/ECCO.jl
+++ b/src/DataWrangling/ECCO/ECCO.jl
@@ -240,9 +240,9 @@ 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 = ECCO_field(ecco_metadata; mask,
+    f = ECCO_field(ECCO_metadata; mask,
                    architecture = arch,
                    kw...)
 

From bccd3a8be91be85d82bff3ceadcb310feca335cf Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 09:39:30 -0500
Subject: [PATCH 22/30] add a `@root` to caching

---
 src/DataWrangling/ECCO/ECCO.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl
index fca755d8..78c09309 100644
--- a/src/DataWrangling/ECCO/ECCO.jl
+++ b/src/DataWrangling/ECCO/ECCO.jl
@@ -213,7 +213,8 @@ function ECCO_field(metadata::ECCOMetadata;
         elapsed = 1e-9 * (time_ns() - start_time)
         @info string(" ... (", prettytime(elapsed), ")")
     
-        if cache_inpainted_data
+        # 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

From f1d3ceda81fb90a657453bf91a404d8f4284b2e2 Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 09:40:11 -0500
Subject: [PATCH 23/30] open with read permission

---
 src/DataWrangling/ECCO/ECCO.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl
index 78c09309..bff90c60 100644
--- a/src/DataWrangling/ECCO/ECCO.jl
+++ b/src/DataWrangling/ECCO/ECCO.jl
@@ -151,7 +151,7 @@ function ECCO_field(metadata::ECCOMetadata;
     inpainted_path = inpainted_metadata_path(metadata)
 
     if !isnothing(inpainting) && isfile(inpainted_path)
-        file = jldopen(inpainted_path)
+        file = jldopen(inpainted_path, "r")
         maxiter = file["inpainting_maxiter"]
 
         # read data if generated with the same inpainting

From 517782ca23a9622d394d0dfede20b008200d044c Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 10:14:05 -0500
Subject: [PATCH 24/30] interpolate on grid beforehand

---
 src/DataWrangling/ECCO/ECCO_restoring.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index 4bb689d4..9037a33b 100644
--- a/src/DataWrangling/ECCO/ECCO_restoring.jl
+++ b/src/DataWrangling/ECCO/ECCO_restoring.jl
@@ -293,8 +293,6 @@ 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.
@@ -339,7 +337,9 @@ function ECCORestoring(architecture::AbstractArchitecture,
     variable_name = metadata.name
     field_name = oceananigans_fieldname[variable_name]
 
-    return ECCORestoring(fts, fts.grid, mask, field_name, rate)
+    fts_grid = isnothing(grid) ? nothing : fts.grid 
+
+    return ECCORestoring(fts, fts_grid, mask, field_name, rate)
 end
 
 # Make sure we can call ECCORestoring with architecture as the first positional argument

From 353966004de81c694a5964583a9feee37ae3034d Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 10:15:26 -0500
Subject: [PATCH 25/30] fix minimum-seaice in the meantime

---
 src/OceanSeaIceModels/minimum_temperature_sea_ice.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl b/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl
index 09e2c155..07dc718d 100644
--- a/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl
+++ b/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl
@@ -43,9 +43,10 @@ end
                                          ocean_temperature)
 
     i, j = @index(Global, NTuple)
-
+    kᴺ = size(grid, 3) 
+    
     @inbounds begin
-        Tₒ = ocean_temperature[i, j, 1]
+        Tₒ = ocean_temperature[i, j, kᴺ]
 
         τx = centered_velocity_fluxes.u
         τy = centered_velocity_fluxes.v

From 235a67c7901679d2fa527763c6da919808270cdf Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 10:16:49 -0500
Subject: [PATCH 26/30] comment

---
 src/DataWrangling/ECCO/ECCO_restoring.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index 9037a33b..e0955591 100644
--- a/src/DataWrangling/ECCO/ECCO_restoring.jl
+++ b/src/DataWrangling/ECCO/ECCO_restoring.jl
@@ -337,9 +337,11 @@ function ECCORestoring(architecture::AbstractArchitecture,
     variable_name = metadata.name
     field_name = oceananigans_fieldname[variable_name]
 
-    fts_grid = isnothing(grid) ? nothing : fts.grid 
+    # 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, fts_grid, mask, field_name, rate)
+    return ECCORestoring(fts, native_grid, mask, field_name, rate)
 end
 
 # Make sure we can call ECCORestoring with architecture as the first positional argument

From 5c23b6e723cadbaa19408b32faeeb6d97afc360e Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 10:18:45 -0500
Subject: [PATCH 27/30] trhow an error

---
 src/DataWrangling/ECCO/ECCO_restoring.jl | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index e0955591..85a9faf7 100644
--- a/src/DataWrangling/ECCO/ECCO_restoring.jl
+++ b/src/DataWrangling/ECCO/ECCO_restoring.jl
@@ -315,7 +315,7 @@ function ECCORestoring(arch::AbstractArchitecture,
     return ECCORestoring(arch, metadata; kw...)
 end
 
-function ECCORestoring(architecture::AbstractArchitecture,
+function ECCORestoring(arch::AbstractArchitecture,
                        metadata::ECCOMetadata;
                        rate,
                        mask = 1,
@@ -325,9 +325,14 @@ function ECCORestoring(architecture::AbstractArchitecture,
                        inpainting = NearestNeighborInpainting(Inf),
                        cache_inpainted_data = false)
 
+    # 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,

From bcfc307e6396065eb55a041d95400e973d9a6786 Mon Sep 17 00:00:00 2001
From: simone-silvestri <silvestri.simone0@gmail.com>
Date: Fri, 15 Nov 2024 10:59:57 -0500
Subject: [PATCH 28/30] correct arch

---
 src/DataWrangling/ECCO/ECCO_restoring.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index 85a9faf7..0757b0c9 100644
--- a/src/DataWrangling/ECCO/ECCO_restoring.jl
+++ b/src/DataWrangling/ECCO/ECCO_restoring.jl
@@ -332,7 +332,7 @@ function ECCORestoring(arch::AbstractArchitecture,
 
     fts = ECCO_field_time_series(metadata; 
                                  grid, 
-                                 arch, 
+                                 architecture = arch, 
                                  time_indices_in_memory, 
                                  time_indexing, 
                                  inpainting,

From 88cc7e02c7c710d50e9ca1ad6793cafbb9ed4462 Mon Sep 17 00:00:00 2001
From: Gregory Wagner <wagner.greg@gmail.com>
Date: Mon, 18 Nov 2024 08:28:14 -0700
Subject: [PATCH 29/30] rm fake sea ice again

---
 .../minimum_temperature_sea_ice.jl            | 67 -------------------
 1 file changed, 67 deletions(-)
 delete mode 100644 src/OceanSeaIceModels/minimum_temperature_sea_ice.jl

diff --git a/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl b/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl
deleted file mode 100644
index 07dc718d..00000000
--- a/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl
+++ /dev/null
@@ -1,67 +0,0 @@
-using Oceananigans.Architectures: architecture
-
-import ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: limit_fluxes_over_sea_ice!
-
-"""
-    struct MinimumTemperatureSeaIce{T}
-
-The minimal possible sea ice representation, providing an "Insulating layer" on the surface.
-Not really a ``model'' per se, however, it is the most simple way to make sure that temperature 
-does not dip below freezing temperature.
-All fluxes are shut down when the surface is below the `minimum_temperature` except for heating.
-
-# Fields
-- `minimum_temperature`: The minimum temperature of water.
-"""
-struct MinimumTemperatureSeaIce{T}
-    minimum_temperature :: T
-end
-
-MinimumTemperatureSeaIce() = MinimumTemperatureSeaIce(-1.8)
-
-function limit_fluxes_over_sea_ice!(grid, kernel_parameters,
-                                    sea_ice::MinimumTemperatureSeaIce,
-                                    centered_velocity_fluxes,
-                                    net_tracer_fluxes,
-                                    ocean_temperature,
-                                    ocean_salinity)
-
-    launch!(architecture(grid), grid, kernel_parameters, _cap_fluxes_on_sea_ice!,
-            centered_velocity_fluxes,
-            net_tracer_fluxes,
-            grid,
-            sea_ice.minimum_temperature,
-            ocean_temperature)
-
-    return nothing
-end
-
-@kernel function _cap_fluxes_on_sea_ice!(centered_velocity_fluxes,
-                                         net_tracer_fluxes,
-                                         grid,
-                                         minimum_temperature,
-                                         ocean_temperature)
-
-    i, j = @index(Global, NTuple)
-    kᴺ = size(grid, 3) 
-    
-    @inbounds begin
-        Tₒ = ocean_temperature[i, j, kᴺ]
-
-        τx = centered_velocity_fluxes.u
-        τy = centered_velocity_fluxes.v
-        Jᵀ = net_tracer_fluxes.T
-        Jˢ = net_tracer_fluxes.S
-
-        sea_ice = Tₒ < minimum_temperature
-        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]) 
-
-        # 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])
-        τx[i, j, 1] = ifelse(sea_ice, zero(grid), τx[i, j, 1])
-        τy[i, j, 1] = ifelse(sea_ice, zero(grid), τy[i, j, 1])
-    end
-end

From 5715bf4c9ffca7b3ed7a18efacd6805cacf310a7 Mon Sep 17 00:00:00 2001
From: Gregory Wagner <wagner.greg@gmail.com>
Date: Mon, 18 Nov 2024 08:35:36 -0700
Subject: [PATCH 30/30] Make default true for caching

---
 src/DataWrangling/ECCO/ECCO_restoring.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl
index bec7b1bf..a19916bd 100644
--- a/src/DataWrangling/ECCO/ECCO_restoring.jl
+++ b/src/DataWrangling/ECCO/ECCO_restoring.jl
@@ -69,7 +69,7 @@ function set!(fts::ECCONetCDFFTS)
     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, cache_inpainted_data = cache_data)
+        set!(fts[t], metadatum; inpainting, cache_inpainted_data=cache_data)
     end
 
     fill_halo_regions!(fts)
@@ -135,7 +135,7 @@ Keyword Arguments
                 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: `false`.
+                          Default: `true`.
 
 """
 function ECCO_field_time_series(metadata::ECCOMetadata;	
@@ -143,7 +143,7 @@ function ECCO_field_time_series(metadata::ECCOMetadata;
                                 time_indices_in_memory = 2,	
                                 time_indexing = Cyclical(),
                                 inpainting = NearestNeighborInpainting(prod(size(metadata))),
-                                cache_inpainted_data = false,
+                                cache_inpainted_data = true,
                                 grid = nothing)	
 
     # Make sure all the required individual files are downloaded
@@ -251,7 +251,7 @@ end
                   rate = 1,
                   grid = nothing,
                   inpainting = NearestNeighborInpainting(prod(size(metadata))),
-                  cache_inpainted_data = false)
+                  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
@@ -300,7 +300,7 @@ Keyword Arguments
           Default: nothing.
 
 - `cache_inpainted_data`: If `true`, the data is cached to disk after inpainting for later retrieving. 
-                          Default: `false`.
+                          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.
@@ -323,7 +323,7 @@ function ECCORestoring(arch::AbstractArchitecture,
                        time_indices_in_memory = 2, # Not more than this if we want to use GPU!
                        time_indexing = Cyclical(),
                        inpainting = NearestNeighborInpainting(Inf),
-                       cache_inpainted_data = false)
+                       cache_inpainted_data = true)
 
     # Validate architecture
     if !isnothing(grid) && architecture(grid) != arch