Added new stomatal model from Buckley et al 2026

plugins/photosynthesis/doc/Photosynthesis.dox

@@ -57,6 +57,8 @@
  <tr><td>Gamma\_CO2</td><td>\f$\Gamma\f$</td><td>\f$\mu\f$mol CO<sub>2</sub>/mol</td><td>\htmlonly<span style="font-family: Courier, monospace; color: green;">float</span>\endhtmlonly</td><td>Farquhar (C<sub>3</sub>)</td><td>Photosynthetic CO<sub>2</sub> compensation point (including "dark respiration").</td></tr>
  <tr><td>limitation\_state</td><td>N/A</td><td>N/A</td><td>\htmlonly<span style="font-family: Courier, monospace; color: green;">int</span>\endhtmlonly</td><td>Farquhar (C<sub>3</sub>), C<sub>4</sub></td><td>Photosynthesis limitation state. For the C<sub>3</sub> Farquhar model: 0 = Rubisco-limited, 1 = electron-transport-limited. For the C<sub>4</sub> model: 1 = enzyme (\f$A_c\f$) limited, 2 = electron-transport (\f$A_j\f$) limited.</td></tr>
  <tr><td>electron\_transport\_ratio</td><td>\f$J/J_{max}\f$</td><td>unitless</td><td>\htmlonly<span style="font-family: Courier, monospace; color: green;">float</span>\endhtmlonly</td><td>Farquhar (C<sub>3</sub>)</td><td>Ratio of actual electron transport rate to maximum electron transport rate. Ranges from 0 (dark) to ~1 (light-saturated). Useful for computing fluorescence quantum yield.</td></tr>
+ <tr><td>electron\_transport\_rate</td><td>\f$J\f$</td><td>\f$\mu\f$mol/m<sup>2</sup>-sec</td><td>\htmlonly<span style="font-family: Courier, monospace; color: green;">float</span>\endhtmlonly</td><td>Farquhar (C<sub>3</sub>)</td><td>Potential electron transport rate at the leaf's light and temperature, the smaller root of \f$\theta J^2-(\alpha Q+J_{max})J+\alpha QJ_{max}=0\f$. Does not depend on \f$C_i\f$. Written automatically when the Buckley et al. (2026) stomatal conductance model is in use.</td></tr>
+ <tr><td>dark\_respiration</td><td>\f$R_d\f$</td><td>\f$\mu\f$mol CO<sub>2</sub>/m<sup>2</sup>-sec</td><td>\htmlonly<span style="font-family: Courier, monospace; color: green;">float</span>\endhtmlonly</td><td>Farquhar (C<sub>3</sub>)</td><td>Day respiration at the leaf temperature. Written automatically when the Buckley et al. (2026) stomatal conductance model is in use.</td></tr>
  <tr><td>Cm</td><td>\f$C_m\f$</td><td>\f$\mu\f$bar</td><td>\htmlonly<span style="font-family: Courier, monospace; color: green;">float</span>\endhtmlonly</td><td>C<sub>4</sub></td><td>Mesophyll cytosolic CO<sub>2</sub> partial pressure, \f$C_m = C_i - A/g_m\f$.</td></tr>
  <tr><td>Vp</td><td>\f$V_p\f$</td><td>\f$\mu\f$mol CO<sub>2</sub>/m<sup>2</sup>/s</td><td>\htmlonly<span style="font-family: Courier, monospace; color: green;">float</span>\endhtmlonly</td><td>C<sub>4</sub></td><td>PEP carboxylation rate.</td></tr>
  </table>

plugins/photosynthesis/include/PhotosynthesisModel.h

@@ -911,7 +911,11 @@ class PhotosynthesisModel {
     //! Add optional output primitive data values to the Context
     /**
      * \param[in] label Name of primitive data. Available labels for the Farquhar model: "Ci" (intercellular CO2), "limitation_state" (0=Rubisco-limited, 1=electron transport-limited), "Gamma_CO2" (CO2 compensation point), "electron_transport_ratio"
-     * (relative light saturation Ja/Je of van der Tol et al. 2014: electron transport used by carbon metabolism divided by its light-limited potential; 1 in the dark and falling as light saturates photosynthesis; drives the RadiationModel SIF fluorescence yield). Additional labels for the C4 model: "Ci", "limitation_state" (1=enzyme-limited, 2=electron-transport-limited), "Cm" (mesophyll cytosolic CO2), "Vp" (PEP carboxylation rate), "electron_transport_ratio".
+     * (relative light saturation Ja/Je of van der Tol et al. 2014: electron transport used by carbon metabolism divided by its light-limited potential; 1 in the dark and falling as light saturates photosynthesis; drives the RadiationModel SIF
+     * fluorescence yield), "electron_transport_rate" (potential electron transport rate J at the leaf's light and temperature, umol/m^2/s), "dark_respiration" (day respiration Rd at the leaf temperature, umol/m^2/s). Additional labels for the C4
+     * model: "Ci", "limitation_state" (1=enzyme-limited, 2=electron-transport-limited), "Cm" (mesophyll cytosolic CO2), "Vp" (PEP carboxylation rate), "electron_transport_ratio".
+     * \note "electron_transport_rate" and "dark_respiration" are also written automatically for the Farquhar model when the StomatalConductanceModel Buckley et al. (2026) model is in use (global data "buckley2026_stomatal_model_enabled" is 1), since
+     * that model needs them.
      */
     void optionalOutputPrimitiveData(const char *label);
 
@@ -992,6 +996,17 @@ class PhotosynthesisModel {
 
     static float respondToTemperature(const PhotosyntheticTemperatureResponseParameters *temperatureResponseParameters, float T);
 
+    //! Potential electron transport rate and day respiration of the Farquhar model at a given light and leaf temperature
+    /**
+     * Neither depends on the intercellular CO2 concentration, so they are evaluated directly rather than by the Ci solver.
+     * \param[in] params Farquhar model coefficients
+     * \param[in] i_PAR Incident PAR photon flux density (umol/m^2/s)
+     * \param[in] TL Leaf temperature (K)
+     * \param[out] J Potential electron transport rate (umol/m^2/s): the smaller root of theta*J^2 - (alpha*Q + Jmax)*J + alpha*Q*Jmax = 0
+     * \param[out] Rd Day respiration (umol/m^2/s)
+     */
+    static void evaluateElectronTransportAndRespiration(const FarquharModelCoefficients &params, float i_PAR, float TL, double &J, float &Rd);
+
     float i_PAR_default;
     float TL_default;
     float CO2_default;

plugins/photosynthesis/src/PhotosynthesisModel.cpp

@@ -949,6 +949,17 @@ void PhotosynthesisModel::run(const std::vector<uint> &lUUIDs) {
     WarningAggregator warnings;
     warnings.setEnabled(message_flag); // Respect existing message flag
 
+    // The electron transport rate and day respiration are written when requested as optional output, and also whenever the StomatalConductanceModel Buckley et al. (2026) model is in use,
+    // which signals this through global data and computes the CO2-saturated electron-transport-limited rate J/4 - Rd from them.
+    bool buckley_inputs_requested = false;
+    if (context->doesGlobalDataExist("buckley2026_stomatal_model_enabled") && context->getGlobalDataType("buckley2026_stomatal_model_enabled") == HELIOS_TYPE_INT) {
+        int buckley_enabled = 0;
+        context->getGlobalData("buckley2026_stomatal_model_enabled", buckley_enabled);
+        buckley_inputs_requested = (buckley_enabled == 1);
+    }
+    const bool write_electron_transport_rate = buckley_inputs_requested || std::find(output_prim_data.begin(), output_prim_data.end(), "electron_transport_rate") != output_prim_data.end();
+    const bool write_dark_respiration = buckley_inputs_requested || std::find(output_prim_data.begin(), output_prim_data.end(), "dark_respiration") != output_prim_data.end();
+
     for (uint UUID: lUUIDs) {
 
         float i_PAR;
@@ -1164,6 +1175,18 @@ void PhotosynthesisModel::run(const std::vector<uint> &lUUIDs) {
                 context->setPrimitiveData(UUID, "Vp", Vp_c4);
             }
         }
+
+        if (model == "farquhar" && (write_electron_transport_rate || write_dark_respiration)) {
+            double J;
+            float Rd;
+            evaluateElectronTransportAndRespiration(getCoefficientsForPrimitive_Farquhar(UUID), i_PAR, TL, J, Rd);
+            if (write_electron_transport_rate) {
+                context->setPrimitiveData(UUID, "electron_transport_rate", static_cast<float>(J));
+            }
+            if (write_dark_respiration) {
+                context->setPrimitiveData(UUID, "dark_respiration", Rd);
+            }
+        }
     }
 
     warnings.report(std::cerr);
@@ -1293,7 +1316,7 @@ float PhotosynthesisModel::evaluateCi_Farquhar(float Ci, std::vector<float> &var
 
     const FarquharModelCoefficients modelcoeffs = *reinterpret_cast<const FarquharModelCoefficients *>(parameters);
 
-    float Vcmax, Jmax, TPU, Rd, alpha, theta;
+    float Vcmax, TPU, Rd, alpha;
 
     float CO2 = variables[0];
     float Q = variables[1];
@@ -1318,17 +1341,9 @@ float PhotosynthesisModel::evaluateCi_Farquhar(float Ci, std::vector<float> &var
         Vcmax = respondToTemperature(&modelcoeffs.VcmaxTempResponse, TL);
     }
 
-    if (modelcoeffs.Jmax > 0) {
-        Jmax = modelcoeffs.Jmax * expf(modelcoeffs.dH_Jmax * (invDiffRT));
-    } else {
-        Jmax = respondToTemperature(&modelcoeffs.JmaxTempResponse, TL);
-    }
-
-    if (modelcoeffs.Rd > 0) {
-        Rd = modelcoeffs.Rd * expf(modelcoeffs.dH_Rd * (invDiffRT));
-    } else {
-        Rd = respondToTemperature(&modelcoeffs.RdTempResponse, TL);
-    }
+    // Electron transport rate J and day respiration Rd do not depend on Ci
+    double J;
+    evaluateElectronTransportAndRespiration(modelcoeffs, Q, TL, J, Rd);
 
     if (modelcoeffs.alpha > 0) {
         alpha = modelcoeffs.alpha;
@@ -1338,7 +1353,6 @@ float PhotosynthesisModel::evaluateCi_Farquhar(float Ci, std::vector<float> &var
 
 
     TPU = respondToTemperature(&modelcoeffs.TPUTempResponse, TL);
-    theta = respondToTemperature(&modelcoeffs.thetaTempResponse, TL);
 
     // Mesophyll conductance gm (mol CO2 / m^2 / s / bar). If gm is infinite/huge,
     // Cc ≡ Ci and we take the original legacy path below. Otherwise we solve a
@@ -1353,14 +1367,6 @@ float PhotosynthesisModel::evaluateCi_Farquhar(float Ci, std::vector<float> &var
     float Ko = exp(c_Ko - dH_Ko / (R * TL));
     float Kco = Kc * (1.f + O / Ko);
 
-    // thetaJ^2 + -(alphaQ+Jmax)J + alphaQJmax = 0
-    double a = std::max(theta, 0.0001f);
-    double ia = 1.000f / a;
-    double b = -(alpha * Q + Jmax);
-    double c = alpha * Q * Jmax;
-    double J = (-b - sqrt(pow(b, 2.000f) - 4.f * a * c)) * 0.5f * ia;
-    // J = Jmax * alpha * Q / (alpha * Q + Jmax);
-
     float A;
     float limitation_state;
 
@@ -1693,6 +1699,40 @@ float PhotosynthesisModel::evaluateC4Model(const C4ModelCoefficients &params, fl
 }
 
 
+void PhotosynthesisModel::evaluateElectronTransportAndRespiration(const FarquharModelCoefficients &params, float i_PAR, float TL, double &J, float &Rd) {
+
+    const float R = 0.0083144598; // molar gas constant (kJ/K/mol)
+    const float invDiffRT = (1.f / 298.15f - 1.f / TL) / R;
+
+    float Jmax;
+    if (params.Jmax > 0) {
+        Jmax = params.Jmax * expf(params.dH_Jmax * (invDiffRT));
+    } else {
+        Jmax = respondToTemperature(&params.JmaxTempResponse, TL);
+    }
+
+    if (params.Rd > 0) {
+        Rd = params.Rd * expf(params.dH_Rd * (invDiffRT));
+    } else {
+        Rd = respondToTemperature(&params.RdTempResponse, TL);
+    }
+
+    float alpha;
+    if (params.alpha > 0) {
+        alpha = params.alpha;
+    } else {
+        alpha = respondToTemperature(&params.alphaTempResponse, TL);
+    }
+    const float theta = respondToTemperature(&params.thetaTempResponse, TL);
+
+    // thetaJ^2 + -(alphaQ+Jmax)J + alphaQJmax = 0
+    double a = std::max(theta, 0.0001f);
+    double ia = 1.000f / a;
+    double b = -(alpha * i_PAR + Jmax);
+    double c = alpha * i_PAR * Jmax;
+    J = (-b - sqrt(pow(b, 2.000f) - 4.f * a * c)) * 0.5f * ia;
+}
+
 float PhotosynthesisModel::respondToTemperature(const PhotosyntheticTemperatureResponseParameters *params, float T_in_Kelvin) {
     float T = T_in_Kelvin;
     float R = 0.0083144598f;
@@ -1779,7 +1819,8 @@ void PhotosynthesisModel::enableMessages() {
 
 void PhotosynthesisModel::optionalOutputPrimitiveData(const char *label) {
 
-    if (strcmp(label, "Ci") == 0 || strcmp(label, "limitation_state") == 0 || strcmp(label, "Gamma_CO2") == 0 || strcmp(label, "electron_transport_ratio") == 0 || strcmp(label, "Cm") == 0 || strcmp(label, "Vp") == 0) {
+    if (strcmp(label, "Ci") == 0 || strcmp(label, "limitation_state") == 0 || strcmp(label, "Gamma_CO2") == 0 || strcmp(label, "electron_transport_ratio") == 0 || strcmp(label, "Cm") == 0 || strcmp(label, "Vp") == 0 ||
+        strcmp(label, "electron_transport_rate") == 0 || strcmp(label, "dark_respiration") == 0) {
         output_prim_data.emplace_back(label);
     } else {
         if (message_flag) {

plugins/photosynthesis/tests/selfTest.cpp

@@ -2357,6 +2357,80 @@ DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar library values match documentati
     DOCTEST_CHECK(pist_male.getVcmaxTempResponse().value_at_25C != doctest::Approx(pist_female.getVcmaxTempResponse().value_at_25C).epsilon(1e-4));
 }
 
+DOCTEST_TEST_CASE("PhotosynthesisModel - electron_transport_rate and dark_respiration optional outputs") {
+    Context context;
+    PhotosynthesisModel photosynthesis(&context);
+    photosynthesis.disableMessages();
+
+    FarquharModelCoefficients coeffs;
+    coeffs.setVcmax(100.f);
+    coeffs.setJmax(150.f);
+    coeffs.setRd(1.5f);
+    coeffs.setQuantumEfficiency_alpha(0.3f);
+    coeffs.setLightResponseCurvature_theta(0.7f);
+    photosynthesis.setModelCoefficients(coeffs);
+
+    const float PAR = 200.f; // W/m^2
+    uint leaf = context.addPatch();
+    context.setPrimitiveData(leaf, "radiation_flux_PAR", PAR);
+    context.setPrimitiveData(leaf, "temperature", 298.15f);
+    context.setPrimitiveData(leaf, "moisture_conductance", 0.3f);
+
+    // Not written unless requested
+    photosynthesis.run();
+    DOCTEST_CHECK(!context.doesPrimitiveDataExist(leaf, "electron_transport_rate"));
+    DOCTEST_CHECK(!context.doesPrimitiveDataExist(leaf, "dark_respiration"));
+
+    photosynthesis.optionalOutputPrimitiveData("electron_transport_rate");
+    photosynthesis.optionalOutputPrimitiveData("dark_respiration");
+    photosynthesis.run();
+
+    // J is the smaller root of theta J^2 - (alpha Q + Jmax) J + alpha Q Jmax = 0, with the temperature-independent coefficients above
+    const double aQ = 0.3 * PAR * 4.57, Jmax = 150.0, theta = 0.7;
+    const double J_expected = ((aQ + Jmax) - std::sqrt((aQ + Jmax) * (aQ + Jmax) - 4.0 * theta * aQ * Jmax)) / (2.0 * theta);
+    float J, Rd;
+    context.getPrimitiveData(leaf, "electron_transport_rate", J);
+    context.getPrimitiveData(leaf, "dark_respiration", Rd);
+    DOCTEST_CHECK(J == doctest::Approx(J_expected).epsilon(1e-4));
+    DOCTEST_CHECK(Rd == doctest::Approx(1.5f).epsilon(1e-6));
+
+    // J and Rd do not depend on Ci, so J/4 - Rd bounds the electron-transport-limited net assimilation from above
+    float A;
+    context.getPrimitiveData(leaf, "net_photosynthesis", A);
+    DOCTEST_CHECK(A < 0.25f * J - Rd);
+}
+
+DOCTEST_TEST_CASE("PhotosynthesisModel - electron_transport_rate and dark_respiration are written when the Buckley et al. (2026) stomatal model requests them") {
+    Context context;
+    PhotosynthesisModel photosynthesis(&context);
+    photosynthesis.disableMessages();
+    photosynthesis.setModelCoefficients(photosynthesis.getFarquharCoefficientsFromLibrary("Almond"));
+
+    uint leaf = context.addPatch();
+    context.setPrimitiveData(leaf, "radiation_flux_PAR", 300.f);
+    context.setPrimitiveData(leaf, "temperature", 303.15f);
+
+    context.setGlobalData("buckley2026_stomatal_model_enabled", 0);
+    photosynthesis.run();
+    DOCTEST_CHECK(!context.doesPrimitiveDataExist(leaf, "electron_transport_rate"));
+
+    context.setGlobalData("buckley2026_stomatal_model_enabled", 1);
+    photosynthesis.run();
+    float J, Rd;
+    DOCTEST_REQUIRE(context.doesPrimitiveDataExist(leaf, "electron_transport_rate"));
+    DOCTEST_REQUIRE(context.doesPrimitiveDataExist(leaf, "dark_respiration"));
+    context.getPrimitiveData(leaf, "electron_transport_rate", J);
+    context.getPrimitiveData(leaf, "dark_respiration", Rd);
+    DOCTEST_CHECK(J > 0.f);
+    DOCTEST_CHECK(Rd > 0.2f); // almond Rd25 = 0.2, higher at 30 C
+
+    // In the dark J is zero
+    context.setPrimitiveData(leaf, "radiation_flux_PAR", 0.f);
+    photosynthesis.run();
+    context.getPrimitiveData(leaf, "electron_transport_rate", J);
+    DOCTEST_CHECK(J == doctest::Approx(0.f).epsilon(1e-6));
+}
+
 int PhotosynthesisModel::selfTest(int argc, char **argv) {
     return helios::runDoctestWithValidation(argc, argv);
 }

plugins/planthydraulics/doc/PlantHydraulics.dox

@@ -41,8 +41,10 @@
 
 | Primitive Data Label | Symbol | Units | Data Type | Description |
 | -------------------- | ------ | ------ | --------- | ----------- |
-| hydraulic_conductance  | \f$K\f$   | mol m\f$^{-2}\f$ s\f$^{-1}\f$ MPa\f$^{-1}\f$ | \htmlonly<span style="color: green;"><tt>float</tt></span>\endhtmlonly    | Hydraulic conductance |
-| hydraulic_capacitance | \f$C\f$    | mol m\f$^{-2}\f$ MPa\f$^{-1}\f$ | \htmlonly<span style="color: green;"><tt>float</tt></span>\endhtmlonly | Hydraulic capacitance |
+| hydraulic_conductance  | \f$K_{leaf}\f$   | mol m\f$^{-2}\f$ s\f$^{-1}\f$ MPa\f$^{-1}\f$ | \htmlonly<span style="color: green;"><tt>float</tt></span>\endhtmlonly    | Leaf hydraulic conductance. Enabled with PlantHydraulicsModel::outputConductancePrimitiveData(). |
+| hydraulic_capacitance | \f$C\f$    | mol m\f$^{-2}\f$ MPa\f$^{-1}\f$ | \htmlonly<span style="color: green;"><tt>float</tt></span>\endhtmlonly | Hydraulic capacitance. Enabled with PlantHydraulicsModel::outputCapacitancePrimitiveData(). |
+| plant_hydraulic_conductance | \f$K_{plant}\f$ | mol m\f$^{-2}\f$ s\f$^{-1}\f$ MPa\f$^{-1}\f$ | \htmlonly<span style="color: green;"><tt>float</tt></span>\endhtmlonly | Whole-plant (soil-to-leaf) hydraulic conductance per unit leaf area, \f$K_{plant}=\left(K_{root}^{-1}+K_{stem}^{-1}+K_{leaf}^{-1}\right)^{-1}\f$. Enabled with PlantHydraulicsModel::outputPlantConductancePrimitiveData(), or automatically when the Buckley et al. (2026) stomatal conductance model is in use. |
+| soil_water_potential | \f$\psi_{soil}\f$ | MPa | \htmlonly<span style="color: green;"><tt>float</tt></span>\endhtmlonly | Soil water potential of the leaf's plant. Written together with plant_hydraulic_conductance, but only for plants whose soil water potential was set with PlantHydraulicsModel::setSoilWaterPotentialOfPlant(), so that a value set directly as primitive data is not replaced by the default placeholder. |
 
 \section PH_Intro Introduction
 The Plant Hydraulics plugin models the water potentials, water contents, hydraulic conductances, and