Fix slope azimuth of exactly east- or west-facing cells

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
 

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.
 

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()