Plant Hydraulics

function [flux, root] = brent_root (func, physcon, atmos, leaf, flux, xa, xb, tol)

% Use Brent's method to find the root of a function, which is known to exist between
% xa and xb. The root is updated until its accuracy is tol. func is the name of the
% function to solve. The variable root is returned as the root of the function. The
% function being evaluated has the definition statement: 
%
% function [flux, fx] = func (physcon, atmos, leaf, flux, x)
%
% The function func is exaluated at x and the returned value is fx. It uses variables
% in the physcon, atmos, leaf, and flux structures. These are passed in as
% input arguments. It also calculates values for variables in the flux structure
% so this must be returned in the function call as an output argument. The matlab
% function feval evaluates func.

% --- Evaluate func at xa and xb and make sure the root is bracketed

a = xa;
b = xb;
[flux, fa] = feval(func, physcon, atmos, leaf, flux, a);
[flux, fb] = feval(func, physcon, atmos, leaf, flux, b);

if ((fa > 0 & fb > 0) | (fa < 0 & fb < 0))
   error('brent_root error: root must be bracketed')
end

% --- Initialize iteration

itmax = 50;      % Maximum number of iterations
eps1 = 1e-08;    % Relative error tolerance

c = b;
fc = fb;

% --- Iterative root calculation

for iter = 1:itmax
   if ((fb > 0 & fc > 0) | (fb < 0 & fc < 0))
      c = a;
      fc = fa;
      d = b - a;
      e = d;
   end
   if (abs(fc) < abs(fb))
      a = b;
      b = c;
      c = a;
      fa = fb;
      fb = fc;
      fc = fa;
   end
   tol1 = 2 * eps1 * abs(b) + 0.5 * tol;
   xm = 0.5 * (c - b);

   % Check to end iteration

   if (abs(xm) <= tol1 | fb == 0)
      break
   end

   if (abs(e) >= tol1 & abs(fa) > abs(fb))
      s = fb / fa;
      if (a == c)
         p = 2 * xm * s;
         q = 1 - s;
      else
         q = fa / fc;
         r = fb / fc;
         p = s * (2 * xm * q * (q - r) - (b - a) * (r - 1));
         q = (q - 1) * (r - 1) * (s - 1);
      end
      if (p > 0)
         q = -q;
      end
      p = abs(p);
      if (2*p < min(3*xm*q-abs(tol1*q), abs(e*q)))
         e = d;
         d = p / q;
      else
         d = xm;
         e = d;
      end
   else
      d = xm;
      e = d;
   end
   a = b;
   fa = fb;
   if (abs(d) > tol1)
      b = b + d;
   else
      if (xm >= 0)
         b = b + abs(tol1);
      else
         b = b - abs(tol1);
      end
   end
   [flux, fb] = feval(func, physcon, atmos, leaf, flux, b);

   % Check to end iteration

   if (fb == 0)
      break
   end

   % Check to see if failed to converge

   if (iter == itmax)
      error('brent_root error: Maximum number of interations exceeded')
   end

end
root = b;

function [flux] = CiFuncOptimization (atmos, leaf, flux)

% Calculate leaf photosynthesis for a specified stomatal conductance.
% Then calculate Ci from the diffusion equation. 
%
% This routine uses a quadratic equation to solve for net photosynthesis (An).
% A general equation for C3 photosynthesis is:
%
%      a*(Ci - Cp)
% An = ----------- - Rd
%       e*Ci + d
%
% where:
%
% An = Net leaf photosynthesis (umol CO2/m2/s)
% Rd = Leaf respiration (umol CO2/m2/s)
% Ci = Intercellular CO2 concentration (umol/mol)
% Cp = CO2 compensation point (umol/mol)
% 
% Rubisco-limited photosynthesis (Ac)
% a  = Vcmax
% e  = 1
% d  = Kc*(1 + Oi/Ko)
%
% RuBP regeneration-limited photosynthesis (Aj)
% a = J
% e = 4
% d = 8*Cp
%
% where:
%
% Vcmax = Maximum carboxylation rate (umol/m2/s)
% Kc    = Michaelis-Menten constant for CO2 (umol/mol)
% Ko    = Michaelis-Menten constant for O2 (mmol/mol)
% Oi    = Intercellular O2 concentration (mmol/mol)
% J     = Electron transport rate (umol/m2/s)
%
% Ci is calculated from the diffusion equation:
%
%                   1.4   1.6
% An = (Ca - Ci) / (--- + ---)
%                   gb    gs
%
%            1.4   1.6
% Ci = Ca - (--- + ---)*An
%            gb    gs
%
% where:
% 
% Ca  = Atmospheric CO2 concentration (umol/mol)
% gb  = Leaf boundary layer conductance (mol H2O/m2/s)
% gs  = Leaf stomatal conductance (mol H2O/m2/s)
% 1.4 = Corrects gb for the diffusivity of CO2 compared with H2O
% 1.6 = Corrects gs for the diffusivity of CO2 compared with H2O
%
% The resulting quadratic equation is: a*An^2 + b*An + c = 0, which
% is solved for An. Correct solution is the smaller of the two roots.
%
% A similar approach is used for C4 photosynthesis.

% ------------------------------------------------------
% Input
%   atmos.o2air         ! Atmospheric O2 (mmol/mol)
%   atmos.co2air        ! Atmospheric CO2 (umol/mol)
%   leaf.c3psn          ! Photosynthetic pathway: 1 = C3. 0 = C4 plant
%   leaf.colim          ! Photosynthesis co-limitation: 0 = no. 1 = yes
%   leaf.colim_c3       ! Empirical curvature parameter for C3 co-limitation
%   leaf.colim_c4a      ! Empirical curvature parameter for C4 co-limitation
%   leaf.colim_c4b      ! Empirical curvature parameter for C4 co-limitation
%   leaf.qe_c4          ! C4: Quantum yield (mol CO2 / mol photons)
%   flux.vcmax          ! Maximum carboxylation rate (umol/m2/s)
%   flux.cp             ! CO2 compensation point (umol/mol)
%   flux.kc             ! Michaelis-Menten constant for CO2 (umol/mol)
%   flux.ko             ! Michaelis-Menten constant for O2 (mmol/mol)
%   flux.je             ! Electron transport rate (umol/m2/s)
%   flux.kp_c4          ! C4: Initial slope of CO2 response curve (mol/m2/s)
%   flux.gs             ! Leaf stomatal conductance (mol H2O/m2 leaf/s)
%   flux.gbc            ! Leaf boundary layer conductance, CO2 (mol CO2/m2 leaf/s)
%   flux.apar           ! Leaf absorbed PAR (umol photon/m2 leaf/s)
%   flux.rd             ! Leaf respiration rate (umol CO2/m2 leaf/s)
%
% Output
%   flux.ac             ! Leaf Rubisco-limited gross photosynthesis (umol CO2/m2 leaf/s)
%   flux.aj             ! Leaf RuBP regeneration-limited gross photosynthesis (umol CO2/m2 leaf/s)
%   flux.ap             ! Leaf product-limited (C3) or CO2-limited (C4) gross photosynthesis (umol CO2/m2 leaf/s)
%   flux.ag             ! Leaf gross photosynthesis (umol CO2/m2 leaf/s)
%   flux.an             ! Leaf net photosynthesis (umol CO2/m2 leaf/s)
%   flux.cs             ! Leaf surface CO2 (umol/mol)
%   flux.ci             ! Leaf intercellular CO2 (umol/mol)
% ------------------------------------------------------

% --- Leaf conductance
% gbc has units mol CO2/m2/s
% gs has units mol H2O/m2/s
% gleaf has units mol CO2/m2/s

gleaf = 1 / (1 / flux.gbc + 1.6 / flux.gs);

% --- Gross assimilation rates

if (leaf.c3psn == 1)

   % C3: Rubisco-limited photosynthesis

   a0 = flux.vcmax;
   e0 = 1;
   d0 = flux.kc * (1 + atmos.o2air / flux.ko);

   aquad = e0 / gleaf;
   bquad = -(e0 * atmos.co2air + d0) - (a0 - e0 * flux.rd) / gleaf;
   cquad = a0 * (atmos.co2air - flux.cp) - flux.rd * (e0 * atmos.co2air + d0);

   pcoeff = [aquad bquad cquad];
   proots = roots(pcoeff);
   flux.ac = min(proots(1), proots(2)) + flux.rd;
 
   % C3: RuBP regeneration-limited photosynthesis

   a0 = flux.je;
   e0 = 4;
   d0 = 8 * flux.cp;

   aquad = e0 / gleaf;
   bquad = -(e0 * atmos.co2air + d0) - (a0 - e0 * flux.rd) / gleaf;
   cquad = a0 * (atmos.co2air - flux.cp) - flux.rd * (e0 * atmos.co2air + d0);

   pcoeff = [aquad bquad cquad];
   proots = roots(pcoeff);
   flux.aj = min(proots(1), proots(2)) + flux.rd;

   % C3: Product-limited photosynthesis

   flux.ap = 0;

else

   % C4: Rubisco-limited photosynthesis

   flux.ac = flux.vcmax;
 
   % C4: RuBP-limited photosynthesis

   flux.aj = leaf.qe_c4 * flux.apar;

   % C4: PEP carboxylase-limited (CO2-limited)

   flux.ap = flux.kp_c4 * (atmos.co2air * gleaf + flux.rd) / (gleaf + flux.kp_c4);

end

% --- Net assimilation as the minimum or co-limited rate

if (leaf.colim == 1) % Use co-limitation

   % First co-limit Ac and Aj. Ai is the intermediate co-limited photosynthesis
   % rate found by solving the polynomial: aquad*Ai^2 + bquad*Ai + cquad = 0 for Ai.
   % Correct solution is the smallest of the two roots.

   if (leaf.c3psn == 1)
      aquad = leaf.colim_c3;
   else
      aquad = leaf.colim_c4a;
   end
   bquad = -(flux.ac + flux.aj);
   cquad = flux.ac * flux.aj;
   pcoeff = [aquad bquad cquad];
   proots = roots(pcoeff);
   ai = min(proots(1), proots(2));

   % Now co-limit again using Ap, but only for C4 plants. Solve the polynomial:
   % aquad*Ag^2 + bquad*Ag + cquad = 0 for Ag. Correct solution is the smallest
   % of the two roots. Ignore the product-limited rate For C3 plants.

   if (leaf.c3psn == 0)
      aquad = leaf.colim_c4b;
      bquad = -(ai + flux.ap);
      cquad = ai * flux.ap;
      pcoeff = [aquad bquad cquad];
      proots = roots(pcoeff);
      flux.ag = min(proots(1), proots(2));
   else
      flux.ag = ai;
   end

elseif (leaf.colim == 0) % No co-limitation

   if (leaf.c3psn == 1)
      flux.ag = min(flux.ac, flux.aj); % C3
   else
      flux.ag = min(flux.ac, flux.aj, flux.ap); % C4
   end

end

% Prevent photosynthesis from ever being negative

flux.ac = max(flux.ac, 0);
flux.aj = max(flux.aj, 0);
flux.ap = max(flux.ap, 0);
flux.ag = max(flux.ag, 0);

% Net CO2 uptake

flux.an = flux.ag - flux.rd;

% --- CO2 at leaf surface

flux.cs = atmos.co2air - flux.an / flux.gbc;
flux.cs = max(flux.cs, 1);

% --- Intercelluar CO2

flux.ci = atmos.co2air - flux.an / gleaf;

function [val] = latvap (tc, mmh2o)

% Latent heat of vaporization (J/mol) at temperature tc (degC)

val = 2501.6 - 2.3773 * tc; % Latent heat of vaporization (J/g)
val = val * 1000;           % Convert from J/g to J/kg
val = val * mmh2o;          % Convert from J/kg to J/mol

function [flux] = LeafBoundaryLayer (physcon, atmos, leaf, flux)

% Leaf boundary layer conductances

% -------------------------------------------------------------------------
% Input
%   physcon.grav      ! Gravitational acceleration (m/s2)
%   physcon.tfrz      ! Freezing point of water (K)
%   physcon.visc0     ! Kinematic viscosity at 0C and 1013.25 hPa (m2/s)
%   physcon.Dh0       ! Molecular diffusivity (heat) at 0C and 1013.25 hPa (m2/s)
%   physcon.Dv0       ! Molecular diffusivity (H2O) at 0C and 1013.25 hPa (m2/s)
%   physcon.Dc0       ! Molecular diffusivity (CO2) at 0C and 1013.25 hPa (m2/s)
%   atmos.patm        ! Atmospheric pressure (Pa)
%   atmos.rhomol      ! Molar density mol/m3)
%   atmos.wind        ! Wind speed (m/s)
%   atmos.tair        ! Air temperature (K)
%   leaf.dleaf        ! Leaf dimension (m)
%   flux.tleaf        ! Leaf temperature (K)
%
% Output
%   flux.gbh          ! Leaf boundary layer conductance, heat (mol/m2 leaf/s)
%   flux.gbv          ! Leaf boundary layer conductance, H2O (mol H2O/m2 leaf/s)
%   flux.gbc          ! Leaf boundary layer conductance, CO2 (mol CO2/m2 leaf/s)
% -------------------------------------------------------------------------

% --- Adjust diffusivity for temperature and pressure

fac = 101325 / atmos.patm * (atmos.tair / physcon.tfrz)^1.81;

visc = physcon.visc0 * fac; % Kinematic viscosity (m2/s)
Dh = physcon.Dh0 * fac;     % Molecular diffusivity, heat (m2/s)
Dv = physcon.Dv0 * fac;     % Molecular diffusivity, H2O (m2/s)
Dc = physcon.Dc0 * fac;     % Molecular diffusivity, CO2 (m2/s)

% --- Dimensionless numbers

Re = atmos.wind * leaf.dleaf / visc; % Reynolds number
Pr = visc / Dh;                      % Prandtl number
Scv = visc / Dv;                     % Schmidt number for H2O
Scc = visc / Dc;                     % Schmidt number for CO2

% Grashof number

Gr = physcon.grav * leaf.dleaf^3 * max(flux.tleaf-atmos.tair, 0) / (atmos.tair * visc * visc);

% --- Empirical correction factor for Nu and Sh

b1 = 1.5;

% --- Nusselt number (Nu) and Sherwood numbers (H2O: Shv, CO2: Shc)

% Forced convection - laminar flow

Nu_lam  = b1 * 0.66 *  Pr^0.33 * Re^0.5;     % Nusselt number
Shv_lam = b1 * 0.66 * Scv^0.33 * Re^0.5;     % Sherwood number, H2O
Shc_lam = b1 * 0.66 * Scc^0.33 * Re^0.5;     % Sherwood number, CO2

% Forced convection - turbulent flow

Nu_turb  = b1 * 0.036 *  Pr^0.33 * Re^0.8;   % Nusselt number
Shv_turb = b1 * 0.036 * Scv^0.33 * Re^0.8;   % Sherwood number, H2O
Shc_turb = b1 * 0.036 * Scc^0.33 * Re^0.8;   % Sherwood number, CO2

% Choose correct flow regime for forced convection

Nu_forced = max(Nu_lam, Nu_turb);
Shv_forced = max(Shv_lam, Shv_turb);
Shc_forced = max(Shc_lam, Shc_turb);

% Free convection

Nu_free  = 0.54 *  Pr^0.25 * Gr^0.25;        % Nusselt number
Shv_free = 0.54 * Scv^0.25 * Gr^0.25;        % Sherwood number, H2O
Shc_free = 0.54 * Scc^0.25 * Gr^0.25;        % Sherwood number, CO2

% Both forced and free convection regimes occur together

Nu = Nu_forced + Nu_free;
Shv = Shv_forced + Shv_free;
Shc = Shc_forced + Shc_free;

% --- Boundary layer conductances (m/s)

flux.gbh = Dh *  Nu / leaf.dleaf;
flux.gbv = Dv * Shv / leaf.dleaf;
flux.gbc = Dc * Shc / leaf.dleaf;

% --- Convert conductance (m/s) to (mol/m2/s)

flux.gbh = flux.gbh * atmos.rhomol;
flux.gbv = flux.gbv * atmos.rhomol;
flux.gbc = flux.gbc * atmos.rhomol;

function [flux] = LeafFluxes (physcon, atmos, leaf, flux)

% --- Leaf temperature, energy fluxes, photosynthesis, and stomatal conductance
% using stomatal optimization

[flux] = StomataOptimization (physcon, atmos, leaf, flux);

function [flux] = LeafPhotosynthesis (physcon, atmos, leaf, flux)

% Calculate leaf photosynthesis for a specified stomatal
% conductance. Then calculate Ci from the diffusion equation.
% This is used with the WUE stomatal optimization.

% ------------------------------------------------------
% Input
%   physcon.tfrz        ! Freezing point of water (K)
%   physcon.rgas        ! Universal gas constant (J/K/mol)
%   atmos.co2air        ! Atmospheric CO2 (umol/mol)
%   atmos.eair          ! Vapor pressure of air (Pa)
%   leaf.c3psn          ! Photosynthetic pathway: 1 = C3. 0 = C4 plant
%   leaf.vcmax25        ! Maximum carboxylation rate at 25C (umol/m2/s)
%   leaf.jmax25         ! Maximum electron transport rate at 25C (umol/m2/s)
%   leaf.rd25           ! Leaf respiration rate at 25C (umol CO2/m2/s)
%   leaf.kc25           ! Michaelis-Menten constant for CO2 at 25C (umol/mol)
%   leaf.ko25           ! Michaelis-Menten constant for O2 at 25C (mmol/mol)
%   leaf.cp25           ! CO2 compensation point at 25C (umol/mol)
%   leaf.kcha           ! Activation energy for Kc (J/mol)
%   leaf.koha           ! Activation energy for Ko (J/mol)
%   leaf.cpha           ! Activation energy for Cp (J/mol)
%   leaf.vcmaxha        ! Activation energy for Vcmax (J/mol)
%   leaf.jmaxha         ! Activation energy for Jmax (J/mol)
%   leaf.rdha           ! Activation energy for Rd (J/mol)
%   leaf.vcmaxhd        ! Deactivation energy for Vcmax (J/mol)
%   leaf.jmaxhd         ! Deactivation energy for Jmax (J/mol)
%   leaf.rdhd           ! Deactivation energy for Rd (J/mol)
%   leaf.vcmaxse        ! Entropy term for Vcmax (J/mol/K)
%   leaf.jmaxse         ! Entropy term for Jmax (J/mol/K)
%   leaf.rdse           ! Entropy term for Rd (J/mol/K)
%   leaf.vcmaxc         ! Vcmax scaling factor for high temperature inhibition (25 C = 1.0)
%   leaf.jmaxc          ! Jmax scaling factor for high temperature inhibition (25 C = 1.0)
%   leaf.rdc            ! Rd scaling factor for high temperature inhibition (25 C = 1.0)
%   leaf.phi_psii       ! Quantum yield of PS II
%   leaf.theta_j        ! Empirical curvature parameter for electron transport rate
%   leaf.kp25_c4        ! C4: Initial slope of CO2 response curve at 25C (mol/m2/s)
%   flux.gbv            ! Leaf boundary layer conductance, H2O (mol H2O/m2 leaf/s)
%   flux.gbc            ! Leaf boundary layer conductance, CO2 (mol CO2/m2 leaf/s)
%   flux.apar           ! Leaf absorbed PAR (umol photon/m2 leaf/s)
%   flux.tleaf          ! Leaf temperature (K)
%   flux.gs             ! Leaf stomatal conductance (mol H2O/m2 leaf/s)
%
% Output
%   flux.vcmax          ! Maximum carboxylation rate (umol/m2/s)
%   flux.jmax           ! Maximum electron transport rate (umol/m2/s)
%   flux.cp             ! CO2 compensation point (umol/mol)
%   flux.kc             ! Michaelis-Menten constant for CO2 (umol/mol)
%   flux.ko             ! Michaelis-Menten constant for O2 (mmol/mol)
%   flux.je             ! Electron transport rate (umol/m2/s)
%   flux.kp_c4          ! C4: Initial slope of CO2 response curve (mol/m2/s)
%   flux.ac             ! Leaf Rubisco-limited gross photosynthesis (umol CO2/m2 leaf/s)
%   flux.aj             ! Leaf RuBP regeneration-limited gross photosynthesis (umol CO2/m2 leaf/s)
%   flux.ap             ! Leaf product-limited (C3) or CO2-limited (C4) gross photosynthesis (umol CO2/m2 leaf/s)
%   flux.ag             ! Leaf gross photosynthesis (umol CO2/m2 leaf/s)
%   flux.an             ! Leaf net photosynthesis (umol CO2/m2 leaf/s)
%   flux.rd             ! Leaf respiration rate (umol CO2/m2 leaf/s)
%   flux.cs             ! Leaf surface CO2 (umol/mol)
%   flux.ci             ! Leaf intercellular CO2 (umol/mol)
%   flux.hs             ! Leaf fractional humidity at surface (dimensionless)
%   flux.vpd            ! Leaf vapor pressure deficit at surface (Pa)
% ------------------------------------------------------

% --- Adjust photosynthetic parameters for temperature

if (leaf.c3psn == 1)

   % C3 temperature response

   ft = @(tl, ha) exp(ha/(physcon.rgas*(physcon.tfrz+25)) * (1-(physcon.tfrz+25)/tl));
   fth = @(tl, hd, se, fc) fc / (1 + exp((-hd+se*tl)/(physcon.rgas*tl)));

   flux.kc = leaf.kc25 * ft(flux.tleaf, leaf.kcha);
   flux.ko = leaf.ko25 * ft(flux.tleaf, leaf.koha);
   flux.cp = leaf.cp25 * ft(flux.tleaf, leaf.cpha);

   t1 = ft(flux.tleaf, leaf.vcmaxha);
   t2 = fth(flux.tleaf, leaf.vcmaxhd, leaf.vcmaxse, leaf.vcmaxc);
   flux.vcmax = leaf.vcmax25 * t1 * t2;

   t1 = ft(flux.tleaf, leaf.jmaxha);
   t2 = fth(flux.tleaf, leaf.jmaxhd, leaf.jmaxse, leaf.jmaxc);
   flux.jmax = leaf.jmax25 * t1 * t2;

   t1 = ft(flux.tleaf, leaf.rdha);
   t2 = fth(flux.tleaf, leaf.rdhd, leaf.rdse, leaf.rdc);
   flux.rd = leaf.rd25 * t1 * t2;

   flux.kp_c4 = 0;

elseif (leaf.c3psn == 0)

   % C4 temperature response

   t1 = 2^((flux.tleaf-(physcon.tfrz+25)) / 10);
   t2 = 1 + exp(0.2*((physcon.tfrz+15)-flux.tleaf));
   t3 = 1 + exp(0.3*(flux.tleaf-(physcon.tfrz+40)));
   flux.vcmax = leaf.vcmax25 * t1 / (t2 * t3);

   t3 = 1 + exp(1.3*(flux.tleaf-(physcon.tfrz+55)));
   flux.rd = leaf.rd25  * t1 / t3;

   flux.kp_c4 = leaf.kp25_c4 * t1;

   flux.kc = 0;
   flux.ko = 0;
   flux.cp = 0;
   flux.jmax = 0;
   flux.je = 0;

end

% --- Electron transport rate for C3 plants

% Solve the polynomial: aquad*Je^2 + bquad*Je + cquad = 0
% for Je. Correct solution is the smallest of the two roots.

if (leaf.c3psn == 1)
   qabs = 0.5 * leaf.phi_psii * flux.apar;
   aquad = leaf.theta_j;
   bquad = -(qabs + flux.jmax);
   cquad = qabs * flux.jmax;
   pcoeff = [aquad bquad cquad];
   proots = roots(pcoeff);
   flux.je = min(proots(1), proots(2));
end

% --- Ci calculation. Calculate photosynthesis for a specified stomatal conductance

[flux] = CiFuncOptimization (atmos, leaf, flux);

% --- Relative humidity and vapor pressure at leaf surface

[esat, desat] = satvap ((flux.tleaf-physcon.tfrz));
flux.hs = (flux.gbv * atmos.eair + flux.gs * esat) / ((flux.gbv + flux.gs) * esat);
flux.vpd = max(esat - flux.hs*esat, 0.1);

% Compare with diffusion equation: An = (ca - ci) * gleaf

an_err = (atmos.co2air - flux.ci) / (1 / flux.gbc + 1.6 / flux.gs);
if (flux.an > 0 & abs(flux.an-an_err) > 0.01)
   fprintf('An = %15.4f\n', flux.an)
   fprintf('An_err = %15.4f\n', an_err)
   error ('LeafPhotosynthesis: failed diffusion error check')
end

function [leaf] = LeafPhysiologyParams (params, physcon, leaf)

% ------------------------------------------------------
% Input
%   params.vis          ! Waveband index for visible radiation
%   params.nir          ! Waveband index for near-infrared radiation
%   physcon.tfrz        ! Freezing point of water (K)
%   physcon.rgas        ! Universal gas constant (J/K/mol)
%   leaf.c3psn          ! Photosynthetic pathway: 1 = C3. 0 = C4 plant
%
% Output
%   leaf.vcmax25        ! Maximum carboxylation rate at 25C (umol/m2/s)
%   leaf.jmax25         ! Maximum electron transport rate at 25C (umol/m2/s)
%   leaf.rd25           ! Leaf respiration rate at 25C (umol CO2/m2/s)
%   leaf.kc25           ! Michaelis-Menten constant for CO2 at 25C (umol/mol)
%   leaf.ko25           ! Michaelis-Menten constant for O2 at 25C (mmol/mol)
%   leaf.cp25           ! CO2 compensation point at 25C (umol/mol)
%   leaf.kcha           ! Activation energy for Kc (J/mol)
%   leaf.koha           ! Activation energy for Ko (J/mol)
%   leaf.cpha           ! Activation energy for Cp (J/mol)
%   leaf.vcmaxha        ! Activation energy for Vcmax (J/mol)
%   leaf.jmaxha         ! Activation energy for Jmax (J/mol)
%   leaf.rdha           ! Activation energy for Rd (J/mol)
%   leaf.vcmaxhd        ! Deactivation energy for Vcmax (J/mol)
%   leaf.jmaxhd         ! Deactivation energy for Jmax (J/mol)
%   leaf.rdhd           ! Deactivation energy for Rd (J/mol)
%   leaf.vcmaxse        ! Entropy term for Vcmax (J/mol/K)
%   leaf.jmaxse         ! Entropy term for Jmax (J/mol/K)
%   leaf.rdse           ! Entropy term for Rd (J/mol/K)
%   leaf.vcmaxc         ! Vcmax scaling factor for high temperature inhibition (25 C = 1.0)
%   leaf.jmaxc          ! Jmax scaling factor for high temperature inhibition (25 C = 1.0)
%   leaf.rdc            ! Rd scaling factor for high temperature inhibition (25 C = 1.0)
%   leaf.phi_psii       ! Quantum yield of PS II
%   leaf.theta_j        ! Empirical curvature parameter for electron transport rate
%   leaf.colim_c3       ! Empirical curvature parameter for C3 co-limitation
%   leaf.colim_c4a      ! Empirical curvature parameter for C4 co-limitation
%   leaf.colim_c4b      ! Empirical curvature parameter for C4 co-limitation
%   leaf.qe_c4          ! C4: Quantum yield (mol CO2 / mol photons)
%   leaf.kp25_c4        ! C4: Initial slope of CO2 response curve at 25C (mol/m2/s)
%   leaf.dleaf          ! Leaf dimension (m)
%   leaf.emiss          ! Leaf emissivity
%   leaf.rho            ! Leaf reflectance for visible (vis) and near-infrared (nir) wavebands
%   leaf.tau            ! Leaf transmittance for visible (vis) and near-infrared (nir) wavebands
%   leaf.iota           ! Stomatal efficiency (umol CO2/ mol H2O)
%   leaf.capac          ! Plant capacitance (mmol H2O/m2 leaf area/MPa)
%   leaf.minlwp         ! Minimum leaf water potential (MPa)
%   leaf.gplant         ! Stem (xylem-to-leaf) hydraulic conductance (mmol H2O/m2 leaf area/s/Mpa)
% ------------------------------------------------------

% --- Vcmax and other parameters (at 25C)

if (leaf.c3psn == 1)
   leaf.vcmax25 = 60;
   leaf.jmax25 = 1.67 * leaf.vcmax25;
   leaf.kp25_c4 = 0;
   leaf.rd25 = 0.015 * leaf.vcmax25;
else
   leaf.vcmax25 = 40;
   leaf.jmax25 = 0;
   leaf.kp25_c4 = 0.02 * leaf.vcmax25;
   leaf.rd25 = 0.025 * leaf.vcmax25;
end

% --- Kc, Ko, Cp at 25C

leaf.kc25 = 404.9;
leaf.ko25 = 278.4;
leaf.cp25 = 42.75;

% --- Activation energy

leaf.kcha = 79430;
leaf.koha = 36380;
leaf.cpha = 37830;
leaf.rdha = 46390;
leaf.vcmaxha = 65330;
leaf.jmaxha  = 43540;

% --- High temperature deactivation

% Deactivation energy (J/mol)

leaf.rdhd = 150000;
leaf.vcmaxhd = 150000;
leaf.jmaxhd  = 150000;

% Entropy term (J/mol/K)

leaf.rdse = 490;
leaf.vcmaxse = 490;
leaf.jmaxse  = 490;

% Scaling factors for high temperature inhibition (25 C = 1.0).
% The factor "c" scales the deactivation to a value of 1.0 at 25C.

fth25 = @(hd, se) 1 + exp((-hd + se*(physcon.tfrz+25)) / (physcon.rgas*(physcon.tfrz+25)));
leaf.vcmaxc = fth25 (leaf.vcmaxhd, leaf.vcmaxse);
leaf.jmaxc  = fth25 (leaf.jmaxhd, leaf.jmaxse);
leaf.rdc    = fth25 (leaf.rdhd, leaf.rdse);

% --- C3 parameters

% Quantum yield of PS II

leaf.phi_psii = 0.85;

% Empirical curvature parameter for electron transport rate

leaf.theta_j = 0.90;

% Empirical curvature parameter for C3 co-limitation

leaf.colim_c3 = 0.98;

% Empirical curvature parameters for C4 co-limitation

leaf.colim_c4a = 0.80;
leaf.colim_c4b = 0.95;

% --- C4: Quantum yield (mol CO2 / mol photons)

leaf.qe_c4 = 0.05;

% --- Stomatal efficiency for optimization (An/E; umol CO2/ mol H2O)

leaf.iota = 750;

% --- Leaf dimension (m)

leaf.dleaf = 0.05;

% --- Leaf emissivity

leaf.emiss = 0.98;

% --- Leaf reflectance and transmittance: visible and near-infrared wavebands

leaf.rho(params.vis) = 0.10;
leaf.tau(params.vis) = 0.10;
leaf.rho(params.nir) = 0.40;
leaf.tau(params.nir) = 0.40;

% --- Plant hydraulic parameters

% Plant capacitance (mmol H2O/m2 leaf area/MPa)

leaf.capac = 2500;

% Minimum leaf water potential (MPa)

leaf.minlwp = -2;

% Stem (xylem-to-leaf) hydraulic conductance (mmol H2O/m2 leaf area/s/Mpa)

leaf.gplant = 4;

function [flux] = LeafTemperature (physcon, atmos, leaf, flux)

% Leaf temperature and energy fluxes

% ------------------------------------------------------
% Input
%   physcon.tfrz     ! Freezing point of water (K)
%   physcon.mmh2o    ! Molecular mass of water (kg/mol)
%   physcon.sigma    ! Stefan-Boltzmann constant (W/m2/K4)
%   atmos.patm       ! Atmospheric pressure (Pa)
%   atmos.cpair      ! Specific heat of air at constant pressure (J/mol/K)
%   atmos.tair       ! Air temperature (K)
%   atmos.eair       ! Vapor pressure of air (Pa)
%   leaf.emiss       ! Leaf emissivity
%   flux.gbh         ! Leaf boundary layer conductance, heat (mol/m2 leaf/s)
%   flux.gbv         ! Leaf boundary layer conductance, H2O (mol H2O/m2 leaf/s)
%   flux.gs          ! Leaf stomatal conductance (mol H2O/m2 leaf/s)
%   flux.qa          ! Leaf radiative forcing (W/m2 leaf)
%
% Input/ouput
%   flux.tleaf       ! Leaf temperature (K)
%
% Output
%   flux.rnet        ! Leaf net radiation (W/m2 leaf)
%   flux.lwrad       ! Longwave radiation emitted from leaf (W/m2 leaf)
%   flux.shflx       ! Leaf sensible heat flux (W/m2 leaf)
%   flux.lhflx       ! Leaf latent heat flux (W/m2 leaf)
%   flux.etflx       ! Leaf transpiration flux (mol H2O/m2 leaf/s)
% ------------------------------------------------------

% --- Latent heat of vaporization (J/mol)

[lambda_val] = latvap ((atmos.tair-physcon.tfrz), physcon.mmh2o);

% --- Newton-Raphson iteration until leaf energy balance is less than
% f0_max or to niter_max iterations

niter = 0;         % Number of iterations
f0 = 1e36;         % Leaf energy balance (W/m2)

niter_max = 100;   % Maximum number of iterations
f0_max = 1e-06;    % Maximum energy imbalance (W/m2)

while (niter <= niter_max & abs(f0) > f0_max)

   % Increment iteration counter

   niter = niter + 1;

   % Saturation vapor pressure (Pa) and temperature derivative (Pa/K)

   [esat, desat] = satvap ((flux.tleaf-physcon.tfrz));

   % Leaf conductance for water vapor (mol H2O/m2/s)

   gleaf = flux.gs * flux.gbv / (flux.gs + flux.gbv);

   % Emitted longwave radiation (W/m2) and temperature derivative (W/m2/K)

   flux.lwrad = 2 * leaf.emiss * physcon.sigma * flux.tleaf^4;
   dlwrad = 8 * leaf.emiss * physcon.sigma * flux.tleaf^3;

   % Sensible heat flux (W/m2) and temperature derivative (W/m2/K)

   flux.shflx = 2 * atmos.cpair * (flux.tleaf - atmos.tair) * flux.gbh;
   dshflx = 2 * atmos.cpair * flux.gbh;

   % Latent heat flux (W/m2) and temperature derivative (W/m2/K)

   flux.lhflx = lambda_val / atmos.patm * (esat - atmos.eair) * gleaf;
   dlhflx = lambda_val / atmos.patm * desat * gleaf;

   % Energy balance (W/m2) and temperature derivative (W/m2/K)

   f0 = flux.qa - flux.lwrad - flux.shflx - flux.lhflx;
   df0 = -dlwrad - dshflx - dlhflx;

   % Change in leaf temperature

   dtleaf = -f0 / df0;

   % Update leaf temperature

   flux.tleaf = flux.tleaf + dtleaf;

end

% --- Net radiation

flux.rnet = flux.qa - flux.lwrad;

% --- Error check

err = flux.rnet - flux.shflx - flux.lhflx;
if (abs(err) > f0_max)
   fprintf('err = %15.3f\n',err)
   fprintf('qa = %15.3f\n',flux.qa)
   fprintf('lwrad = %15.3f\n',flux.lwrad)
   fprintf('sh = %15.3f\n',flux.shflx)
   fprintf('lh = %15.3f\n',flux.lhflx)
   error ('LeafTemperature error')
end

% Water vapor flux: W/m2 -> mol H2O/m2/s

flux.etflx = flux.lhflx / lambda_val;

function [leafwp] = LeafWaterPotential (physcon, leaf, flux)

% Calculate leaf water potential for the current transpiration rate

% -------------------------------------------------------------------------
% Input
%   physcon.grav     ! Gravitational acceleration (m/s2)
%   physcon.denh2o   ! Density of liquid water (kg/m3)
%   leaf.capac       ! Plant capacitance (mmol H2O/m2 leaf area/MPa)
%   flux.height      ! Leaf height (m)
%   flux.lsc         ! Leaf-specific conductance (mmol H2O/m2 leaf/s/MPa)
%   flux.psi_leaf    ! Leaf water potential at beginning of time step (MPa)
%   flux.psi_soil    ! Weighted soil water potential (MPa)
%   flux.etflx       ! Leaf transpiration flux (mol H2O/m2 leaf/s)
%   flux.dt          ! Model time step (s)
% Output
%   leafwp           ! Leaf water potential at end of time step (MPa)
% -------------------------------------------------------------------------

% --- Head of pressure (MPa/m)

head = physcon.denh2o * physcon.grav * 1e-06;

% --- Change in leaf water potential is: dy / dt = (a - y) / b. The integrated change 
% over a full model timestep is: dy = (a - y0) * (1 - exp(-dt/b))

y0 = flux.psi_leaf;
a = flux.psi_soil - head * flux.height - 1000 * flux.etflx / flux.lsc;
b = leaf.capac / flux.lsc;
dy = (a - y0) * (1 - exp(-flux.dt/b));
leafwp = y0 + dy;

function [rootvar] = RootParams

% Fine root radius (m)

rootvar.radius = 0.29e-03;

% Fine root density (g biomass / m3 root)

rootvar.density = 0.31e06;

% Hydraulic resistivity of root tissue (MPa.s.g/mmol H2O)

rootvar.resist = 25;

function [esat, desat] = satvap (tc)

% Compute saturation vapor pressure and change in saturation vapor pressure
% with respect to temperature. Polynomial approximations are from:
% Flatau et al. (1992) Polynomial fits to saturation vapor pressure.
% Journal of Applied Meteorology 31:1507-1513. Input temperature is Celsius.

% --- For water vapor (temperature range is 0C to 100C)
 
a0 =  6.11213476;        b0 =  0.444017302;
a1 =  0.444007856;       b1 =  0.286064092e-01;
a2 =  0.143064234e-01;   b2 =  0.794683137e-03;
a3 =  0.264461437e-03;   b3 =  0.121211669e-04;
a4 =  0.305903558e-05;   b4 =  0.103354611e-06;
a5 =  0.196237241e-07;   b5 =  0.404125005e-09;
a6 =  0.892344772e-10;   b6 = -0.788037859e-12;
a7 = -0.373208410e-12;   b7 = -0.114596802e-13;
a8 =  0.209339997e-15;   b8 =  0.381294516e-16;
 
% --- For ice (temperature range is -75C to 0C)
 
c0 =  6.11123516;        d0 =  0.503277922;
c1 =  0.503109514;       d1 =  0.377289173e-01;
c2 =  0.188369801e-01;   d2 =  0.126801703e-02;
c3 =  0.420547422e-03;   d3 =  0.249468427e-04;
c4 =  0.614396778e-05;   d4 =  0.313703411e-06;
c5 =  0.602780717e-07;   d5 =  0.257180651e-08;
c6 =  0.387940929e-09;   d6 =  0.133268878e-10;
c7 =  0.149436277e-11;   d7 =  0.394116744e-13;
c8 =  0.262655803e-14;   d8 =  0.498070196e-16;

% --- Limit temperature to -75C to 100C
 
tc = min(tc, 100);
tc = max(tc, -75);

% --- Saturation vapor pressure (esat, mb) and derivative (desat, mb)

if (tc >= 0)
   esat  = a0 + tc*(a1 + tc*(a2 + tc*(a3 + tc*(a4 ...
         + tc*(a5 + tc*(a6 + tc*(a7 + tc*a8)))))));
   desat = b0 + tc*(b1 + tc*(b2 + tc*(b3 + tc*(b4 ...
         + tc*(b5 + tc*(b6 + tc*(b7 + tc*b8)))))));
else
   esat  = c0 + tc*(c1 + tc*(c2 + tc*(c3 + tc*(c4 ...
         + tc*(c5 + tc*(c6 + tc*(c7 + tc*c8)))))));
   desat = d0 + tc*(d1 + tc*(d2 + tc*(d3 + tc*(d4 ...
         + tc*(d5 + tc*(d6 + tc*(d7 + tc*d8)))))));
end

% --- Convert from mb to Pa

esat  = esat  * 100;
desat = desat * 100;

function [soil] = SoilParams (soil)

% Set soil hydraulic parameters, soil depth, and rooting fraction

% -------------------------------------------------------------------------
% Input
%   soil.texture     ! Soil texture class
% Output
%   soil.nlevsoi     ! Number of soil layers
%   soil.dz          ! Soil layer thickness (m)
%   soil.rootfr      ! Fraction of roots in each soil layer (-)
%   soil.watsat      ! Soil layer volumetric water content at saturation (porosity)
%   soil.psisat      ! Soil layer matric potential at saturation (mm)
%   soil.hksat       ! Soil layer hydraulic conductivity at saturation (mm H2O/s)
%   soil.bsw         ! Soil layer Clapp and Hornberger "b" parameter
% -------------------------------------------------------------------------

% --- Soil texture classes
%  1: sand
%  2: loamy sand
%  3: sandy loam
%  4: silt loam
%  5: loam
%  6: sandy clay loam
%  7  silty clay loam
%  8: clay loam
%  9: sandy clay
% 10: silty clay
% 11: clay

% Volumetric water content at saturation (porosity)

theta_sat = [0.395, 0.410, 0.435, 0.485, 0.451, 0.420, 0.477, 0.476, 0.426, 0.492, 0.482];

% Matric potential at saturation (mm)

psi_sat = [-121, -90, -218, -786, -478, -299, -356, -630, -153, -490, -405];

% Clapp and Hornberger "b" parameter

b = [4.05, 4.38, 4.90, 5.30, 5.39, 7.12, 7.75, 8.52, 10.40, 10.40, 11.40];

% Hydraulic conductivity at saturation (cm/min -> mm/s)

k_sat = [1.056, 0.938, 0.208, 0.0432, 0.0417, 0.0378, 0.0102, 0.0147, 0.0130, 0.0062, 0.0077];
k_sat = k_sat * 10 / 60;

% --- Soil layer thickness (total depth is 2.5 m)

soil.nlevsoi = 11;
soil.dz = [0.05, 0.05, 0.10, 0.10, 0.20, 0.20, 0.20, 0.30, 0.40, 0.40, 0.50];

% --- Root profile

  beta_param = 0.90;  % root profile parameter: shallow profile
% beta_param = 0.97;  % root profile parameter: deep profile

for j = 1:soil.nlevsoi
   if (j == 1)
      z2 = soil.dz(j) * 100;
      soil.rootfr(j) = 1 - beta_param^z2;
   else
      z1 = z2;
      z2 = z1 + soil.dz(j) * 100;
      soil.rootfr(j) = beta_param^z1 - beta_param^z2;
   end
end

% --- Soil texture to process

k = soil.texture;

% --- Soil layer values

for j = 1:soil.nlevsoi
   soil.watsat(j) = theta_sat(k);
   soil.psisat(j) = psi_sat(k);
   soil.hksat(j) = k_sat(k);
   soil.bsw(j) = b(k);
end

function [flux] = SoilResistance (physcon, leaf, rootvar, soil, flux)

% Calculate soil hydraulic resistance, weighted soil water potential and
% water uptake from each soil layer

% ----------------------------------------------------------------------------------
% Input
%   physcon.grav      ! Gravitational acceleration (m/s2)
%   physcon.denh2o    ! Density of liquid water (kg/m3)
%   physcon.mmh2o     ! Molecular mass of water (kg/mol)
%   leaf.minlwp       ! Minimum leaf water potential (MPa)
%   rootvar.biomass   ! Fine root biomass (g biomass / m2)
%   rootvar.radius    ! Fine root radius (m)
%   rootvar.density   ! Fine root density (g biomass / m3 root)
%   rootvar.resist    ! Hydraulic resistivity of root tissue (MPa.s.g/mmol H2O)
%   soil.nlevsoi      ! Number of soil layers
%   soil.h2osoi_vol   ! Soil layer volumetric water content (m3/m3)
%   soil.psi          ! Soil layer matric potential (mm)
%   soil.watsat       ! Soil layer volumetric water content at saturation (porosity)
%   soil.hksat        ! Soil layer hydraulic conductivity at saturation (mm H2O/s)
%   soil.bsw          ! Soil layer Clapp and Hornberger "b" parameter
%   soil.rootfr       ! Fraction of roots in each soil layer (-)
%   soil.dz           ! Soil layer thickness (m)
%   flux.lai          ! Canopy leaf area index (m2/m2)
% Output
%   flux.rsoil        ! Soil hydraulic resistance (MPa.s.m2/mmol H2O)
%   flux.psi_soil     ! Weighted soil water potential (MPa)
%   flux.et_loss      ! Fraction of total transpiration from each soil layer (-)
% ----------------------------------------------------------------------------------

% --- Head of pressure  (MPa/m)

head = physcon.denh2o * physcon.grav * 1e-06;

% --- Root cross-sectional area (m2 root)

root_cross_sec_area = pi * rootvar.radius^2;

% --- Soil and root resistances for each layer

flux.rsoil = 0;
for j = 1:soil.nlevsoi

   % Hydraulic conductivity for each layer (mmol/m/s/MPa)

   s = max(min(soil.h2osoi_vol(j)/soil.watsat(j), 1), 0.01);
   hk = soil.hksat(j) * s^(2 * soil.bsw(j) + 3);     % mm/s
   hk = hk * 1e-03 / head;                           % mm/s -> m/s -> m2/s/MPa
   hk = hk * physcon.denh2o / physcon.mmh2o * 1000;  % m2/s/MPa -> mmol/m/s/MPa

   % Matric potential for each layer (MPa)

   psi_mpa(j) = soil.psi(j) * 1e-03 * head;          % mm -> m -> MPa

   % Root biomass density (g biomass / m3 soil)

   root_biomass_density = rootvar.biomass * soil.rootfr(j) / soil.dz(j);
   root_biomass_density = max(root_biomass_density, 1e-10);

   % Root length density (m root per m3 soil)

   root_length_density = root_biomass_density / (rootvar.density * root_cross_sec_area);

   % Distance between roots (m)

   root_dist = sqrt(1 / (root_length_density * pi));

   % Soil-to-root resistance (MPa.s.m2/mmol H2O)

   soilr1 = log(root_dist/rootvar.radius) / (2 * pi * root_length_density * soil.dz(j) * hk);

   % Root-to-stem resistance (MPa.s.m2/mmol H2O)

   soilr2 = rootvar.resist / (root_biomass_density * soil.dz(j));

   % Belowground resistance (MPa.s.m2/mmol H2O) 

   soilr = soilr1 + soilr2;

   % Total belowground resistance. First sum the conductances (1/soilr)
   % for each soil layer and then convert back to a resistance after the
   % summation

   flux.rsoil = flux.rsoil + 1 / soilr;

   % Maximum transpiration for each layer (mmol H2O/m2/s)

   evap(j) = (psi_mpa(j) - leaf.minlwp) / soilr;
   evap(j) = max (evap(j), 0);

end

% Belowground resistance: resistance = 1 / conductance

flux.rsoil = flux.lai / flux.rsoil;

% Weighted soil water potential (MPa)

flux.psi_soil = 0;
for j = 1:soil.nlevsoi
   flux.psi_soil = flux.psi_soil + psi_mpa(j) * evap(j);
end

totevap = sum(evap);      % Total maximum transpiration
if (totevap > 0)
   flux.psi_soil = flux.psi_soil / totevap;
else
   flux.psi_soil = leaf.minlwp;
end

% Fractional transpiration uptake from soil layers

for j = 1:soil.nlevsoi
   if (totevap > 0)
      flux.et_loss(j) = evap(j) / totevap;
   else
      flux.et_loss(j) = 1 / soil.nlevsoi;
   end
end

function [flux, val] = StomataEfficiency (physcon, atmos, leaf, flux, gs_val)

% Stomata water-use efficiency check and cavitation check to determine maximum gs. 
% For the stomatal conductance gs_val, calculate photosynthesis and leaf
% water potential for an increase in stomatal conductance equal to "delta".
% The returned value is positive if this increase produces a change in
% photosynthesis > iota*vpd*delta or if the leaf water potential is > minlwp.
% The returned value is negative if the increase produces a change in
% photosynthesis < iota*vpd*delta or if the leaf water potential is < minlwp.

% --- Leaf boundary layer conductances

[flux] = LeafBoundaryLayer (physcon, atmos, leaf, flux);

% --- Specify "delta" as a small difference in gs (mol H2O/m2/s)

delta = 0.001;

% --- Calculate photosynthesis at lower gs (gs_val - delta), but first need leaf temperature for this gs
% gs2 - lower value for gs (mol H2O/m2/s)
% an2 - leaf photosynthesis at gs2 (umol CO2/m2/s)

gs2 = gs_val - delta;
flux.gs = gs2;
[flux] = LeafTemperature (physcon, atmos, leaf, flux);
[flux] = LeafPhotosynthesis (physcon, atmos, leaf, flux);
an2 = flux.an;

% --- Calculate photosynthesis at higher gs (gs_val), but first need leaf temperature for this gs
% gs1 - higher value for gs (mol H2O/m2/s)
% an1 - leaf photosynthesis at gs1 (umol CO2/m2/s)

gs1 = gs_val;
flux.gs = gs1;
[flux] = LeafTemperature (physcon, atmos, leaf, flux);
[flux] = LeafPhotosynthesis (physcon, atmos, leaf, flux);
an1 = flux.an;

% --- Efficiency check: wue < 0 when d(An) / d(gs) < iota * vpd

wue = (an1 - an2) - leaf.iota * delta * (flux.vpd / atmos.patm);

% --- Cavitation check: minpsi < 0 when leafwp < minlwp

[leafwp] = LeafWaterPotential (physcon, leaf, flux);
minpsi = leafwp - leaf.minlwp;

% --- Return the minimum of the two checks

val = min(wue, minpsi);

function [flux] = StomataOptimization (physcon, atmos, leaf, flux)

% Leaf temperature, energy fluxes, photosynthesis, and stomatal conductance
% with water-use efficiency stomatal optimization

% --- Low and high initial estimates for gs (mol H2O/m2/s)

gs1 = 0.002;
gs2 = 2.0;

% --- Check for minimum stomatal conductance linked to low light or drought stress based
% on the water-use efficiency and cavitation checks for gs1 and gs2 (check1, check2)

[flux, check1] = StomataEfficiency (physcon, atmos, leaf, flux, gs1);
[flux, check2] = StomataEfficiency (physcon, atmos, leaf, flux, gs2);

if (check1 * check2 < 0)

   % Calculate gs using the function StomataEfficiency to iterate gs
   % to an accuracy of tol (mol H2O/m2/s)

   tol = 0.004;
   func_name = 'StomataEfficiency';
   [flux, root] = brent_root (func_name, physcon, atmos, leaf, flux, gs1, gs2, tol);
   flux.gs = root;

else

   % Low light or drought stress. Set gs to minimum conductance

   flux.gs = 0.002;

end

% --- Leaf fluxes for this gs

[flux] = LeafBoundaryLayer (physcon, atmos, leaf, flux);
[flux] = LeafTemperature (physcon, atmos, leaf, flux);
[flux] = LeafPhotosynthesis (physcon, atmos, leaf, flux);
[flux.psi_leaf] = LeafWaterPotential (physcon, leaf, flux);