diff --git a/CHANGELOG.md b/CHANGELOG.md index df0a7ac4..43b0a3e1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,7 @@ SPDX-License-Identifier: CC-BY-4.0 * In a run on several MPI ranks, only rank 0 logs to the console at `--log-level-console`; the other ranks only log warnings and errors there, instead of repeating every message. Log files (`--log-file`, one per rank) are unchanged. The current model time is logged instead of printed. https://github.com/EBFMorg/EBFM/pull/204 * A MATLAB grid (`--matlab-mesh`) can now run on several MPI ranks, e.g. `mpirun -np 4 ebfm --matlab-mesh examples/dem_and_mask.mat --no-shading`: each rank computes a contiguous block of the columns, and rank 0 writes the output file, which is identical to that of a run on a single rank. Shading is not supported on several ranks and stops the run with an error. https://github.com/EBFMorg/EBFM/pull/202 * A usage error, e.g. a rejected command line option, now ends the whole MPI job with its exit status when `MPI.COMM_WORLD` has more than one process. EBFM previously exited on its own, and the other processes, e.g. a coupled model in an MPMD run, waited for it forever. A run on a single process still exits without `MPI_ABORT`. https://github.com/EBFMorg/EBFM/pull/206 +* Fixed the slope azimuth of grid cells on MATLAB and NetCDF meshes that face exactly east or west (north-south gradient exactly zero): they were treated as facing the opposite direction, so the incoming shortwave radiation of an east-facing cell peaked in the afternoon instead of the morning. Results change in these cells. https://github.com/EBFMorg/EBFM/pull/194 # v0.8.0 diff --git a/src/ebfm/core/INIT.py b/src/ebfm/core/INIT.py index cd868639..203e5cf2 100644 --- a/src/ebfm/core/INIT.py +++ b/src/ebfm/core/INIT.py @@ -386,19 +386,7 @@ def init_grid(grid: GridDict, io, config: GridConfig, comm): # Calculate slope_beta and slope_gamma (defining the tilt and orientation of a sloping surface) grid["slope_beta"] = np.arctan(grid["slope"]) - grid["slope_gamma"] = np.zeros_like(grid["slope"]) - grid["slope_gamma"][grid["slope_y"] > 0] = np.arctan( - -grid["slope_x"][grid["slope_y"] > 0] / grid["slope_y"][grid["slope_y"] > 0] - ) - cond = np.logical_and(grid["slope_y"] < 0, grid["slope_x"] > 0) - grid["slope_gamma"][cond] = -np.pi + np.arctan(-grid["slope_x"][cond] / grid["slope_y"][cond]) - cond = np.logical_and(grid["slope_y"] < 0, grid["slope_x"] < 0) - grid["slope_gamma"][cond] = np.pi + np.arctan(-grid["slope_x"][cond] / grid["slope_y"][cond]) - grid["slope_gamma"][(grid["slope_x"] == 0) & (grid["slope_y"] < 0)] = np.pi - grid["slope_gamma"][(grid["slope_x"] == 0) & (grid["slope_y"] == 0)] = 0 - grid["slope_gamma"][(grid["slope_x"] > 0) & (grid["slope_y"] == 0)] = np.pi / 2 - grid["slope_gamma"][(grid["slope_x"] < 0) & (grid["slope_y"] == 0)] = -np.pi / 2 - grid["slope_gamma"] = -grid["slope_gamma"] + grid["slope_gamma"] = slope_azimuth(grid["slope_x"], grid["slope_y"]) # ----------------------------------------------------------------------------------------------------- # Pre-compute maximum grid elevation angle for various azimuth angles (needed for shading calculation) @@ -422,6 +410,22 @@ def init_grid(grid: GridDict, io, config: GridConfig, comm): return grid +def slope_azimuth(slope_x: np.ndarray, slope_y: np.ndarray) -> np.ndarray: + """ + Azimuth of the direction a sloping surface faces, as used for the incidence angle in LOOP_EBM_SWin. + + 0 is south-facing, -pi/2 east-facing, pi/2 west-facing and +-pi north-facing. Flat cells get 0. + + Parameters: + slope_x (np.ndarray): Surface gradient towards east (dz/dx). + slope_y (np.ndarray): Surface gradient towards north (dz/dy). + """ + gamma = np.arctan2(slope_x, slope_y) + # arctan2(0, -0.0) is pi + gamma[(slope_x == 0) & (slope_y == 0)] = 0.0 + return gamma + + def _precompute_shading_matlab(grid: GridDict, mask_2D) -> None: """Pre-compute the shading look-up table of a MATLAB grid. diff --git a/tests/core/test_init.py b/tests/core/test_init.py index 9381f61c..6647799e 100644 --- a/tests/core/test_init.py +++ b/tests/core/test_init.py @@ -475,5 +475,47 @@ def test_valid_orography_of_inactive_cell_is_kept(self): self.assertEqual(data["z"][1, 1], 500.0) +class TestSlopeAzimuth(unittest.TestCase): + """ + slope_gamma for slopes along the grid axes, where one gradient component is exactly zero. + + An exactly east- or west-facing cell used to get the azimuth of the opposite direction. + """ + + def test_every_sign_combination(self): + # (slope_x, slope_y): expected azimuth; a surface faces downhill + expected = { + (-1.0, -1.0): -3 * np.pi / 4, # north-east + (-1.0, 0.0): -np.pi / 2, # east + (-1.0, 1.0): -np.pi / 4, # south-east + (0.0, -1.0): np.pi, # north + (0.0, 0.0): 0.0, # flat + (0.0, 1.0): 0.0, # south + (1.0, -1.0): 3 * np.pi / 4, # north-west + (1.0, 0.0): np.pi / 2, # west + (1.0, 1.0): np.pi / 4, # south-west + } + for (slope_x, slope_y), gamma in expected.items(): + with self.subTest(slope_x=slope_x, slope_y=slope_y): + result = INIT.slope_azimuth(np.array([slope_x]), np.array([slope_y])) + np.testing.assert_allclose(result, gamma, rtol=0, atol=1e-15) + + def test_matches_arctan2(self): + slope_x = np.array([0.3, -0.3, 0.3, -0.3, 0.3, -0.3, 0.0, 0.0, -0.0, 0.3, -0.3]) + slope_y = np.array([0.2, 0.2, -0.2, -0.2, 0.0, 0.0, 0.2, -0.2, -0.2, -0.0, -0.0]) + np.testing.assert_array_equal(INIT.slope_azimuth(slope_x, slope_y), np.arctan2(slope_x, slope_y)) + + def test_east_and_west_facing_match_their_limits(self): + for slope_x in (0.3, -0.3): + with self.subTest(slope_x=slope_x): + eps = 1e-12 + gamma = INIT.slope_azimuth(np.array([slope_x] * 4), np.array([eps, 0.0, -0.0, -eps])) + np.testing.assert_allclose(gamma, np.sign(slope_x) * np.pi / 2, rtol=0, atol=1e-10) + + def test_flat_is_zero(self): + gamma = INIT.slope_azimuth(np.array([0.0, -0.0, 0.0, -0.0]), np.array([0.0, 0.0, -0.0, -0.0])) + np.testing.assert_array_equal(gamma, 0.0) + + if __name__ == "__main__": unittest.main()