From e43d3a7e9c4de9eea3a3e76aa526586355f544de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Milan=20Kl=C3=B6wer?= Date: Tue, 11 Feb 2025 23:19:54 +0100 Subject: [PATCH] Grid cell average bug fixes (#673) * grid cell average bug fixes * changelog updated --- CHANGELOG.md | 1 + src/RingGrids/interpolation.jl | 12 ++++++++---- src/physics/land_sea_mask.jl | 2 ++ test/grids/interpolation.jl | 19 +++++++++++++++++++ 4 files changed, 30 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6c61369dc..98449851d 100644 --- 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 index 9e352a96b..d4a3ec3a9 100644 --- a/src/RingGrids/interpolation.jl +++ b/src/RingGrids/interpolation.jl @@ -533,7 +533,7 @@ 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 @@ -541,8 +541,12 @@ function grid_cell_average!( # 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 @@ -550,7 +554,7 @@ function grid_cell_average!( 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 index b5b470108..97e4c870d 100644 --- 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 index cda75cca6..eccf49533 100644 --- 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 \ No newline at end of file