Implicit Flux–Profile Solution

function [fluxvar, root] = brent_root (func, physcon, forcvar, surfvar, fluxvar, 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 [fluxvar, fx] = func (physcon, forcvar, surfvar, fluxvar, x)
%
% The function func is exaluated at x and the returned value is fx. It uses variables
% in the physcon, forcvar, surfvar, and fluxvar structures. These are passed in as
% input arguments. It also calculates values for variables in the fluxvar 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;
[fluxvar, fa] = feval(func, physcon, forcvar, surfvar, fluxvar, a);
[fluxvar, fb] = feval(func, physcon, forcvar, surfvar, fluxvar, 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
   [fluxvar, fb] = feval(func, physcon, forcvar, surfvar, fluxvar, 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 [fluxvar] = canopy_turbulence (dt, physcon, forcvar, surfvar, leafvar, soilvar, fluxvar)

% Canopy turbulence, aeorodynamic conductances, and wind/temperature/water vapor
% profiles using above- and within-canopy coupling with a roughness sublayer
% (RSL) parameterization 

% ------------------------------------------------------------------------------------
% Input
%   dt                  ! Model time step (s)
%   physcon.vkc         ! von Karman constant
%   physcon.grav        ! Gravitational acceleration (m/s2)
%   physcon.mmh2o       ! Molecular mass of water (kg/mol)
%   physcon.tfrz        ! Freezing point of water (K)
%   physcon.hvap        ! Latent heat of evaporation (J/kg)
%
%   forcvar.zref        ! Reference height (m)
%   forcvar.uref        ! Wind speed at reference height (m/s)
%   forcvar.thref       ! Potential temperature at reference height (K)
%   forcvar.thvref      ! Virtual potential temperature at reference height (K)
%   forcvar.qref        ! Water vapor at reference height (mol/mol)
%   forcvar.pref        ! Atmospheric pressure (Pa)
%   forcvar.mmair       ! Molecular mass of air at reference height (kg/mol)
%   forcvar.cpair       ! Specific heat of air at constant pressure (J/mol/K)
%   forcvar.rhomol      ! Molar density (mol/m3)
%
%   surfvar.p           ! Index for grid point to process
%   surfvar.hc          ! Canopy height (m)
%   surfvar.pai         ! Canopy plant area index (m2/m2)
%   surfvar.nlev        ! Index for top level
%   surfvar.ntop        ! Index for top leaf layer
%   surfvar.nsoi        ! First canopy layer is soil
%   surfvar.zs          ! Canopy height for scalar concentration and source (m)
%   surfvar.zw          ! Canopy height at layer interfaces (m)
%   surfvar.dpai        ! Layer plant area index (m2/m2)
%   surfvar.fwet        ! Fraction of plant area index that is wet
%   surfvar.fdry        ! Fraction of plant area index that is green and dry
%   surfvar.fracsun     ! Sunlit fraction of canopy layer
%   surfvar.fracsha     ! Shaded fraction of canopy layer
%
%   leafvar.nleaf       ! Number of leaf types (sunlit and shaded)
%   leafvar.isun        ! Sunlit leaf index
%   leafvar.isha        ! Shaded leaf index
%   leafvar.gbh         ! Leaf boundary layer conductance, heat (mol/m2 leaf/s)
%   leafvar.gbv         ! Leaf boundary layer conductance, H2O (mol H2O/m2 leaf/s)
%   leafvar.gs          ! Leaf stomatal conductance (mol H2O/m2 leaf/s)
%   leafvar.cpleaf      ! Leaf heat capacity (J/m2 leaf/K)
%   leafvar.rnleaf      ! Leaf net radiation (W/m2 leaf)
%
%   soilvar.tk          ! Soil thermal conductivity (W/m/K)
%   soilvar.dz          ! Soil layer depth (m)
%   soilvar.tsoi        ! Soil temperature (K)
%   soilvar.resis       ! Soil evaporative resistance (s/m)
%   soilvar.rhg         ! Relative humidity of airspace at soil surface (fraction)
%
%   fluxvar.rnsoi       ! Net radiation at ground (W/m2)
%   fluxvar.tair_old    ! Air temperature profile for previous timestep (K)
%   fluxvar.qair_old    ! Water vapor profile for previous timestep (mol/mol)
%   fluxvar.tveg_old    ! Vegetation temperature profile for previous timestep (K)
%
% Input/output
%   fluxvar.taf         ! Air temperature at canopy top (K)
%   fluxvar.qaf         ! Water vapor at canopy top (mol/mol)
%
% Output (specific to canopy_turbulence)
%   fluxvar.Lc          ! Canopy density length scale (m)
%   fluxvar.wind        ! Wind speed profile (m/s)
%   fluxvar.uaf         ! Wind speed at canopy top (m/s)
%   fluxvar.ga_prof     ! Canopy layer aerodynamic conductance for scalars (mol/m2/s)
%
% Output (from obukhov_function)
%   fluxvar.c1m         ! Roughness sublayer c1 parameter for momentum (dimensionless)
%   fluxvar.c1c         ! Roughness sublayer c1 parameter for scalars (dimensionless)
%   fluxvar.c2          ! Roughness sublayer depth scale multiplier (dimensionless)
%   fluxvar.disp        ! Displacement height (m)
%   fluxvar.beta        ! u* / u(hc)
%   fluxvar.PrSc        ! Prandtl (Schmidt) number at canopy top
%   fluxvar.ustar       ! Friction velocity (m/s)
%   fluxvar.tstar       ! Temperature scale (K)
%   fluxvar.qstar       ! Water vapor scale (mol/mol)
%   fluxvar.gac         ! Aerodynamic conductance for a scalar above canopy (mol/m2/s)
%   fluxvar.obu_ustar   ! Obukhov length used for u* (m)
%   fluxvar.obu         ! Value for Obukhov length (m)
%
% Output (from scalar_profile)
%   fluxvar.tair        ! Air temperature profile (K)
%   fluxvar.qair        ! Water vapor profile (mol/mol)
%   fluxvar.tveg        ! Vegetation temperature profile (K)
%   fluxvar.shsoi       ! Ground sensible heat flux, ground (W/m2)
%   fluxvar.etsoi       ! Ground evaporation flux (mol H2O/m2/s)
%   fluxvar.gsoi        ! Soil heat flux (W/m2)
%   fluxvar.tg          ! Soil surface temperature (K)
%   fluxvar.shveg       ! Leaf sensible heat flux (W/m2 leaf)
%   fluxvar.etveg       ! Leaf evapotranspiration flux (mol H2O/m2 leaf/s)
%   fluxvar.stveg       ! Leaf storage heat flux (W/m2 leaf)
%   fluxvar.shair       ! Canopy air sensible heat flux (W/m2)
%   fluxvar.etair       ! Canopy air water vapor flux (mol H2O/m2/s)
%   fluxvar.stair       ! Canopy air storage heat flux (W/m2)
% ------------------------------------------------------------------------------------

% --- Index for grid point to process

p = surfvar.p;

% --- Leaf drag coefficient (dimensionless)

cd = 0.25;

% --- Canopy density length scale (m)

fluxvar.Lc(p) = surfvar.hc(p) / (cd * surfvar.pai(p));

% --- Calculate the Obukhov length

% Calculate the Obukhov length (obu) for the current surface temperature
% and surface vapor pressure using Harman & Finnigan (2007, 2008) roughness 
% sublayer (RSL) theory. Use the function "obukhov_function" to iterate obu
% until the change in obu is less than tol.

obu0 = 100;                     % Initial estimate for Obukhov length (m)
obu1 = -100;                    % Initial estimate for Obukhov length (m)
tol = 0.01;                     % Accuracy tolerance for Obukhov length (m)
func_name = 'obukhov_function'; % The function is the file obukhov_function.m

% Solve for the Obukhov length. Do not use final returned value for obu.
% Instead, use the value used to calculate u*

[fluxvar, oburoot] = hybrid_root (func_name, physcon, forcvar, surfvar, fluxvar, obu0, obu1, tol);
fluxvar.obu(p) = fluxvar.obu_ustar(p);

% --- Wind profile (m/s)

% Above-canopy wind speed: defined at zs

h_minus_d = surfvar.hc(p) - fluxvar.disp(p);
[psi_m_hc] = psi_m_monin_obukhov (h_minus_d / fluxvar.obu(p));
[psi_m_rsl_hc] = psi_m_rsl (h_minus_d, h_minus_d, fluxvar.obu(p), fluxvar.c1m(p), fluxvar.c2);

for ic = surfvar.ntop(p)+1:surfvar.nlev(p)
   z_minus_d = surfvar.zs(p,ic) - fluxvar.disp(p);
   [psi_m_zs] = psi_m_monin_obukhov (z_minus_d / fluxvar.obu(p));
   [psi_m_rsl_zs] = psi_m_rsl (z_minus_d, h_minus_d, fluxvar.obu(p), fluxvar.c1m(p), fluxvar.c2);
   psim = -psi_m_zs + psi_m_hc + psi_m_rsl_zs - psi_m_rsl_hc + physcon.vkc / fluxvar.beta(p);
   fluxvar.wind(p,ic) = fluxvar.ustar(p) / physcon.vkc * (log(z_minus_d/h_minus_d) + psim);
end

% Wind speed at top of canopy

fluxvar.uaf(p) = fluxvar.ustar(p) / fluxvar.beta(p);

% Within-canopy wind speed: defined at zs. Limit to > 0.1 m/s

lm = 2 * fluxvar.beta(p)^3 * fluxvar.Lc(p);
lm_over_beta = lm / fluxvar.beta(p);
for ic = surfvar.nsoi(p)+1:surfvar.ntop(p)
   fluxvar.wind(p,ic) = fluxvar.uaf(p) * exp((surfvar.zs(p,ic) - surfvar.hc(p)) / lm_over_beta);
   fluxvar.wind(p,ic) = max(fluxvar.wind(p,ic), 0.1);
end

% Wind speed at ground

fluxvar.wind(p,surfvar.nsoi(p)) = 0;

% --- Aerodynamic conductances (mol/m2/s) between zs(i) and zs(i+1)

% Above-canopy conductances

h_minus_d = surfvar.hc(p) - fluxvar.disp(p);
[psi_c_hc] = psi_c_monin_obukhov (h_minus_d / fluxvar.obu(p));
[psi_c_rsl_hc] = psi_c_rsl (h_minus_d, h_minus_d, fluxvar.obu(p), fluxvar.c1c(p), fluxvar.c2);

for ic = surfvar.ntop(p)+1:surfvar.nlev(p)-1

   % Lower height zs(i)

   z_minus_d = surfvar.zs(p,ic) - fluxvar.disp(p);
   [psi_c_zs] = psi_c_monin_obukhov (z_minus_d / fluxvar.obu(p));
   [psi_c_rsl_zs] = psi_c_rsl (z_minus_d, h_minus_d, fluxvar.obu(p), fluxvar.c1c(p), fluxvar.c2);
   psic1 = -psi_c_zs + psi_c_hc + psi_c_rsl_zs - psi_c_rsl_hc;

   % Upper height zs(i+1)

   z_minus_d = surfvar.zs(p,ic+1) - fluxvar.disp(p);
   [psi_c_zs] = psi_c_monin_obukhov (z_minus_d / fluxvar.obu(p));
   [psi_c_rsl_zs] = psi_c_rsl (z_minus_d, h_minus_d, fluxvar.obu(p), fluxvar.c1c(p), fluxvar.c2);
   psic2 = -psi_c_zs + psi_c_hc + psi_c_rsl_zs - psi_c_rsl_hc;

   % Conductance. Note that psic = psic2 - psic1 is equivalent to:
   % -psi_c_z2 + psi_c_z1 + psi_c_rsl_z2 - psi_c_rsl_z1

   psic = psic2 - psic1;
   zlog = log((surfvar.zs(p,ic+1)-fluxvar.disp(p)) / (surfvar.zs(p,ic)-fluxvar.disp(p)));
   fluxvar.ga_prof(p,ic) = forcvar.rhomol(p) * physcon.vkc * fluxvar.ustar(p) / (zlog + psic);
end

% Special case for the top layer to the reference height

ic = surfvar.nlev(p);
z_minus_d = surfvar.zs(p,ic) - fluxvar.disp(p);
[psi_c_zs] = psi_c_monin_obukhov (z_minus_d / fluxvar.obu(p));
[psi_c_rsl_zs] = psi_c_rsl (z_minus_d, h_minus_d, fluxvar.obu(p), fluxvar.c1c(p), fluxvar.c2);
psic1 = -psi_c_zs + psi_c_hc + psi_c_rsl_zs - psi_c_rsl_hc;

z_minus_d = forcvar.zref(p) - fluxvar.disp(p);
[psi_c_zs] = psi_c_monin_obukhov (z_minus_d / fluxvar.obu(p));
[psi_c_rsl_zs] = psi_c_rsl (z_minus_d, h_minus_d, fluxvar.obu(p), fluxvar.c1c(p), fluxvar.c2);
psic2 = -psi_c_zs + psi_c_hc + psi_c_rsl_zs - psi_c_rsl_hc;

psic = psic2 - psic1;
zlog = log((forcvar.zref(p)-fluxvar.disp(p)) / (surfvar.zs(p,ic)-fluxvar.disp(p)));
fluxvar.ga_prof(p,ic) = forcvar.rhomol(p) * physcon.vkc * fluxvar.ustar(p) / (zlog + psic);

% Within-canopy aerodynamic conductances

for ic = surfvar.nsoi(p)+1:surfvar.ntop(p)-1
   zl = surfvar.zs(p,ic) - surfvar.hc(p);
   zu = surfvar.zs(p,ic+1) - surfvar.hc(p);
   res = fluxvar.PrSc(p) / (fluxvar.beta(p) * fluxvar.ustar(p)) * (exp(-zl/lm_over_beta) - exp(-zu/lm_over_beta));
   fluxvar.ga_prof(p,ic) = forcvar.rhomol(p) / res;
end

% Special case for top canopy layer: conductance from zs to hc

ic = surfvar.ntop(p);
zl = surfvar.zs(p,ic) - surfvar.hc(p);
zu = surfvar.hc(p) - surfvar.hc(p);
res = fluxvar.PrSc(p) / (fluxvar.beta(p) * fluxvar.ustar(p)) * (exp(-zl/lm_over_beta) - exp(-zu/lm_over_beta));
ga_below_hc = forcvar.rhomol(p) / res;

% Now include additional conductance from hc to first atmospheric layer

z_minus_d = surfvar.hc(p) - fluxvar.disp(p);
[psi_c_zs] = psi_c_monin_obukhov (z_minus_d / fluxvar.obu(p));
[psi_c_rsl_zs] = psi_c_rsl (z_minus_d, h_minus_d, fluxvar.obu(p), fluxvar.c1c(p), fluxvar.c2);
psic1 = -psi_c_zs + psi_c_hc + psi_c_rsl_zs - psi_c_rsl_hc;

z_minus_d = surfvar.zs(p,ic+1) - fluxvar.disp(p);
[psi_c_zs] = psi_c_monin_obukhov (z_minus_d / fluxvar.obu(p));
[psi_c_rsl_zs] = psi_c_rsl (z_minus_d, h_minus_d, fluxvar.obu(p), fluxvar.c1c(p), fluxvar.c2);
psic2 = -psi_c_zs + psi_c_hc + psi_c_rsl_zs - psi_c_rsl_hc;

psic = psic2 - psic1;
zlog = log((surfvar.zs(p,ic+1)-fluxvar.disp(p)) / (surfvar.hc(p)-fluxvar.disp(p)));
ga_above_hc = forcvar.rhomol(p) * physcon.vkc * fluxvar.ustar(p) / (zlog + psic);
fluxvar.ga_prof(p,ic) = 1 / (1 / ga_below_hc + 1 / ga_above_hc);

% Make sure above-canopy aerodynamic resistances sum to 1/gac

sumres = 1 / ga_above_hc;
for ic = surfvar.ntop(p)+1:surfvar.nlev(p)
   sumres = sumres + 1 / fluxvar.ga_prof(p,ic);
end

if (abs(1/sumres - fluxvar.gac(p)) > 1e-06)
   error('canopy_turbulence: above-canopy aerodynamic conductance error')
end

% Aerodynamic conductance at ground

ic = surfvar.nsoi(p);
ic_plus_one = surfvar.nsoi(p)+1;

z0m_g = 0.01;                         % Roughness length of ground (m)
z0c_g = 0.1 * z0m_g;                  % Roughness length for scalars
zlog_m = log(surfvar.zs(p,ic_plus_one)/z0m_g);
zlog_c = log(surfvar.zs(p,ic_plus_one)/z0c_g);
ustar_g = fluxvar.wind(p,ic_plus_one) * physcon.vkc / zlog_m;
ustar_g = max(ustar_g, 0.01);
res = zlog_c / (physcon.vkc * ustar_g);
fluxvar.ga_prof(p,ic) = forcvar.rhomol(p) / res;

% --- Limit resistances to < 500 s/m

for ic = surfvar.nsoi(p):surfvar.nlev(p)
   res = min (forcvar.rhomol(p)/fluxvar.ga_prof(p,ic), 500);
   fluxvar.ga_prof(p,ic) = forcvar.rhomol(p) / res;
end

% --- Calculate within-canopy scalar profiles for temperature and vapor pressure

[fluxvar] = scalar_profile (dt, p, physcon, forcvar, surfvar, leafvar, soilvar, fluxvar);

% --- Temperature and water vapor at top of canopy

fluxvar.taf(p) = fluxvar.tair(p,surfvar.ntop(p));
fluxvar.qaf(p) = fluxvar.qair(p,surfvar.ntop(p));

function [fluxvar, root] = hybrid_root (func, physcon, forcvar, surfvar, fluxvar, xa, xb, tol)

% Solve for the root of a function using the secant and Brent's methods given
% initial estimates 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 [fluxvar, fx] = func (physcon, forcvar, surfvar, fluxvar, x)
%
% The function func is exaluated at x and the returned value is fx. It uses variables
% in the physcon, forcvar, surfvar, and fluxvar structures. These are passed in as
% input arguments. It also calculates values for variables in the fluxvar 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 see if this is the root

x0 = xa;
[fluxvar, f0] = feval(func, physcon, forcvar, surfvar, fluxvar, x0);
if (f0 == 0)
   root = x0;
   return
end

% --- Evaluate func at xb and see if this is the root

x1 = xb;
[fluxvar, f1] = feval(func, physcon, forcvar, surfvar, fluxvar, x1);
if (f1 == 0)
   root = x1;
   return
end

% --- Order initial root estimates correctly

if (f1 < f0)
   minx = x1;
   minf = f1;
else
   minx = x0;
   minf = f0;
end

% --- Iterative root calculation. Use the secant method, with Brent's method as a backup

itmax = 40;
for iter = 1:itmax
   dx = -f1 * (x1 - x0) / (f1 - f0);
   x = x1 + dx;

   % Check if x is the root. If so, exit the iteration

   if (abs(dx) < tol)
      x0 = x;
      break
   end

   % Evaluate the function at x

   x0 = x1;
   f0 = f1;
   x1 = x;
   [fluxvar, f1] = feval(func, physcon, forcvar, surfvar, fluxvar, x1);
   if (f1 < minf)
      minx = x1;
      minf = f1;
   end

   % If a root zone is found, use Brent's method for a robust backup strategy
   % and exit the iteration

   if (f1 * f0 < 0)
      [fluxvar, x] = brent_root (func, physcon, forcvar, surfvar, fluxvar, x0, x1, tol);
      x0 = x;
      break
   end

   % In case of failing to converge within itmax iterations stop at the minimum function

   if (iter == itmax)
      [fluxvar, f1] = feval(func, physcon, forcvar, surfvar, fluxvar, minx);
      x0 = minx;
   end

end

root = x0;

function [fluxvar, fx] = obukhov_function (physcon, forcvar, surfvar, fluxvar, x)

% Use Harman & Finnigan (2007, 2008) roughness sublayer (RSL) theory to obtain
% the Obukhov length (obu). This is the function to solve for the Obukhov
% length. For the current estimate of the Obukhov length (x), calculate
% u*, T*, and q* and then the new length (obu). The function value is the
% change in Obukhov length: fx = x - obu.

% -------------------------------------------------------------------------
% Input
%   x                   ! Current estimate for Obukhov length (m)
%   physcon.vkc         ! von Karman constant
%   physcon.grav        ! Gravitational acceleration (m/s2)
%   physcon.mmh2o       ! Molecular mass of water (kg/mol)
%   forcvar.zref        ! Reference height (m)
%   forcvar.uref        ! Wind speed at reference height (m/s)
%   forcvar.thref       ! Potential temperature at reference height (K)
%   forcvar.thvref      ! Virtual potential temperature at reference height (K)
%   forcvar.qref        ! Water vapor at reference height (mol/mol)
%   forcvar.rhomol      ! Molar density (mol/m3)
%   forcvar.mmair       ! Molecular mass of air at reference height (kg/mol)
%   fluxvar.taf         ! Air temperature at canopy top (K)
%   fluxvar.qaf         ! Water vapor at canopy top (mol/mol)
%   fluxvar.Lc          ! Canopy density length scale (m)
%   surfvar.hc          ! Canopy height (m)
%   surfvar.p           ! Index for grid point to process
%
% Output
%   fluxvar.c1m         ! Roughness sublayer c1 parameter for momentum (dimensionless)
%   fluxvar.c1c         ! Roughness sublayer c1 parameter for scalars (dimensionless)
%   fluxvar.c2          ! Roughness sublayer depth scale multiplier (dimensionless)
%   fluxvar.disp        ! Displacement height (m)
%   fluxvar.beta        ! u* / u(hc)
%   fluxvar.PrSc        ! Prandtl (Schmidt) number at canopy top
%   fluxvar.ustar       ! Friction velocity (m/s)
%   fluxvar.tstar       ! Temperature scale (K)
%   fluxvar.qstar       ! Water vapor scale (mol/mol)
%   fluxvar.gac         ! Aerodynamic conductance for a scalar above canopy (mol/m2/s)
%   fluxvar.obu_ustar   ! Obukhov length used for u* (m)
%   fluxvar.obu         ! Value for Obukhov length (m)
%   fx                  ! Change in Obukhov length (x - obu)
% -------------------------------------------------------------------------

% --- Index for grid point to process

p = surfvar.p;

% --- Prevent near-zero values of Obukhov length

if (abs(x) <= 0.1)
   x = 0.1;
end

% --- Determine beta_val = u* / u(hc) for the current Obukhov length

% Neutral value for beta = u* / u(hc)

beta_neutral = 0.35;

% Lc/obu

LcL = fluxvar.Lc(p)/x;

if (LcL <= 0)

   % The unstable case is a quadratic equation for beta^2 at LcL

   a = 1;
   b = 16 * LcL * beta_neutral^4;
   c = -beta_neutral^4;
   beta_val = sqrt((-b + sqrt(b^2 - 4 * a * c)) / (2 * a));

   % Error check

   y = beta_val^2 * LcL;
   fy = (1 - 16 * y)^(-0.25);
   err = beta_val * fy - beta_neutral;
   if (abs(err) > 1e-10)
      error('obukhov_function: unstable case - error in beta')
   end

else

   % The stable case is a cubic equation for beta at LcL

   a = 5 * LcL;
   b = 0;
   c = 1;
   d = -beta_neutral;
   q = (2*b^3 - 9*a*b*c + 27*(a^2)*d)^2 - 4*(b^2 - 3*a*c)^3;
   q = sqrt(q);
   r = 0.5 * (q + 2*b^3 - 9*a*b*c + 27*(a^2)*d);
   r = r^(1/3);
   beta_val = -(b+r)/(3*a) - (b^2 - 3*a*c)/(3*a*r);

   % Error check

   y = beta_val^2 * LcL;
   fy = 1 + 5 * y;
   err = beta_val * fy - beta_neutral;
   if (abs(err) > 1e-10)
      error('obukhov_function: stable case - error in beta')
   end

end

% Place limits on beta

beta_val = min(0.5, max(beta_val,0.2));
fluxvar.beta(p) = beta_val;

% --- For current beta = u*/u(hc) determine displacement height

dp = beta_val^2 * fluxvar.Lc(p);               % dp = hc - disp
fluxvar.disp(p) = max(surfvar.hc(p) - dp, 0);  % Displacement height (m)

% Save reference height and canopy height (relative to displacement height)
% because these are used many times

z_minus_d = forcvar.zref(p) - fluxvar.disp(p);
h_minus_d = surfvar.hc(p) - fluxvar.disp(p);

% --- Turbulent Prandlt (Schmidt) number (PrSc) at canopy height

PrSc0 = 0.5;        % Neutral value for Pr (Sc)
PrSc1 = 0.3;        % Magnitude of variation of Pr (Sc) with stability
PrSc2 = 2.0;        % Scale of variation of Pr (Sc) with stability

fluxvar.PrSc(p) = PrSc0 + PrSc1 * tanh(PrSc2*fluxvar.Lc(p)/x);

% --- Calculate the parameters c1 and c2 needed for the RSL function phi_hat

% Evaluate Monin-Obukhov phi functions at (hc-disp)/obu

[phi_m_hc] = phi_m_monin_obukhov (h_minus_d / x);
[phi_c_hc] = phi_c_monin_obukhov (h_minus_d / x);

% Roughness sublayer depth scale multiplier (dimensionless)

fluxvar.c2 = 0.5;

% c1 for momentum and scalars (dimensionless)

fluxvar.c1m(p) = (1 -                 physcon.vkc / (2 * beta_val * phi_m_hc)) * exp(fluxvar.c2/2);
fluxvar.c1c(p) = (1 - fluxvar.PrSc(p)*physcon.vkc / (2 * beta_val * phi_c_hc)) * exp(fluxvar.c2/2);

% --- Evaluate the roughness sublayer psi_hat functions for momentum and scalars

% These are calculated at the reference height and at the canopy height. Note that
% here the heights are adjusted for the displacement height before the integration.

[psi_m_rsl_zref] = psi_m_rsl (z_minus_d, h_minus_d, x, fluxvar.c1m(p), fluxvar.c2);  % momentum at (zref-disp)
[psi_m_rsl_hc]   = psi_m_rsl (h_minus_d, h_minus_d, x, fluxvar.c1m(p), fluxvar.c2);  % momentum at (hc-disp)

[psi_c_rsl_zref] = psi_c_rsl (z_minus_d, h_minus_d, x, fluxvar.c1c(p), fluxvar.c2);  % scalars at (zref-disp)
[psi_c_rsl_hc]   = psi_c_rsl (h_minus_d, h_minus_d, x, fluxvar.c1c(p), fluxvar.c2);  % scalars at (hc-disp)

% --- Evaluate the Monin-Obukhov psi functions for momentum and scalars

% These are calculated at the reference height and at the canopy height

[psi_m_zref] = psi_m_monin_obukhov (z_minus_d / x);    % momentum at (zref-disp)/obu
[psi_m_hc]   = psi_m_monin_obukhov (h_minus_d / x);    % momentum at (hc-disp)/obu

[psi_c_zref] = psi_c_monin_obukhov (z_minus_d / x);    % scalars at (zref-disp)/obu
[psi_c_hc]   = psi_c_monin_obukhov (h_minus_d / x);    % scalars at (hc-disp)/obu

% --- Calculate u* (m/s), T* (K), q* (mol/mol), and Tv* (K)

zlog = log(z_minus_d / h_minus_d);
psim = -psi_m_zref + psi_m_hc + psi_m_rsl_zref - psi_m_rsl_hc + physcon.vkc / beta_val;
psic = -psi_c_zref + psi_c_hc + psi_c_rsl_zref - psi_c_rsl_hc;

fluxvar.ustar(p) = forcvar.uref(p) * physcon.vkc / (zlog + psim);
fluxvar.tstar(p) = (forcvar.thref(p) - fluxvar.taf(p)) * physcon.vkc / (zlog + psic);
fluxvar.qstar(p) = (forcvar.qref(p) - fluxvar.qaf(p)) * physcon.vkc / (zlog + psic);
fluxvar.obu_ustar(p) = x;

% --- Aerodynamic conductance to canopy height

fluxvar.gac(p) = forcvar.rhomol(p) * physcon.vkc * fluxvar.ustar(p) / (zlog + psic);

% --- Calculate new Obukhov length (m)

tvstar = fluxvar.tstar(p) + 0.61 * forcvar.thref(p) * fluxvar.qstar(p) * (physcon.mmh2o / forcvar.mmair(p));
fluxvar.obu(p) = fluxvar.ustar(p)^2 * forcvar.thvref(p) / (physcon.vkc * physcon.grav * tvstar);
fx = x - fluxvar.obu(p);