Fix calc_zostoga: reference thickness, temperature-dependent alpha, optional temp_ref (#288)

* fix: update calc_zostoga to use reference thickness and add temp_ref handling

* fix: convert pot_temp and temp_ref from K to °C in calc_zostoga

The mapping declares pot_temp in K, but calc_zostoga was passing the
raw value directly into the Gill 1982 α(T) formula and the (T − T_ref)
anomaly, both of which require °C. With values of ~280–300 K the
thermal expansion coefficient and steric height were substantially wrong.

Changes:
- Convert pot_temp K→°C internally (pot_temp_c = pot_temp - 273.15).
- Change temp_ref unit contract from °C to K (matching pot_temp) and
  convert temp_ref K→°C internally (temp_ref_c = temp_ref - 273.15).
  This eliminates the asymmetry where a caller building temp_ref as a
  time-mean of pot_temp would naturally have it in K and would silently
  pass the wrong units.
- Update fallback default from 4.0 °C → 277.15 K (= 4 °C).
- Update docstring and warning message to reflect the new K contract.
- Update unit tests: pot_temp values +273.15 K, temp_ref=277.15 K.

* fix: add area_t to model_variables and calculation args in ACCESS-ESM1.6 mappings

Author: rbeucher

src/access_moppy/derivations/calc_ocean.py

@@ -1,5 +1,6 @@
 #!/usr/bin/env python
 import logging
+import warnings
 
 import xarray as xr
 
@@ -81,63 +82,99 @@ def calc_rsdoabsorb(sw_heat: xr.DataArray, swflux: xr.DataArray) -> xr.DataArray
     return rsdoabsorb
 
 
-def calc_zostoga(pot_temp, dzt, areacello, depth_coord="st_ocean"):
+def calc_zostoga(pot_temp, dzt_ref, areacello, temp_ref=None, depth_coord="st_ocean"):
     """Calculate Global Average Thermosteric Sea Level Change.
 
-    This function computes thermosteric sea level change by comparing
-    in-situ density to reference density at 4°C, integrating over depth,
-    and computing the global average.
+    Computes thermosteric sea level using a temperature-dependent thermal
+    expansion coefficient integrated against a reference layer-thickness field.
+    Using a fixed reference thickness (rather than the time-varying model dzt)
+    is required for Boussinesq models such as MOM5 to isolate the thermosteric
+    signal from the barotropic free-surface contribution.
 
-    Uses simplified density calculations suitable for thermosteric computations.
+    All operations are dask-lazy: no ``.compute()`` or ``.values`` calls are
+    made, so large datasets can be processed out-of-core.
 
     Parameters
     ----------
     pot_temp : xarray.DataArray
-        Potential temperature in degrees Celsius
+        Potential temperature in Kelvin (K), as provided by the model.
+        The function converts to degrees Celsius internally before applying
+        the Gill 1982 thermal expansion formula.
         Dimensions: (time, depth, lat, lon)
-    dzt : xarray.DataArray
-        Model level thickness with same dimensions as pot_temp
+    dzt_ref : xarray.DataArray
+        Reference (time-invariant) model level thickness.
+        Dimensions: (depth, lat, lon) or (depth,)
         Units: m
     areacello : xarray.DataArray
-        Ocean grid cell areas
+        Ocean grid cell areas.
         Dimensions: (lat, lon)
         Units: m²
+    temp_ref : float or xarray.DataArray or None, optional
+        Reference temperature in Kelvin (K), matching the units of ``pot_temp``.
+        If None (default), falls back to 277.15 K (= 4 °C) with a scientific
+        warning.  Using a scalar 277.15 K is a temporary approximation — it
+        computes an absolute steric height relative to 4 °C rather than a
+        temporal anomaly, which is not CMIP-compliant.  The result will carry
+        a large, physically meaningless baseline offset whose magnitude depends
+        on how far the mean ocean temperature departs from 4 °C.  For a
+        physically correct result, pass a 3-D reference-period mean temperature
+        field with dimensions (depth, lat, lon) in K (e.g. the piControl
+        year-1 mean of ``pot_temp``).
     depth_coord : str, optional
         Name of the depth coordinate, default 'st_ocean'
 
     Returns
     -------
     zostoga : xarray.DataArray
-        Global Average Thermosteric Sea Level Change
+        Global Average Thermosteric Sea Level Change.
         Dimensions: (time,)
         Units: m
 
     Notes
     -----
-    Uses simplified seawater density approximation:
-    - Reference salinity: 35 PSU
-    - Reference temperature: 4°C
-    - Linear thermal expansion coefficient: ~2e-4 /°C
+    Uses a temperature-dependent thermal expansion coefficient (Gill 1982):
+        α(T) ≈ 5.27×10⁻⁵ + 7.1×10⁻⁶·T − 4×10⁻⁸·T²  [°C⁻¹]
+    valid for S ≈ 35 PSU at surface pressure. This avoids the large errors
+    introduced by a constant α in cold deep and polar waters.
     """
+    if temp_ref is None:
+        warnings.warn(
+            "calc_zostoga: 'temp_ref' was not provided.  Falling back to a "
+            "scalar reference temperature of 277.15 K (= 4 °C).\n"
+            "Scientific implications:\n"
+            "  * zostoga will be an ABSOLUTE steric height relative to 4 °C, "
+            "not a temporal anomaly as required by CMIP.  The result will "
+            "carry a large, physically meaningless baseline offset whose "
+            "magnitude depends on how far the mean ocean temperature departs "
+            "from 4 °C.\n"
+            "  * Differences between time steps are still meaningful, but the "
+            "absolute values and any multi-model comparison will be incorrect.\n"
+            "To fix this, pass a 3-D reference-period mean temperature field "
+            "in K (e.g. the piControl year-1 mean of pot_temp) as 'temp_ref'.",
+            UserWarning,
+            stacklevel=2,
+        )
+        temp_ref = 277.15
 
-    # Simple approximation for seawater density temperature dependence
-    # ρ(T) ≈ ρ₀[1 - α(T - T₀)] where α ≈ 2e-4 /°C for seawater
-    rho_0 = 1025.0  # Reference density at 4°C, kg/m³
-    temp_ref = 4.0  # Reference temperature, °C
-    alpha = 2e-4  # Thermal expansion coefficient, /°C
+    # Convert both pot_temp and temp_ref from Kelvin to Celsius.
+    # All inputs are expected in K (matching the model/mapping convention);
+    # the Gill 1982 α(T) formula and the (T − T_ref) anomaly both require °C.
+    pot_temp_c = pot_temp - 273.15
+    temp_ref_c = temp_ref - 273.15
 
-    # Calculate density at in-situ temperature
-    rho_insitu = rho_0 * (1.0 - alpha * (pot_temp - temp_ref))
+    # Temperature-dependent thermal expansion coefficient (Gill 1982, simplified)
+    # α(T) ≈ 5.27e-5 + 7.1e-6·T − 4e-8·T²  [°C⁻¹], valid for S≈35 PSU
+    alpha = 5.27e-5 + 7.1e-6 * pot_temp_c - 4e-8 * pot_temp_c**2
 
-    # Calculate reference density at 4°C
-    rho_ref = rho_0  # At reference temperature
+    # Thermosteric height contribution: α(T) × (T − T_ref) × dz_ref
+    # dzt_ref is time-invariant so it does not carry the free-surface
+    # barotropic signal that is present in the model's time-varying dzt.
+    thermo_height = alpha * (pot_temp_c - temp_ref_c) * dzt_ref
 
-    # Calculate thermosteric height change and integrate over depth
-    # (1 - ρ/ρ₄) gives fractional density difference
-    thermo_height = (1.0 - rho_insitu / rho_ref) * dzt
+    # Integrate over depth (lazy with dask)
     integrated_height = thermo_height.sum(dim=depth_coord, skipna=True)
 
-    # Calculate area-weighted global average
+    # Area-weighted global average (lazy with dask)
     zostoga = integrated_height.weighted(areacello).mean(dim=["yt_ocean", "xt_ocean"])
 
     return zostoga

src/access_moppy/mappings/ACCESS-ESM1.6_mappings.json

@@ -5570,20 +5570,23 @@
             "positive": null,
             "model_variables": [
                 "pot_temp",
-                "dzt"
+                "dzt",
+                "area_t"
             ],
             "calculation": {
                 "type": "formula",
                 "operation": "calc_zostoga",
                 "args": [
                     "pot_temp",
-                    "dzt"
+                    "dzt",
+                    "area_t"
                 ]
             },
             "CF standard Name": "",
             "model_variable_units": {
                 "pot_temp": "K",
-                "dzt": "m"
+                "dzt": "m",
+                "area_t": "m2"
             }
         }
     },

tests/unit/test_derivations_calc_ocean.py

@@ -1,5 +1,6 @@
 """Tests for access_moppy.derivations.calc_ocean."""
 
+import dask.array as da
 import numpy as np
 import pytest
 import xarray as xr
@@ -165,55 +166,94 @@ class TestCalcZostoga:
     @pytest.mark.unit
     def test_returns_dataarray(self):
         pot_temp = xr.DataArray(
-            np.ones((NT, NZ, NY, NX)) * 10.0,
+            np.ones((NT, NZ, NY, NX)) * 283.15,  # 10 °C in K
             dims=["time", "st_ocean", "yt_ocean", "xt_ocean"],
         )
-        dzt = xr.DataArray(
-            np.ones((NT, NZ, NY, NX)) * 100.0,
-            dims=["time", "st_ocean", "yt_ocean", "xt_ocean"],
+        dzt_ref = xr.DataArray(
+            np.ones((NZ, NY, NX)) * 100.0,
+            dims=["st_ocean", "yt_ocean", "xt_ocean"],
         )
         areacello = xr.DataArray(
             np.ones((NY, NX)) * 1e10,
             dims=["yt_ocean", "xt_ocean"],
         )
-        result = calc_zostoga(pot_temp, dzt, areacello)
+        result = calc_zostoga(pot_temp, dzt_ref, areacello)
         assert isinstance(result, xr.DataArray)
 
     @pytest.mark.unit
     def test_time_dim_preserved(self):
         pot_temp = xr.DataArray(
-            np.ones((NT, NZ, NY, NX)) * 10.0,
+            np.ones((NT, NZ, NY, NX)) * 283.15,  # 10 °C in K
             dims=["time", "st_ocean", "yt_ocean", "xt_ocean"],
         )
-        dzt = xr.DataArray(
-            np.ones((NT, NZ, NY, NX)),
-            dims=["time", "st_ocean", "yt_ocean", "xt_ocean"],
+        dzt_ref = xr.DataArray(
+            np.ones((NZ, NY, NX)),
+            dims=["st_ocean", "yt_ocean", "xt_ocean"],
         )
         areacello = xr.DataArray(
             np.ones((NY, NX)),
             dims=["yt_ocean", "xt_ocean"],
         )
-        result = calc_zostoga(pot_temp, dzt, areacello)
+        result = calc_zostoga(pot_temp, dzt_ref, areacello)
         assert "time" in result.dims
 
     @pytest.mark.unit
     def test_zero_thermosteric_at_reference_temp(self):
-        """At the reference temperature (4°C), thermosteric change should be ~0."""
+        """When pot_temp == temp_ref, thermosteric change should be ~0."""
         pot_temp = xr.DataArray(
-            np.ones((NT, NZ, NY, NX)) * 4.0,  # exactly at reference
+            np.ones((NT, NZ, NY, NX)) * 277.15,  # 4 °C in K
             dims=["time", "st_ocean", "yt_ocean", "xt_ocean"],
         )
-        dzt = xr.DataArray(
-            np.ones((NT, NZ, NY, NX)) * 10.0,
-            dims=["time", "st_ocean", "yt_ocean", "xt_ocean"],
+        dzt_ref = xr.DataArray(
+            np.ones((NZ, NY, NX)) * 10.0,
+            dims=["st_ocean", "yt_ocean", "xt_ocean"],
         )
         areacello = xr.DataArray(
             np.ones((NY, NX)),
             dims=["yt_ocean", "xt_ocean"],
         )
-        result = calc_zostoga(pot_temp, dzt, areacello)
+        # Pass temp_ref in K (same units as pot_temp) to avoid the fallback warning
+        result = calc_zostoga(pot_temp, dzt_ref, areacello, temp_ref=277.15)
         np.testing.assert_allclose(result.values, 0.0, atol=1e-10)
 
+    @pytest.mark.unit
+    def test_warns_when_temp_ref_not_provided(self):
+        """A UserWarning must be raised when temp_ref is omitted."""
+        pot_temp = xr.DataArray(
+            np.ones((NT, NZ, NY, NX)) * 283.15,  # 10 °C in K
+            dims=["time", "st_ocean", "yt_ocean", "xt_ocean"],
+        )
+        dzt_ref = xr.DataArray(
+            np.ones((NZ, NY, NX)) * 10.0,
+            dims=["st_ocean", "yt_ocean", "xt_ocean"],
+        )
+        areacello = xr.DataArray(
+            np.ones((NY, NX)),
+            dims=["yt_ocean", "xt_ocean"],
+        )
+        with pytest.warns(UserWarning, match="temp_ref.*not provided"):
+            calc_zostoga(pot_temp, dzt_ref, areacello)
+
+    @pytest.mark.unit
+    def test_dask_lazy(self):
+        """Result should remain dask-backed (lazy) when inputs are dask arrays."""
+        pot_temp = xr.DataArray(
+            da.from_array(
+                np.ones((NT, NZ, NY, NX)) * 283.15, chunks=(1, NZ, NY, NX)
+            ),  # 10 °C in K
+            dims=["time", "st_ocean", "yt_ocean", "xt_ocean"],
+        )
+        dzt_ref = xr.DataArray(
+            da.from_array(np.ones((NZ, NY, NX)) * 10.0, chunks=(NZ, NY, NX)),
+            dims=["st_ocean", "yt_ocean", "xt_ocean"],
+        )
+        areacello = xr.DataArray(
+            da.from_array(np.ones((NY, NX)), chunks=(NY, NX)),
+            dims=["yt_ocean", "xt_ocean"],
+        )
+        result = calc_zostoga(pot_temp, dzt_ref, areacello)
+        assert isinstance(result.data, da.Array)
+
 
 # ---------------------------------------------------------------------------
 # calc_overturning_streamfunction