From 09e13eef1a61759cf7bcb684f2629bb8c75b3e48 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Wed, 5 Aug 2026 14:13:04 -0400
Subject: [PATCH] face_areas fesom test, & warning if total too big

---
 test/grid/grid/test_areas.py |  7 +++++++
 uxarray/grid/grid.py         | 20 +++++++++++++++++---
 2 files changed, 24 insertions(+), 3 deletions(-)

diff --git a/test/grid/grid/test_areas.py b/test/grid/grid/test_areas.py
index 4e720ef04..f17681b34 100644
--- a/test/grid/grid/test_areas.py
+++ b/test/grid/grid/test_areas.py
@@ -115,3 +115,10 @@ def test_latlon_bounds_populate_bounds_MPAS(gridpath):
     """Test bounds population with MPAS grid."""
     uxgrid = ux.open_grid(gridpath("mpas", "QU", "oQU480.231010.nc"))
     bounds_xarray = uxgrid.bounds
+
+
+def test_face_areas_fesom(gridpath):
+    """Ensure correct total area for FESOM grid (~8.3780 sr). Regression test for #425."""
+    uxgrid = ux.open_grid(gridpath("ugrid", "fesom", "fesom.mesh.diag.nc"))
+    total_area = uxgrid.calculate_total_face_area()
+    nt.assert_almost_equal(total_area, 8.3780, decimal=4)
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 66ba087bd..27ef00864 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -1958,7 +1958,12 @@ def calculate_total_face_area(
         """Calculate the total surface area of all the faces in a mesh.
 
         Equivalent to ``self.compute_face_areas(...).sum()``; provided as a
-        convenience.
+        convenience. (Note: for HEALPix grids, when called with default arguments,
+        this method actually returns ``self.face_areas.sum()`` instead,
+        which respects HEALPix equal-area property.)
+
+        Additionally, raises a warning if the result is larger than
+        the total area of a sphere (4 * pi * self.sphere_radius**2).
 
         Parameters
         ----------
@@ -1986,9 +1991,9 @@ def calculate_total_face_area(
             and order == 4
             and not latitude_adjusted_area
         ):
-            return np.sum(self.face_areas.values)
+            result = np.sum(self.face_areas.values)
 
-        return np.sum(
+        result = np.sum(
             self.compute_face_areas(
                 quadrature_rule=quadrature_rule,
                 order=order,
@@ -1996,6 +2001,15 @@ def calculate_total_face_area(
             )
         )
 
+        RTOL = 1e-6  # 1e-7 had warnings in existing CI tests (as of 2026-08-05), 1e-6 did not.
+        if result > 4 * np.pi * self.sphere_radius**2 * (1 + RTOL):
+            warnings.warn(
+                f"Total face area (={result}) exceeds the surface area of the whole sphere "
+                f"(={4 * np.pi * self.sphere_radius**2}) (with sphere_radius={self.sphere_radius}).",
+            )
+
+        return result
+
     def compute_face_areas(
         self,
         quadrature_rule: str = "triangular",