Grid cell average bug fixes (#673)

* grid cell average bug fixes * changelog updated
SpeedyWeather · Feb 11, 2025 · e43d3a7 · e43d3a7
1 parent 898f997
commit e43d3a7
Show file tree

Hide file tree

Showing 4 changed files with 30 additions and 4 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Unreleased
 
+- Fix grid_cell_average! bugs for some grids and resolutions [#673](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/673)
 - Interfaces for interpolation of AbstractGridArray [#671](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/671)
 - Test folder sorted into subfolders [#671](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/671)
 - Land model modularised + land netCDF output [#671](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/671)

diff --git a/src/RingGrids/interpolation.jl b/src/RingGrids/interpolation.jl
@@ -533,24 +533,28 @@ function grid_cell_average!(
         j1 = findlast( θ -> θ <  θ1, colat_in)
 
         for ij in ring
-            ϕ0 = mod(lons_out[ij] - Δϕ/2, 2π)        # western edge
+            ϕ0 = mod(lons_out[ij] - Δϕ/2, 2π)       # western edge
             ϕ1 = ϕ0 + Δϕ                            # eastern edge
             # the last line does not require a mod and in fact throws
             # an error if, as long as the offset from prime meridian
             # is at most Δϕ/2 (true for HEALPixGrids, offset 0 for
             # Gaussian and Clenshaw grids)
 
             # matrix indices for input grid that lie in output grid cell
-            a = findfirst(ϕ -> ϕ >= ϕ0, lon_in)
-            b = findlast( ϕ -> ϕ <  ϕ1, lon_in)
+            a = findlast( ϕ -> ϕ <  ϕ0, lon_in)     # western edge
+            b = findfirst(ϕ -> ϕ >= ϕ1, lon_in)     # eastern edge
+
+            # map around prime meridian if coordinates outside of range
+            a = isnothing(a) ? nlon_in : a
+            b = isnothing(b) ? 1 : b
 
             # in most cases we will have 1 <= a < b <=n, then loop i in a:b (b:a isn't looping)
             # however at the edges we have a < 1 or n < b the mod turns this into
             # 1 <= b < a <= n, for which we want to loop 1:b AND a:n
             ab, ba = b < a ? (1:b, a:nlon_in) : (a:b, b:a)
 
             sum_of_weights = 0
-            @inbounds for j′ in j0:j1
+            for j′ in j0:j1
 
                 for i in ab
                     w = coslat[j′]

diff --git a/src/physics/land_sea_mask.jl b/src/physics/land_sea_mask.jl
@@ -95,6 +95,8 @@ function initialize!(land_sea_mask::LandSeaMask, model::PrimitiveEquation)
     RingGrids.grid_cell_average!(land_sea_mask.mask, lsm_highres)
 
     # TODO this shoudln't be necessary, but at the moment grid_cell_average! can return values > 1
+    # lo, hi = extrema(land_sea_mask.mask)
+    # (lo < 0 || hi > 1) && @warn "Land-sea mask has values in [$lo, $hi], clamping to [0, 1]."
     clamp!(land_sea_mask.mask, 0, 1)
 end
 

diff --git a/test/grids/interpolation.jl b/test/grids/interpolation.jl
@@ -210,4 +210,23 @@ end
     RingGrids.interpolate!(C, A, interpolator)
 
     @test B == C
+end
+
+@testset "Grid cell average" begin
+    for Grid in (   FullGaussianGrid,
+                    FullClenshawGrid,
+                    OctahedralGaussianGrid,
+                    OctahedralClenshawGrid,
+                    OctaminimalGaussianGrid,
+                    HEALPixGrid,
+                    OctaHEALPixGrid)
+
+        for trunc in (31, 42, 63)
+
+            spectral_grid = SpectralGrid(; trunc, Grid)
+
+            land_sea_mask = LandSeaMask(spectral_grid)
+            initialize!(land_sea_mask, PrimitiveDryModel(spectral_grid))
+        end
+    end
 end