Repository: ImperialCollegeLondon/pyrealm

Issue:
Update FaparLimitation to allow alternative methods
With the reference code for the Zhu methods, we can now build out `FaparLimitation` to support alternative approaches. The main differences here:

* The current code (Cai method) uses an `f0` cost that varies with aridity and then just takes the minimum of the energy and water limited fAPAR.
* The Zhu method has a fixed `f0` cost and then has a more complex function that compute maximum fapar from energy and water limited fAPAR. It doesn't use aridity data.

Soooo... I think the solution is probably to

* Update FaparLimitation to handle optional variables like `PModelEnvironment` currently does.
* Follow the practice elsewhere in `pyrealm` of breaking out method specific differences into subclasses of a new `FaparMethodABC`. The class should:
    * Check any required variables are provided in `__init__`
    * Provide a `set_f0_and_z` method that is called on `__init__` to populate those parameters
    * Provide a method that returns the maximum fAPAR.

Discussion:

Fix (patch):
diff --git a/pyrealm/constants/phenology_const.py b/pyrealm/constants/phenology_const.py
index db227c3a..d53eef22 100644
--- a/pyrealm/constants/phenology_const.py
+++ b/pyrealm/constants/phenology_const.py
@@ -61,3 +61,41 @@ def calculate_f0(self, aridity_index: NDArray[np.floating]) -> NDArray[np.floati
 
         a, b, c = self.f0_coefficients
         return a * np.exp(-b * np.log(aridity_index / c) ** 2)
+
+
+@dataclass(frozen=True)
+class PhenologyConstNew(ConstantsClass):
+    r"""Model constants for the phenology module class.
+
+    This data class defines constants used in the phenology module.
+    """
+
+    cai_z: float = 12.227
+    r"""The constant :math:`z` (mol C m^{-2} year^{-1}) accounts for the carbon costs of
+    building and maintaining leaves and the total below-ground allocation required to
+    support the nutrient demand of those leaves. The default is an empirical estimate
+    from global data."""
+
+    k: float = 0.5
+    """The canopy light extinction coefficient (unitless)"""
+
+    cai_f0_coefficients: tuple[float, float, float] = (0.65, 0.604169, 1.9)
+    r"""Coefficients :math:`a,b,c` to calculate :math:`f_0` from the local aridity index
+    (AI), where :math:`a` is the maximum value, :math:`b` is the slope of the
+    relationship with AI and :math:`c` is the AI value at which the maximum value
+    :math:`a` is reached.
+    """
+
+    cai_sigma: float = 0.771
+    r"""The :math:`\sigma` parameter captures the ability of plants to maintain maximum
+    :term:`LAI` throughout the growing season. If a canopy with maximum LAI grows
+    instantaneously and is present until the end of the growing season, then the
+    potential LAI profile forms a square wave (:math:`\sigma = 1`). However real plants
+    require time to deploy a canopy at the start of the growing season, through bud
+    formation, budburst and leaf growth. Time is also needed at the end of the season
+    for processes such as nutrient resorption and leaf sensescence. These processes
+    round off that square wave and this is captured by values of :math:`\sigma < 1`."""
+
+    zhu_f0: float = 0.5
+    zhu_z: float = 17
+    zhu_budyko: float = 4
diff --git a/pyrealm/phenology/fapar_limitation_new.py b/pyrealm/phenology/fapar_limitation_new.py
new file mode 100644
index 00000000..4709639e
--- /dev/null
+++ b/pyrealm/phenology/fapar_limitation_new.py
@@ -0,0 +1,624 @@
+"""Class to compute the fAPAR_max and annual peak Leaf Area Index (LAI).
+
+The :class:`FaparLimitation` class and the :meth:`FaparLimitation.from_pmodel` are
+designed to work with inputs that can have multiple dimensions. The first axis is
+_always_ assumed to represent a time series of annual observations. If the inputs
+are one dimensional, then this is a time series for a single site; if they are three
+dimensional then these are observations for a grid of sites. Usually all array inputs
+will have the same shape but note the following instances where you might need to
+take care with array broadcasting.
+
+* Growing season length might well be constant across sites. If so - for example
+    with 3D data - the input would need shape `(N, 1, 1)` to broadcast N years of
+    data over the array of sites.
+* The climatological aridity index is very likely to be constant through time.
+    If so, the array should have a singleton first dimension to broadcast an
+    observation per site across observations. For example, with `(10, 3, 3)` data
+    (10 years for a 3x3 grid of sites), the aridity index could be provided as
+    `(1,3,3)` to broadcast the aridity index across each year. It could also use a
+    single scalar value to use the same aridity index for all sites.
+"""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+
+import numpy as np
+from numpy.typing import NDArray
+
+from pyrealm.constants.phenology_const import PhenologyConstNew
+from pyrealm.core.experimental import warn_experimental
+from pyrealm.core.time_series import AnnualValueCalculator
+from pyrealm.core.utilities import (
+    check_input_shapes,
+    summarize_attrs,
+)
+from pyrealm.pmodel.pmodel import PModel, PModelABC, SubdailyPModel
+
+FAPAR_LIMITATION_METHOD_CLASS_REGISTRY: dict[str, type[FaparLimitationMethodABC]] = {}
+"""A registry for optimal chi calculation classes.
+
+Different implementations of the calculation of optimal chi must all be subclasses of
+:class:`~pyrealm.pmodel.optimal_chi.OptimalChiABC` abstract base class.
+This dictionary is used as a registry for defined subclasses and a method name
+is used to retrieve a particular implementation from this registry. For example:
+
+.. code:: python
+
+    prentice14_opt_chi = OPTIMAL_CHI_CLASS_REGISTRY['prentice14']
+"""
+
+
+class FaparLimitationMethodABC(ABC):
+    r"""Abstract base class for methods computing maximum annual fAPAR and LAI."""
+
+    __experimental__ = True
+
+    method: str
+    """A short method name used to identify the class in
+    :data:`FAPAR_LIMITATION_METHOD_CLASS_REGISTRY`.
+    """
+
+    attrs: tuple[tuple[str, str], ...]
+    """A tuple of attributes to be reported in the FaparLimitation.summarize() output
+    for the method."""
+
+    requires: tuple[str, ...]
+    """A tuple of any additional variables required to create a FaparLimitation
+    instance using the method."""
+
+    def __init__(
+        self, fapar_limitation: FaparLimitationNew, phenology_const: PhenologyConstNew
+    ):
+        """Initialise the method instance.
+
+        Subclass should provide a method specific instance that calls
+        `super().__init__(...)`, carries out method specific validation and then runs
+        ``self.set_z_and_f0``.
+        """
+        self.fapar_limitation = fapar_limitation
+        self.phenology_const = phenology_const
+
+        self.f0: NDArray[np.floating]
+        self.z: NDArray[np.floating]
+
+        # Check for required variables
+        self._check_required_variables()
+
+    def _check_required_variables(self) -> None:
+        """Check required variables.
+
+        Checks that any required variables for the method have been passed to the
+        FaparLimitation constructor.
+        """
+        for var in self.requires:
+            if not hasattr(self.fapar_limitation, var):
+                raise ValueError(
+                    f"Values for {var} are required to use the {self.method} method "
+                    "with FaparLimitation."
+                )
+
+    @abstractmethod
+    def set_z_and_f0(self) -> None:
+        """Sets the f0 and z values for the method."""
+        pass
+
+    @abstractmethod
+    def calculate_maximum_fapar(
+        self, energy_limited_fapar: NDArray, water_limited_fapar: NDArray
+    ) -> NDArray[np.floating]:
+        """Calculate the maximum fAPAR.
+
+        Provides the method specific calculation of maximum fAPAR from the energy and
+        water limited maximum values and should return the calculated maximum fAPAR.
+
+        Args:
+            energy_limited_fapar: The maximum fAPAR given energy limitation.
+            water_limited_fapar: The maximum fAPAR given water limitation.
+        """
+        pass
+
+    @classmethod
+    def __init_subclass__(cls, method: str) -> None:
+        """Initialise a subclass deriving from this ABC."""
+
+        cls.method = method
+        FAPAR_LIMITATION_METHOD_CLASS_REGISTRY[cls.method] = cls
+
+
+class FaparLimitationCai(FaparLimitationMethodABC, method="cai"):
+    r"""Compute maximum annual fAPAR and LAI using the method of :cite:`cai:2025a`.
+
+    TBD - what is specific here versus general in FaparLimitation.
+
+    """
+
+    __experimental__ = True
+
+    attrs = (
+        ("lai_max", "-"),
+        ("fapar_max", "-"),
+        ("lai_to_gpp_ratio_m", "-"),
+    )
+
+    requires = ("aridity_index",)
+
+    def __init__(
+        self, fapar_limitation: FaparLimitationNew, phenology_const: PhenologyConstNew
+    ):
+        """Initialise a FaparLimitationMethod instance using the Cai approach."""
+
+        # Run the superclass init method.
+        super().__init__(
+            fapar_limitation=fapar_limitation, phenology_const=phenology_const
+        )
+
+        # This is only set as a side effect of the calculate_maximum_fapar method being
+        # called, which is a little bit hacky, but at the moment just preserving the
+        # attribute _somewhere_. Will see how this class evolves.
+        self.energy_limited: NDArray[np.bool_]
+        """Boolean array showing if annual :math:`fAPAR_{max}` is water or energy
+        limited."""
+
+        # Make sure the aridity index is not zero
+        self.aridity_index = getattr(self.fapar_limitation, "aridity_index")
+        if np.any(self.aridity_index <= 0):
+            raise ValueError("The aridity index has to be positive.")
+
+        # Set z and f0
+        self.set_z_and_f0()
+
+    def set_z_and_f0(self) -> None:
+        r"""Set the :math:`z` and :math:`f_0` parameters.
+
+        The value :math:`f_0` is the ratio of annual total transpiration of annual total
+        precipitation. It is calculated from site specific estimates of the
+        climatological aridity index, calculated as the long term (typically 20 years)
+        total PET over total precipitation (:math:`AI`, unitless) as:
+
+        .. math::
+
+                f_0 = a \exp{\left(-b \left(\frac{AI}{c}\right)^2\right)}
+
+        where :math:`a,b,c` are defined in the
+        :attr:`~pyrealm.constants.phenology_const.PhenologyConstNew.cai_f0_coefficients`
+        attribute.
+        """
+
+        a, b, c = self.phenology_const.cai_f0_coefficients
+
+        self.f0 = a * np.exp(-b * np.log(self.aridity_index / c) ** 2)
+        self.z = np.array([self.phenology_const.cai_z])
+
+    def calculate_maximum_fapar(
+        self, energy_limited_fapar: NDArray, water_limited_fapar: NDArray
+    ) -> NDArray:
+        """Calculate the maximum fAPAR.
+
+        The Cai method uses the simple minimum of the water and energy limited fAPAR
+        values as the maximum possible fAPAR.
+
+        Args:
+            energy_limited_fapar: The maximum fAPAR given energy limitation.
+            water_limited_fapar: The maximum fAPAR given water limitation.
+        """
+
+        # Calculate fAPAR max and record whether the location is energy or water limited
+        fapar_max = np.minimum(water_limited_fapar, energy_limited_fapar)
+        self.energy_limited = energy_limited_fapar < water_limited_fapar
+
+        # self.lai_to_gpp_ratio_m = (
+        #     self.phenology_const.cai_sigma
+        #     * self.annual_growing_season_length
+        #     * self.lai_max
+        # ) / (self.annual_total_potential_gpp * self.fapar_max)
+        # """The steady state ratio of leaf area index to potential GPP (:math:`m`)"""
+
+        return fapar_max
+
+
+class FaparLimitationZhu(FaparLimitationMethodABC, method="zhu"):
+    r"""Compute maximum annual fAPAR and LAI using the method of Zhu .
+
+    TBD.
+
+    """
+
+    __experimental__ = True
+
+    attrs = (
+        ("lai_max", "-"),
+        ("fapar_max", "-"),
+    )
+
+    requires = tuple()
+
+    def __init__(
+        self, fapar_limitation: FaparLimitationNew, phenology_const: PhenologyConstNew
+    ):
+        """Initialise a FaparLimitationMethod instance using the Zhu approach."""
+
+        # Run the superclass init method.
+        super().__init__(
+            fapar_limitation=fapar_limitation, phenology_const=phenology_const
+        )
+
+        # Set z and f0
+        self.set_z_and_f0()
+
+    def set_z_and_f0(self) -> None:
+        r"""Set the :math:`z` and :math:`f_0` parameters.
+
+        This method has fixed values for :math:`z` (see
+        :attr:`~pyrealm.constants.phenology_const.PhenologyConstNew.zhu_f0`) and
+        :math:`f_0` (see
+        :attr:`~pyrealm.constants.phenology_const.PhenologyConstNew.zhu_z`)
+        """
+
+        self.f0 = np.array([self.phenology_const.zhu_f0])
+        self.z = np.array([self.phenology_const.zhu_z])
+
+    def calculate_maximum_fapar(
+        self, energy_limited_fapar: NDArray, water_limited_fapar: NDArray
+    ) -> NDArray:
+        """Calculate the maximum fAPAR.
+
+        The Zhu method calculates the maximum fAPAR as a function of the energy and
+        water limited fAPAR:
+
+        TODO: equation.
+
+        Args:
+            energy_limited_fapar: The maximum fAPAR given energy limitation.
+            water_limited_fapar: The maximum fAPAR given water limitation.
+        """
+        # Calculate the ratio of energy limited to water limited fAPAR, safeguarding
+        # against divide by zero
+        cw_ratio = energy_limited_fapar / (
+            np.clip(water_limited_fapar, a_min=np.finfo(float).eps, a_max=None)
+        )
+
+        # Calculate the maximum fapar.
+        fapar_max = (
+            (1 + cw_ratio)
+            - (1 + cw_ratio**self.phenology_const.zhu_budyko)
+            ** (1 / self.phenology_const.zhu_budyko)
+        ) * water_limited_fapar
+
+        return fapar_max
+
+
+class FaparLimitationNew:
+    r"""Compute maximum annual fAPAR and LAI.
+
+    This class calculates maximum annual fAPAR, which can be limited either by the
+    availability of light energy ($f_{APAR_{c}}$) or by the availability of water
+    ($f_{APAR_{w}}$). The equations for these two variables, following :cite:`cai2025a:
+    are:
+
+    .. math::
+        :nowrap:
+
+        \[    
+          \begin{align*}
+            f_{APAR_{c}} &= 1 - \frac{z}{k A_0}\\
+            f_{APAR_{w}} &= \left(\frac{ c_a \left( 1 - \chi \right)}{ 1.6 D }\right)
+                            \left(\frac{ f_0 P }{ A_0 }\right) \\
+          \end{align*}
+        \]
+
+    The maximum fAPAR is then calculated as a function of :math:`f_{APAR_{c}}` and
+    :math:`f_{APAR_{w}}`. In these equations:
+
+    * :math:`z` accounts for the growth and maintenance costs of leaves.
+    * :math:`f_0` accounts for water limitation on annual assimulation and is is the
+      ratio of annual total transpiration to annual total precipitation.
+
+    The other variables are the required arguments to the class defined below. 
+    
+    There are different approaches to estimating :math:`z` and :math:`f_0` and to
+    calculating the maximum fAPAR from the two inputs. For details see:
+
+    * ``method=cai``; :class:`FaparLimitationMethodCai`
+    * ``method=zhu``; :class:`FaparLimitationMethodZhu`
+    
+    The maximum annual LAI can then be calculated using Beer's law as:
+
+    .. math::
+
+        L_{max} = - ( 1 / k ) \ln {1 -f_{APAR_{max}}}
+
+    The most common source of most of the variables needed to calculate maximum fAPAR is
+    a P Model, and the
+    :meth:`~pyrealm.phenology.fapar_limitation.FaparLimitation.from_pmodel` method can
+    be used to estimate maximum fAPAR directly from a fitted P Model.
+
+    Args:
+        annual_total_potential_gpp: The annual sum of potential GPP (:math:`A_0,
+            \text{mol C m}^{-2} \text{year}^{-1}`)
+        annual_mean_ca: The ambient CO2 partial pressure during the growing season
+            (:math:`c_a`, Pa)
+        annual_mean_chi: The annual mean ratio of ambient to leaf CO2 partial during the
+            growing season (:math:`\chi`, Pa)
+        annual_mean_vpd: The annual mean vapour pressure deficit during the growing
+            season (:math:`D`, Pa)
+        annual_total_precip: The annual total precipitation (:math:`P, \text{mol m}^{-2}
+            \text{year}^{-1}`)
+        annual_growing_season_length: The length of the growing season in days for each
+            year (:math:`G`, days)
+        years: An array of year datetimes for the observations.
+        method: The method to be applied when calculating maximum fAPAR, defaulting to
+        ``cai``. 
+        phenology_const: An instance of
+            :class:`~pyrealm.constants.phenology_const.PhenologyConstNew`
+        **kwargs: Any additional variables required by specific method choices.
+    """
+
+    __experimental__ = True
+
+    def __init__(
+        self,
+        annual_total_potential_gpp: NDArray[np.floating],
+        annual_mean_ca: NDArray[np.floating],
+        annual_mean_chi: NDArray[np.floating],
+        annual_mean_vpd: NDArray[np.floating],
+        annual_total_precip: NDArray[np.floating],
+        annual_growing_season_length: NDArray[np.floating],
+        years: NDArray[np.datetime64],
+        method: str = "cai",
+        phenology_const: PhenologyConstNew = PhenologyConstNew(),
+        **kwargs: NDArray[np.floating],
+    ) -> None:
+        # Experimental class
+        warn_experimental("FaparLimitation")
+
+        # Validate the input shapes.
+        self.shape: tuple[int, ...] = check_input_shapes(
+            annual_total_potential_gpp,
+            annual_mean_ca,
+            annual_mean_chi,
+            annual_mean_vpd,
+            annual_total_precip,
+            annual_growing_season_length,
+            *kwargs.values(),
+        )
+
+        # Check the years values - must be datetime64[Y] and be one dimensional,
+        # matching the first dimension of the other inputs
+        # TODO - this is a bit stringent, but is more robust
+        if not years.dtype == "<M8[Y]":
+            raise ValueError("The years argument must provide np.datetime64[Y] values")
+
+        if years.shape != (self.shape[0],):
+            raise ValueError(
+                "The years argument must be one dimensional and match the length "
+                "of the first axis of the other arguments"
+            )
+
+        self.years = years
+        r"""The year of each observation."""
+        self.annual_total_potential_gpp = annual_total_potential_gpp
+        r"""The annual sum of potential GPP 
+        (:math:`A_0, \text{mol C m}^{-2} \text{year}^{-1}`)"""
+        self.annual_mean_ca = annual_mean_ca
+        r"""Ambient CO2 partial pressure during the growing season (:math:`c_a`, Pa)"""
+        self.annual_mean_chi = annual_mean_chi
+        r"""Annual mean ratio of ambient to leaf CO2 partial during the 
+        growing season (:math:`\chi`, Pa)"""
+        self.annual_mean_vpd = annual_mean_vpd
+        r"""Annual mean vapour pressure deficit during the growing season (:math:`D`,
+        Pa)"""
+        self.annual_total_precip = annual_total_precip
+        r"""Annual total precipitation
+        (:math:`P, \text{mol m}^{-2} \text{year}^{-1}`)"""
+        self.annual_growing_season_length = annual_growing_season_length
+        r"""Annual growing season length (:math:`G`, days)"""
+
+        # Additional variables -  add them to the instance
+        for var_name, var_values in kwargs.items():
+            setattr(self, var_name, var_values)
+
+        self._additional_vars: tuple[str, ...] = tuple(kwargs.keys())
+        """A tuple containing the attribute names of additional variables passed to the
+        FaparLimitation instance."""
+
+        # Constants used for phenology computations
+        self.phenology_const: PhenologyConstNew = phenology_const
+
+        # Set up the fAPAR limitation method to be used.
+        self.method: str = method
+
+        if method not in FAPAR_LIMITATION_METHOD_CLASS_REGISTRY:
+            raise ValueError(f"Unknown FaparLimitation method: {method}")
+
+        # Get an instance of the method class from the registry
+        self.limitation_method: FaparLimitationMethodABC = (
+            FAPAR_LIMITATION_METHOD_CLASS_REGISTRY[method](
+                fapar_limitation=self, phenology_const=self.phenology_const
+            )
+        )
+
+        # Calculate the energy and water limited terms.
+        energy_limited_fapar = 1.0 - self.limitation_method.z / (
+            self.phenology_const.k * annual_total_potential_gpp
+        )
+        water_limited_fapar = (
+            self.limitation_method.f0
+            * annual_total_precip
+            * annual_mean_ca
+            * (1 - annual_mean_chi)
+            / (1.6 * annual_mean_vpd * annual_total_potential_gpp)
+        )
+
+        # Calculate the maximum fapar using the limitation method
+        self.fapar_max: NDArray[np.floating] = (
+            self.limitation_method.calculate_maximum_fapar(
+                energy_limited_fapar=energy_limited_fapar,
+                water_limited_fapar=water_limited_fapar,
+            )
+        )
+
+        # """Estimated annual maximum fAPAR (unitless)."""
+        # self.energy_limited: NDArray[np.bool_] = fapar_energylim < fapar_waterlim
+        # """Boolean array showing if annual :math:`fAPAR_{max}` is water or energy
+        # limited."""
+
+        self.lai_max: NDArray[np.floating] = -(1 / self.phenology_const.k) * np.log(
+            1.0 - self.fapar_max
+        )
+        """Estimated annual maximum LAI (unitless)"""
+
+        # self.lai_to_gpp_ratio_m = (
+        #     self.phenology_const.sigma
+        #     * self.annual_growing_season_length
+        #     * self.lai_max
+        # ) / (self.annual_total_potential_gpp * self.fapar_max)
+        # """The steady state ratio of leaf area index to potential GPP (:math:`m`)"""
+
+    def __repr__(self) -> str:
+        """Simple representation of class instance."""
+        return f"FaparLimitationNew(shape={self.shape})"
+
+    def summarize(self, dp: int = 2) -> None:
+        """Print summary of estimates of fAPAR limitation.
+
+        Prints a summary of the calculated values in a FaparLimitation instance
+        including the mean, range and number of nan values.
+
+        Args:
+            dp: The number of decimal places used in rounding summary stats.
+        """
+
+        summarize_attrs(self, self.limitation_method.attrs, dp=dp)
+
+    @classmethod
+    def from_pmodel(
+        cls,
+        pmodel: PModelABC,
+        growing_season: NDArray[np.bool],
+        precip: NDArray[np.floating],
+        datetimes: NDArray[np.datetime64] | None = None,
+        gpp_penalty_factor: NDArray[np.floating] | None = None,
+        method: str = "cai",
+        phenology_const: PhenologyConstNew = PhenologyConstNew(),
+        **kwargs: NDArray[np.floating],
+    ) -> FaparLimitationNew:
+        r"""Create a FaparLimitation instance from a P Model and other inputs.
+
+        The annual summary values of :math:`A_0, c_a, \chi` and :math:`D` used by the
+        :meth:`~pyrealm.phenology.fapar_limitation.FaparLimitation` class can be taken
+        directly from the predictions of a P Model. This method automatically extracts
+        the required data from a fitted P Model and returns a ``FaparLimitation``
+        instance.
+
+        .. NOTE::
+
+          The calculation of fAPAR limitation requires estimates of **potential** GPP,
+          so the :class:`~pyrealm.pmodel.pmodel_environment.PModelEnvironment` instance
+          used to fit the model **must** set ``fapar`` to be one.
+
+        Some additional information is needed:
+
+        * The calculation requires annual summaries of variables, so the ``datetimes``
+          argument must be used to provide an array of datetimes for each observation.
+
+        * The annual mean values :math:`c_a, \chi` and :math:`D` should be estimated
+          during the growing season, so the ``growing_season`` argument must be used to
+          provide an array of boolean values indicating which observations should be
+          treated as in the growing season.
+
+        * The calculation requires estimates of precipitation, so the ``precipitation``
+          argument must provide estimates of total precipitation during each
+          observations in moles of water per metre squared.
+
+        * The calculation of the :math:`f_0` parameter requires estimates of site
+          specific aridity index.
+
+        The method accepts both standard and subdaily P Models and automatically uses
+        the actual time intervals between observations to calculate the required
+        weighted annual means and sums. This might lead to unexpected values: the yearly
+        mean of monthly values :math:`1, 2, \dots, 12` would not be 6.5 because the
+        monthly values are weighted according to the length of the month.
+
+        Lastly, potential GPP is taken directly from the P Model instance. If you want
+        to apply a post-hoc penalty factor to GPP (e.g. a water limitation factor), then
+        you can optionally provide per-observation penalty estimates and they will be
+        applied when calculating annual total potential assimilation.
+
+        Args:
+            pmodel: A :class:`pyrealm.pmodel.pmodel.PModel` or
+                :class:`pyrealm.pmodel.pmodel.SubdailyPModel` instance, fitted with
+                ``fapar`` fixed at one.
+            datetimes: An array giving the datetimes of observations.
+            growing_season: A boolean array indicating which observations are to be
+                considered as part of the growing season.
+            precip: An array of precipitation for each observation.
+            gpp_penalty_factor: A post-hoc penalty factor to be applied to estimated
+                GPP.
+            method: The method to be used in calculating maximum fAPAR, defaulting to
+                `cai`.
+            phenology_const: An instance of
+                :class:`~pyrealm.constants.phenology_const.PhenologyConstNew`
+            **kwargs: Any additional variables required by specific method choices.
+        """
+
+        # Check the datetimes - should they be taken from the AcclimationModel of the
+        # SubdailyPModel or are they required for standard PModels?
+        if isinstance(pmodel, SubdailyPModel):
+            if datetimes is not None:
+                raise ValueError(
+                    "Observation datetimes are not required with SubdailyPModel "
+                    "inputs, the acclimation model datetimes are used."
+                )
+            datetimes = pmodel.acclim_model.datetimes
+
+        elif isinstance(pmodel, PModel):
+            if datetimes is None:
+                raise ValueError(
+                    "Observation datetimes are required with PModel inputs."
+                )
+
+        # Create the annual value calculator
+        # - the code above guards against datetimes being None
+        avc = AnnualValueCalculator(
+            data_shape=pmodel.shape,
+            timing=datetimes,  # type: ignore [arg-type]
+            subset_mask=growing_season,
+        )
+
+        # Get the total GPP for each observation
+        # - also need to handle missing values, easier to take _mean_ annual value
+        #   and scale it up to an annual total
+        # - TODO - handle incompleteness - when do we stop estimating annual values from
+        #   partial years (or at least warn about it)
+
+        # Extract GPP and apply any observation level penalty factor
+        total_gpp = pmodel.gpp
+        if gpp_penalty_factor is not None:
+            total_gpp *= gpp_penalty_factor
+
+        # Calculate annual mean potential GPP and scale up to the year
+        annual_mean_potential_gpp = avc.get_annual_means(total_gpp)
+        annual_total_potential_gpp = (
+            annual_mean_potential_gpp * (avc.year_n_days) * 86400 * 1e-6
+        ) / pmodel.core_const.k_c_molmass
+
+        # Calculate annual mean ca, chi and VPD within growing season
+        annual_mean_ca = avc.get_annual_means(pmodel.env.ca, within_subset=True)
+        annual_mean_chi = avc.get_annual_means(pmodel.optchi.chi, within_subset=True)
+        annual_mean_vpd = avc.get_annual_means(pmodel.env.vpd, within_subset=True)
+
+        # Calculate total annual precipitation
+        annual_total_precip = avc.get_annual_totals(precip)
+
+        return cls(
+            annual_total_potential_gpp=annual_total_potential_gpp,
+            annual_mean_ca=annual_mean_ca,
+            annual_mean_chi=annual_mean_chi,
+            annual_mean_vpd=annual_mean_vpd,
+            annual_total_precip=annual_total_precip,
+            annual_growing_season_length=avc.year_n_days_subset,
+            years=avc.years,
+            method=method,
+            phenology_const=phenology_const,
+        )
diff --git a/pyrealm/pmodel/pmodel_environment.py b/pyrealm/pmodel/pmodel_environment.py
index 43261768..0824e99c 100644
--- a/pyrealm/pmodel/pmodel_environment.py
+++ b/pyrealm/pmodel/pmodel_environment.py
@@ -194,7 +194,7 @@ def __init__(
 
         self._additional_vars: tuple[str, ...] = tuple(kwargs.keys())
         """A tuple containing the attribute names of additional variables passed to the
-        PModelEnivronment."""
+        PModelEnvironment."""
 
     def __repr__(self) -> str:
         """Generates a string representation of PModelEnvironment instance."""
diff --git a/pyrealm_build_data/phenology/zhu_method/zhu_method.py b/pyrealm_build_data/phenology/zhu_method/zhu_method.py
index c80d4fbd..12ab55c2 100644
--- a/pyrealm_build_data/phenology/zhu_method/zhu_method.py
+++ b/pyrealm_build_data/phenology/zhu_method/zhu_method.py
@@ -20,7 +20,9 @@ def cal_fapar_actually_in_python(
     fapar_carbon: NDArray[np.floating], fapar_water: NDArray[np.floating], budyko=4
 ) -> NDArray[np.floating]:
     """Refactor of plmodel_timeseries.cal_fapar using array calculations."""
-    cw_ratio = fapar_carbon / (np.clip(fapar_water, min=np.finfo(float).eps, max=None))
+    cw_ratio = fapar_carbon / (
+        np.clip(fapar_water, a_min=np.finfo(float).eps, a_max=None)
+    )
     return ((1 + cw_ratio) - (1 + cw_ratio**budyko) ** (1 / budyko)) * fapar_water
 
 


Test changes:
diff --git a/tests/broadcasting/utils.py b/tests/broadcasting/utils.py
index 9ce4a79d..d803e836 100644
--- a/tests/broadcasting/utils.py
+++ b/tests/broadcasting/utils.py
@@ -61,7 +61,9 @@
     PlantFunctionalTypeStrict,
     StemTraits,
 )
-from pyrealm.phenology.fapar_limitation import FaparLimitation, PhenologyConst
+from pyrealm.phenology.fapar_limitation import FaparLimitation
+
+# from pyrealm.phenology.fapar_limitation_new import FaparLimitationNew
 
 ## Lists / functions to manually define arguments, methods or outputs to ignore, etc.
 
@@ -105,11 +107,15 @@
 ]
 
 
-# Ignore these outputs, they are not expected to be equal.
 # Formats: [fn name] for function results, [class]:[attr] for class attributes
 IGNORE_OUTPUTS = [
+    # Ignore these outputs, they are not expected to be equal.
     "Cohorts:_cohort_id",
     "Calendar:n_dates",
+    # Causes infinite recursion, because this is an attribute of FaparLimitationNew, but
+    # contains a reference to the parent FaparLimitationNew, so
+    # is_equal/compare_instances runs away.
+    "FaparLimitationNew:limitation_method",
 ]
 
 # The REQUIRES dictionary provides a way of populating required keywords arguments for
@@ -130,7 +136,7 @@ def kwarg_params(names: tuple[str, ...]) -> dict[str, Parameter]:
     }
 
 
-REQUIRES: dict[tuple[str, tuple[str, ...]], dict[str, Parameter]] = {
+REQUIRES: dict[tuple[str, tuple[str, ...]], dict[str, Any]] = {
     (
         "PModelEnvironment",
         (
@@ -177,6 +183,17 @@ def kwarg_params(names: tuple[str, ...]) -> dict[str, Parameter]:
             "QuantumYieldSandoval",
         ),
     ): kwarg_params(("aridity_index", "mean_growth_temperature")),
+    (
+        "FaparLimitationNew",
+        (
+            "FaparLimitationCai.calculate_maximum_fapar",
+            "FaparLimitationCai",
+        ),
+    ): kwarg_params(("aridity_index",)),
+    (
+        "FaparLimitationNew",
+        ("FaparLimitationNew",),
+    ): kwarg_params(("aridity_index",)),  # Default method cai requires aridity_index
 }
 
 
@@ -206,18 +223,16 @@ def defined_method_args(argument: str, ctx: Context) -> Any | None:
     n_heights = 2
     pft_names = [f"Tree{i + 1}" for i in range(n_pft)]
 
+    # Fapar Limitation for two years
     fapar_limitation = FaparLimitation(
-        annual_total_potential_gpp=np.ones(48),
-        annual_mean_ca=np.ones(48),
-        annual_mean_chi=np.ones(48),
-        annual_mean_vpd=np.ones(48),
-        annual_total_precip=np.ones(48),
-        annual_growing_season_length=np.ones(48),
-        aridity_index=np.ones(48),
-        years=np.zeros((48,), dtype="datetime64[Y]"),
-        phenology_const=PhenologyConst(
-            z=12.227, k=0.5, f0_coefficients=(0.65, 0.604169, 1.9), sigma=0.771
-        ),
+        annual_total_potential_gpp=np.ones(2),
+        annual_mean_ca=np.ones(2),
+        annual_mean_chi=np.ones(2),
+        annual_mean_vpd=np.ones(2),
+        annual_total_precip=np.ones(2),
+        annual_growing_season_length=np.ones(2),
+        aridity_index=np.ones(2),
+        years=np.array(["2000", "2001"], dtype="datetime64[Y]"),
     )
 
     method_arguments_list: dict[str, dict] = {
@@ -271,8 +286,12 @@ def defined_method_args(argument: str, ctx: Context) -> Any | None:
         ## Phenology
         "FaparLimitation": {"years": np.ones(bcast_shape[0], dtype="datetime64[Y]")},
         "Phenology": {
-            "daily_gpp": np.full((48,), 0.5),
-            "datetimes": np.arange(0, 48, dtype="datetime64[D]"),
+            "daily_gpp": np.full((31,), 0.5),
+            "datetimes": np.arange(
+                np.datetime64("2000-01-01"),
+                np.datetime64("2000-02-01"),
+                np.timedelta64(1, "D"),
+            ),
             "fapar_limitation": fapar_limitation,
         },
     }
@@ -326,6 +345,7 @@ def defined_method_args(argument: str, ctx: Context) -> Any | None:
             nTime = 366
         else:
             nTime = splashDatesLen
+
         splashShape = (1 if shape[0] == 1 else nTime, *shape[1:])
         arguments = {
             "dates": Calendar(np.arange(0, nTime, dtype="datetime64[D]")),
@@ -336,18 +356,30 @@ def defined_method_args(argument: str, ctx: Context) -> Any | None:
             "pn": np.full(splashShape, 10),
         }
 
-    if ctx.name == "PModelEnvironment":
-        if ctx.parents and ctx.parents[-1] == "SubdailyPModel":
-            # SubdailyPModel needs more than 1 day (uses 48 hourly times)
-            envShape = (1 if shape[0] == 1 else 48, *shape[1:])
-            arguments = {
-                "tc": np.full(envShape, 20),
-                "vpd": np.full(envShape, 40),
-                "co2": np.full(envShape, 1000),
-                "patm": np.full(envShape, 101325),
-                "fapar": np.full(envShape, 1),
-                "ppfd": np.full(envShape, 800),
-            }
+    if (  # SubdailyPModel needs more than 1 day (uses 48 hourly times)
+        ctx.name == "PModelEnvironment"
+        and ctx.parents
+        and ctx.parents[-1] == "SubdailyPModel"
+    ):
+        envShape = (1 if shape[0] == 1 else 48, *shape[1:])
+        arguments = {
+            "tc": np.full(envShape, 20),
+            "vpd": np.full(envShape, 40),
+            "co2": np.full(envShape, 1000),
+            "patm": np.full(envShape, 101325),
+            "fapar": np.full(envShape, 1),
+            "ppfd": np.full(envShape, 800),
+        }
+
+    # FaparLimitationNew requires years as Datetime64[Y]
+    if ctx.name == "FaparLimitationNew":
+        arguments = {
+            "years": np.ones(bcast_shape[0], dtype="datetime64[Y]"),
+        }
+
+        # Need to enforce the method selection for FaparLimitationZhu
+        if ctx.parents and ctx.parents[-1] == "FaparLimitationZhu":
+            arguments["method"] = "zhu"
 
     return arguments.get(argument)
 
diff --git a/tests/regression/phenology/test_fapar_limitation_new.py b/tests/regression/phenology/test_fapar_limitation_new.py
new file mode 100644
index 00000000..49bd96f4
--- /dev/null
+++ b/tests/regression/phenology/test_fapar_limitation_new.py
@@ -0,0 +1,104 @@
+"""Test the FaparLimitation class."""
+
+import json
+from importlib import resources
+
+import pandas as pd
+import pytest
+from numpy.testing import assert_allclose
+
+
+def dataframe_to_dict_of_nparrays(data):
+    """Utility function to preconvert a dataframe of pd.Series to np.array."""
+
+    data_dict = {k: v.to_numpy() for k, v in data.items()}
+
+    return data_dict
+
+
+@pytest.fixture
+def site_data():
+    """Load the site data."""
+
+    datafile = (
+        resources.files("pyrealm_build_data.phenology.inputs.source")
+        / "DE-GRI_site_data.json"
+    )
+
+    with open(str(datafile)) as json_src:
+        site_data = json.load(json_src)
+
+    return site_data
+
+
+@pytest.fixture()
+def data_fapar_limitation(timescale, method_predictions_dir):
+    """Load the input data from csv file."""
+
+    inputs = pd.read_csv(
+        str(
+            resources.files(f"pyrealm_build_data.phenology.inputs.{timescale}")
+            / "annual_inputs.csv"
+        )
+    )
+
+    predictions = pd.read_csv(
+        str(
+            resources.files(f"pyrealm_build_data.phenology.{method_predictions_dir}")
+            / "fapar_max_predictions.csv"
+        )
+    )
+
+    data = inputs.merge(predictions)
+    return dataframe_to_dict_of_nparrays(data)
+
+
+@pytest.mark.parametrize(
+    argnames="timescale,timescale_abbr,assim_var",
+    argvalues=(
+        pytest.param("fortnightly", "ft", "annual_total_A0", id="fortnightly"),
+        pytest.param("subdaily", "hh", "annual_total_A0_smstress", id="subdaily"),
+    ),
+)
+@pytest.mark.parametrize(
+    argnames="method_predictions_dir, method",
+    argvalues=(
+        pytest.param("cai_zhou_method", "cai", id="cai"),
+        pytest.param("zhu_method", "zhu", id="zhu"),
+    ),
+)
+def test_faparlimitation(
+    site_data,
+    data_fapar_limitation,
+    timescale,  # parameterises data_fapar_limitation fixture
+    timescale_abbr,
+    assim_var,
+    method_predictions_dir,  # parameterises data_fapar_limitation fixture
+    method,
+):
+    """Regression test for FaparLimitation constructor with fortnightly data."""
+
+    from pyrealm.phenology.fapar_limitation_new import FaparLimitationNew
+
+    faparlim = FaparLimitationNew(
+        annual_total_potential_gpp=data_fapar_limitation[assim_var],
+        annual_mean_ca=data_fapar_limitation["annual_mean_ca_in_GS"],
+        annual_mean_chi=data_fapar_limitation["annual_mean_chi_in_GS"],
+        annual_mean_vpd=data_fapar_limitation["annual_mean_VPD_in_GS"],
+        annual_total_precip=data_fapar_limitation["annual_precip_molar"],
+        annual_growing_season_length=data_fapar_limitation["N_growing_days"],
+        years=data_fapar_limitation["year"].astype(str).astype("datetime64[Y]"),
+        method=method,
+        aridity_index=site_data["AI_from_cruts"],  # Not used by zhu method.
+    )
+
+    assert_allclose(
+        data_fapar_limitation[f"fapar_max_{timescale_abbr}"],
+        faparlim.fapar_max,
+        rtol=1e-6,
+    )
+    assert_allclose(
+        data_fapar_limitation[f"lai_max_{timescale_abbr}"],
+        faparlim.lai_max,
+        rtol=1e-6,
+    )