Do not restrict bathymetry (#335)

* Update Bathymetry.jl

* add a fill halo

* go ahead

* the bathymetry is supposed to be on centers not on nodes

* this should work

* bugfix

Author: simone-silvestri

src/Bathymetry.jl

@@ -126,29 +126,10 @@ function regrid_bathymetry(target_grid;
     close(dataset)
 
     # Diagnose target grid information
-    arch = architecture(target_grid)
-    φ₁, φ₂ = y_domain(target_grid)
-    λ₁, λ₂ = x_domain(target_grid)
-
-    if λ₁ < 0 || λ₂ > 360
-        throw(ArgumentError("Cannot regrid bathymetry between λ₁ = $(λ₁) and λ₂ = $(λ₂).
-                             Bathymetry data is defined on longitudes spanning λ = (0, 360)."))
-    end
-
-    # Calculate limiting indices on the bathymetry grid
-    i₁ = searchsortedfirst(λ_data, λ₁)
-    i₂ = searchsortedfirst(λ_data, λ₂) - 1
-    ii = i₁:i₂
-
-    j₁ = searchsortedfirst(φ_data, φ₁)
-    j₂ = searchsortedfirst(φ_data, φ₂) - 1
-    jj = j₁:j₂
-
-    # Restrict bathymetry coordinate_data to region of interest
-    λ_data = λ_data[ii] |> Array{BigFloat}
-    φ_data = φ_data[jj] |> Array{BigFloat}
-    z_data = z_data[ii, jj] 
-
+    arch   = architecture(target_grid)
+    λ_data = λ_data |> Array{BigFloat}
+    φ_data = φ_data |> Array{BigFloat}
+    
     if !isnothing(height_above_water)
         # Overwrite the height of cells above water.
         # This has an impact on reconstruction. Greater height_above_water reduces total
@@ -170,6 +151,7 @@ function regrid_bathymetry(target_grid;
     Nyn = length(φ_data)
     Nzn = 1
 
+  
     native_grid = LatitudeLongitudeGrid(arch, Float32;
                                         size = (Nxn, Nyn, Nzn),
                                         latitude  = (φ₁_data, φ₂_data),
@@ -179,7 +161,8 @@ function regrid_bathymetry(target_grid;
 
     native_z = Field{Center, Center, Nothing}(native_grid)
     set!(native_z, z_data)
-
+    fill_halo_regions!(native_z)
+    
     target_z = interpolate_bathymetry_in_passes(native_z, target_grid; 
                                                 passes = interpolation_passes)
 
@@ -206,9 +189,15 @@ function interpolate_bathymetry_in_passes(native_z, target_grid;
     Nλt, Nφt = Nt = size(target_grid)
     Nλn, Nφn = Nn = size(native_z)
 
+    resxt = minimum_xspacing(target_grid)
+    resyt = minimum_yspacing(target_grid)
+
+    resxn = minimum_xspacing(native_z.grid)
+    resyn = minimum_yspacing(native_z.grid)
+
     # Check whether we are coarsening the grid in any directions.
     # If so, skip interpolation passes.
-    if Nλt > Nλn || Nφt > Nφn
+    if resxt > resxn || resyt > resyn
         target_z = Field{Center, Center, Nothing}(target_grid)
         interpolate!(target_z, native_z)
         @info string("Skipping passes for interpolating bathymetry of size $Nn ", '\n',

test/test_distributed_utils.jl

@@ -78,9 +78,9 @@ end
 end
 
 @testset "Distributed Bathymetry interpolation" begin
-    # We start by building a fake bathyemtry on rank 0 and save it to file
+    # We start by building a fake bathymetry on rank 0 and save it to file
     @root begin
-        λ = -180:0.1:180
+        λ = -179.5:0.1:179.5
         φ = 0:0.1:50
 
         Nλ = length(λ)

test/test_simulations.jl

@@ -2,7 +2,7 @@ include("runtests_setup.jl")
 
 using CUDA
 using OrthogonalSphericalShellGrids
-using ClimaOcean.OceanSeaIceModels: adjust_freezing_ocean_temperature!
+using ClimaOcean.OceanSeaIceModels: above_freezing_ocean_temperature!
 using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
 
 @testset "GPU time stepping test" begin
@@ -53,7 +53,7 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
         # Set the ocean temperature and salinity
         set!(ocean.model, T=temperature_metadata[1], S=salinity_metadata[1])
 
-        adjust_freezing_ocean_temperature!(ocean, sea_ice)
+        above_freezing_ocean_temperature!(ocean, sea_ice)
 
         # Test that ocean temperatures are above freezing
         T = on_architecture(CPU(), ocean.model.T)