More one degree simulation (#260)

* one change

* new changes

* close file

* Inf inpainting in the inpaint mask

* correct stuff

* last bugfix

* revert back

* revert back

* revert back

* remove circular dependency of mask and field

* bugfix

* last bugfix in mask

* another bugfix

* fix mask

experiments/one_degree_simulation/one_degree_simulation.jl

@@ -10,7 +10,6 @@ using CUDA: @allowscalar, device!
 
 using Oceananigans.Grids: znode
 
-device!(3)
 arch = GPU()
 
 #####
@@ -47,7 +46,7 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_he
 
 gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000)
 catke = ClimaOcean.OceanSimulations.default_ocean_closure()
-viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=10000)
+viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=2000)
 
 closure = (gm, catke, viscous_closure)
 
@@ -65,8 +64,8 @@ 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))
+FT = ECCORestoring(arch, temperature; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(50))
+FS = ECCORestoring(arch, salinity;    grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(50))
 forcing = (T=FT, S=FS)
 
 #####
@@ -95,9 +94,7 @@ atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20))
 ##### Coupled simulation
 #####
 
-sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce()
-coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
-
+coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) 
 simulation = Simulation(coupled_model; Δt=15minutes, stop_time=2*365days)
 
 #####

src/DataWrangling/ECCO/ECCO.jl

@@ -6,7 +6,7 @@ export ECCORestoring, LinearlyTaperedPolarMask
 
 using ClimaOcean
 using ClimaOcean.DataWrangling
-using ClimaOcean.DataWrangling: inpaint_mask!
+using ClimaOcean.DataWrangling: inpaint_mask!, NearestNeighborInpainting
 using ClimaOcean.InitialConditions: three_dimensional_regrid!, interpolate!
 
 using Oceananigans
@@ -138,16 +138,11 @@ within the specified `mask`. `mask` is set to `ECCO_mask` for non-nothing
 """
 function ECCO_field(metadata::ECCOMetadata;
                     architecture = CPU(),
-                    inpainting = nothing,
+                    inpainting = NearestNeighborInpainting(Inf),
                     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)
 
@@ -162,6 +157,8 @@ function ECCO_field(metadata::ECCOMetadata;
             copyto!(parent(field), data)
             return field
         end
+
+        close(file)
     end
 
     download_dataset(metadata)
@@ -201,6 +198,11 @@ function ECCO_field(metadata::ECCOMetadata;
     fill_halo_regions!(field)
 
     if !isnothing(inpainting)
+        # Respect user-supplied mask, but otherwise build default ECCO mask.
+        if isnothing(mask)
+            mask = ECCO_mask(metadata, architecture; data_field=field)
+        end
+
         # Make sure all values are extended properly
         name = string(metadata.name)
         date = string(metadata.dates)

src/DataWrangling/ECCO/ECCO_mask.jl

@@ -10,11 +10,11 @@ import ClimaOcean: stateindex
 A boolean field where `true` represents a missing value in the ECCO dataset.
 """
 function ECCO_mask(metadata, architecture = CPU(); 
+                   data_field = ECCO_field(metadata; architecture, inpainting=nothing),
                    minimum_value = Float32(-1e5),
                    maximum_value = Float32(1e5))
 
-    field = ECCO_field(metadata; architecture)
-    mask  = Field{location(field)...}(field.grid, Bool)
+    mask  = Field{location(data_field)...}(data_field.grid, Bool)
 
     # ECCO4 has zeros in place of the missing values, while
     # ECCO2 expresses missing values with values < -1e5
@@ -25,7 +25,7 @@ function ECCO_mask(metadata, architecture = CPU();
     end
 
     # Set the mask with zeros where field is defined
-    launch!(architecture, field.grid, :xyz, _set_mask!, mask, field, minimum_value, maximum_value)
+    launch!(architecture, data_field.grid, :xyz, _set_mask!, mask, data_field, minimum_value, maximum_value)
 
     return mask
 end

src/DataWrangling/ECCO/ECCO_restoring.jl

@@ -54,7 +54,7 @@ Base.summary(backend::ECCONetCDFBackend) = string("ECCONetCDFBackend(", backend.
 const ECCONetCDFFTS{N} = FlavorOfFTS{<:Any, <:Any, <:Any, <:Any, <:ECCONetCDFBackend{N}} where N
 
 new_backend(b::ECCONetCDFBackend{native, cache_data}, start, length) where {native, cache_data} =
-    ECCONetCDFBackend{native}(start, length, b.inpainting, b.metadata)
+    ECCONetCDFBackend{native, cache_data}(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

src/DataWrangling/inpaint_mask.jl

@@ -140,7 +140,7 @@ Arguments
 - `inpainting`: The inpainting algorithm to use. For the moment, the only option is `NearestNeighborInpainting(maxiter)`, 
                 where an average of the valid surrounding values is used `maxiter` times.
 """
-function inpaint_mask!(field, mask; inpainting=NearestNeighborInpainting(10))
+function inpaint_mask!(field, mask; inpainting=NearestNeighborInpainting(Inf))
 
     if inpainting isa Int
         inpainting = NearestNeighborInpainting(inpainting)